Relativistic second-order dissipative hydrodynamics at finite chemical potential

Starting from the Boltzmann equation in the relaxation time approximation and employing a Chapman-Enskog like expansion for the distribution function close to equilibrium, we derive second-order evolution equations for the shear stress tensor and the dissipative charge current for a system of massless quarks and gluons. The transport coefficients are obtained exactly using quantum statistics for the phase space distribution functions at non-zero chemical potential. We show that, within the relaxation time approximation, the second-order evolution equations for the shear stress tensor and the dissipative charge current can be decoupled. We find that, for large values of the ratio of chemical potential to temperature, the charge conductivity is small compared to the coefficient of shear viscosity. Moreover, we show that in the relaxation-time approximation, the limiting behaviour of the ratio of heat conductivity to shear viscosity is qualitatively similar to that obtained for a strongly coupled conformal plasma.


I. INTRODUCTION
High-energy heavy ion collisions at the BNL Relativistic Heavy Ion Collider (RHIC) [1,2] and the CERN Large Hadron Collider (LHC) [3][4][5] create strongly interacting matter under extreme conditions of high temperature and density as it is believed to have existed in the very early universe [6,7]. At such conditions, quarks and gluons are deconfined to form a new state of matter, the quark-gluon plasma (QGP). The QGP behaves as a strongly coupled plasma having the smallest shear viscosity-to-entropy density ratio, η/s [8][9][10][11][12][13]. Relativistic hydrodynamics has been applied quite successfully to describe the space-time evolution of the QGP formed in high-energy heavy ion collisions and estimate its transport coefficients [14].
In applications of hydrodynamics it is rather straightforward to employ the ideal (Euler) equations. The inclusion of dissipative effects in the evolution of the QGP started only a few years ago. However, most of the studies have focused on exploring the effects of the shear viscosity on the QGP evolution and extracting its magnitude from experimental measurements. Nevertheless, there are other sources of dissipation such as bulk viscous pressure and dissipative charge current that may have a significant effect on the hydrodynamic evolution of the QGP. While the effects of bulk viscous pressure has been studied in some details [15][16][17][18][19][20], the dissipative charge current has been largely ignored. This may be attributed to the fact that at very high energies, baryon number and its corresponding chemical potential are negligible. However, at lower collision energies such as those probed in the RHIC low-energy scan or at the upcoming experiments at the Facility for Antiproton and Ion Research (FAIR), baryon number can no longer be ignored and therefore charge diffusion may play an important role.
The earliest theoretical formulation of relativistic dissipative hydrodynamics are due to Eckart [21] and Landau-Lifshitz [22]. However these formulations, collectively called relativistic Navier-Stokes theory, involve only firstorder gradients and suffer from acausality and numerical instability due to the parabolic nature of the equations. Second order or extended theories by Grad [23], Müller [24] and Israel and Stewart (IS) [25] were introduced to restore causality. Therefore it is imperative that second order dissipative hydrodynamic equations should be employed in order to correctly describe the evolution of the QGP. However, the IS formulation of a causal theory of relativistic hydrodynamics from kinetic theory, contains several inconsistencies and approximations, the resolution of which is currently an active research area [26][27][28][29][30][31][32][33][34][35][36][37][38].
In order to formulate a causal theory of relativistic dissipative hydrodynamics from kinetic theory, it is desirable to first specify the form of the non-equilibrium phase-space distribution function. For a system close to local thermodynamic equilibrium, the non-equilibrium corrections to the distribution function can be obtained using either (i) Grad's moment method [23] or (ii) the Chapman-Enskog method [39]. Although both methods involve expanding the distribution function around its equilibrium value, it has been demonstrated that the Chapman-Enskog method in the relaxation time approximation results in a better agreement with microscopic Boltzmann simulations [32,33] as well as with exact solutions of the Boltzmann equation in the relaxation-time approximation [32][33][34][35][36].
In the absence of conserved charges, the Chapman-Enskog method has been used to compute the secondorder transport coefficients for vanishing [32][33][34] as well as finite particle masses [35,36]. On the other hand, in the presence of conserved charges but for classical particles with vanishing masses, the second-order transport coefficients corresponding to charge diffusion (or alternatively heat conduction) have been obtained by employing the moment method [40,41]. However, they still remain to be determined for quantum statistics. Here, we employ the Chapman-Enskog method to achieve this.
In this Letter, we present the derivation of secondorder evolution equations for shear stress tensor and dissipative charge current for a system consisting of massless quarks and gluons. In order to obtain the form of the non-equilibrium distribution function, we employ a Chapman-Enskog like expansion to iteratively solve the Boltzmann equation in the relaxation time approximation [32]. Using this expansion, we derive the first-order constitutive relations and subsequently the second-order evolution equations for the dissipative quantities. The transport coefficients are obtained exactly using quantum statistics for the quark and gluon phase-space distribution functions with a non-vanishing quark chemical potential. Moreover, we show that, up to second-order in the gradient expansion, the evolution equations for the shear stress tensor and the dissipative charge current can be decoupled. We also find that, for large values of the ratio of chemical potential to temperature, the charge conductivity is small compared to the coefficient of shear viscosity. Finally we demonstrate that the limiting behaviour of the heat conductivity to shear viscosity ratio, obtained here in the relaxation-time approximation, is qualitatively identical to that of a conformal fluid in the strong coupling limit.

II. RELATIVISTIC HYDRODYNAMICS
In the case of massless partons, i.e., massless quarks and gluons, the conserved energy-momentum tensor and the net-quark current can be expressed in terms of the single particle phase-space distribution function as [42] T µν = dp p µ p ν [g q (f q + fq) where dp = dp/[(2π) 3 |p|], p µ is the particle four momenta, and g q and g g are the quark and gluon degeneracy factor, respectively. Here f q , fq, and f g are the phasespace distribution functions for quarks, anti-quarks, and gluons. In the tensor decompositions, , P , and n are the energy density, pressure, and the net quark number density. The projection operator ∆ µν = g µν − u µ u ν is orthogonal to the hydrodynamic four-velocity u µ defined in the Landau frame: T µν u ν = u µ . We work with the Minkowskian metric tensor g µν ≡ diag(+, −, −, −).
The dissipative quantities in Eqs. (1) and (2) are the shear stress tensor π µν and the particle diffusion current n µ . With the definition of the energy-momentum tensor in Eq. (1), the bulk viscous pressure vanishes in the massless case. The energy-momentum conservation, ∂ µ T µν = 0, and particle four-current conservation, ∂ µ N µ = 0, yields the fundamental evolution equations for , u µ and n, as˙ Here we use the standard notationȦ = u µ ∂ µ A for comoving derivatives, θ ≡ ∂ µ u µ for the expansion scalar, θ∆ µν for the velocity stress tensor, and ∇ α = ∆ µα ∂ µ for space-like derivatives.
In the following, we briefly outline the thermodynamic properties of a QGP in equilibrium. In this case, the phase-space distribution functions for quarks, antiquarks and gluons are given by respectively, where u · p ≡ u µ p µ , β = 1/T is the inverse temperature and α = µ/T is the ratio of the quark chemical potential to temperature. We consider vanishing chemical potential for gluons because they are unconstrained by the conservation laws. The temperature, T , and chemical potential, µ, of the system is determined by the matching condition = 0 and n = n 0 , where 0 and n 0 is the energy density and the net quark number density in equilibrium. The energy density, pressure and the net quark number density for a system of massless quarks and gluons in equilibrium is given by The equilibrium entropy density then becomes The above expressions for 0 , P 0 , n 0 , and s 0 can also be obtained directly from the partition function of an ideal QGP [42], where V is the volume of the system. Indeed, using the thermodynamic relations one recovers Eqs. (9)- (12). The matching conditions = 0 and n = n 0 allows us to define thermodynamic quantities like temperature and chemical potential of a dissipative system. The pressure P can then be obtained from the equation of state of the system.
Even if the equation of state, relating , P and n is provided, Eqs. (3)-(5) are not closed unless the dissipative quantities π µν and n µ are specified. However, before we derive the evolution equations for the dissipative quantities, we need to obtain expressions for the derivatives of β and α. Using Eqs. (3)-(5) and Eqs. (9)-(11), we geṫ where O(δ 2 ) represents terms which are of second or higher order in derivatives. Since dissipative forces are caused by thermodynamic gradients present in a nonideal system, π µν and n µ are at least linear in the gradient expansion. Note that while Eq. (15) is terminated at first-order (sufficient for the present work), Eq. (16) is exact.
The QGP is a strongly coupled system and is conjectured to be formed close to local thermodynamic equilibrium. Therefore, the phase-space distribution function can be split into equilibrium and non-equilibrium parts, 1. Hence, from Eqs. (1) and (2), the shear stress tensor π µν and the particle diffusion current n µ can be expressed in terms of δf as π µν = ∆ µν αβ dp p α p β [g q (δf q + δfq) + g g δf g ] , is a traceless symmetric projection operator orthogonal to u µ and ∆ µν . In the following, we obtain δf up to first order by using the iterative solution of the Boltzmann equation in the relaxation-time approximation and then derive secondorder evolution equations for the dissipative quantities.

III. DISSIPATIVE EVOLUTION EQUATIONS
The determination of the form of the non-equilibrium phase-space distribution function is a central problem in statistical physics. This can be achieved by solving a kinetic equation like the Boltzmann equation. The Boltzmann equation governs the evolution of the distribution function which provides a complete description of the microscopic dynamics of a system in the dilute limit. The relativistic Boltzmann equation with the collision term written in the relaxation-time approximation is [43], where τ R is the relaxation time. Note, that for different species of particles, with inter-and intra-species interactions, the relaxation-times are usually distinct. Thus, in general, one should consider the QGP as a true multicomponent system. In the following, we consider the special case with a common relaxation time for all particle species in the QGP. The general case with a different relaxation time is left for future work.
We employ iterative solution of the Boltzmann equation (19) to derive the dissipative equations [31][32][33]. The first-order expressions for shear stress tensor and dissipative charge current are obtained as [32], Here β π and β n are the first-order transport coefficients obtained after performing the momentum integrations in Eqs. (17) and (18). For a system of massless partons, as considered here, we find where, Heref (0) ≡ 1 − rf (0) , where r = 1 for Fermions (quarks and anti-quarks) and r = −1 for Bosons (gluons). Using Eq. (20), we also obtain the first-order dissipative corrections to the distribution function, Note, that the tensorial form of dissipative corrections to the distribution function, as given in the above equations, is analogous to that of Grad's 14-moment approximation [23]. However, the coefficients of these terms are different which has interesting implications in the context of relativistic heavy-ion collisions [38].
To obtain the second-order evolution equations, we follow the procedure discussed in Ref. [26]. We consider the comoving derivative of Eqs. (17) and (18), and rewrite Eq. (19) in favour of δḟ . Using Eqs. (23)- (25) and performing the momentum integrations, we finally obtain the second-order evolution equation for π µν and n µ , Here ω µν ≡ (∇ µ u ν − ∇ ν u µ )/2 is the anti-symmetric vorticity tensor and we have ignored terms higher than quadratic order in the gradients [32]. Note that in the relaxation-time approximation, the Boltzmann relaxation time τ R is the time scale for evolution of both π µν and n µ , i.e., τ π = τ n = τ R . By comparing the first-order equations, Eq. (20), with the relativistic Navier-Stokes equations for the dissipative quantities [22], the dissipative relaxation times can be related to the firstorder transport coefficients τ π = η/β π and τ n = κ n /β n . It is interesting to note that terms like ∆ µ ν ∂ γ π νγ and π µνu ν do not appear in Eq. (27) because their coefficients vanish. Note also that, with β π = ( + P )/5 calculated using the appropriate equation of state, Eq. (26) is valid even in a hadronic phase dominated by massless pions.
The last term in Eq. (27) couples the evolution of dissipative charge current with the shear stress tensor. This type of coupling leads to disagreement with transport results as shown in Ref. [40]. We observe, however, that using Eq. (20), the last term in Eq. (27) is, up to second-order in the gradient expansion, equivalent to −(6/5)n ν σ νµ . Thus, the evolution equation, Eq. (27), for the dissipative charge current becomeṡ This is the main result of the present work. The compact form of the above equation makes it straightforward for direct implementation in a viscous hydrodynamic code. Note that while it is possible to formally rewrite any second-order equation only in terms of gradients using the first-order expressions, Eq. (20), it does not usually imply decoupling of the dissipative evolution equations [44]. However, we have ensured that the second-order terms in Eq. (29) is product of a dissipative quantity and a gradient, as inherent in Eq. (26). It is important to note that the gradient expansion converges only for small deviations from equilibrium. This implies that the second-order scheme is, strictly speaking, justified only if the deviations from the first-order  20), are of second order, or smaller. This is true, in general, to ensure the convergence of gradient expansion in the formulation of dissipative hydrodynamics. Hence, for consistency, the initial conditions for the second-order evolution equations, Eqs. (26) and (29), should be chosen such that the constitutive relations, Eq. (20), are satisfied. Nevertheless, it was shown in Refs. [33][34][35][36] that the solutions of the second-order dissipative hydrodynamic equations rapidly converge to the exact solution of the Boltzmann equation, irrespective of the choices of initial conditions for the dissipative quantities. Thus, the use of Eq. (20), to modify the secondorder terms in the evolution equation for the dissipative charge current, is tenable even if the initial conditions are chosen such that these relations are violated initially.

IV. TRANSPORT COEFFICIENTS
The dimensionless ratio κ n T /η is a measure of the relative importance of the charge conductivity and the shear viscosity. This quantity, in the relaxation-time approximation, is given by κ n T /η = β n T /β π , which can be studied as a function µ/T . To quantify this ratio, one still need to specify the appropriate quark and gluon degeneracy factors, g q and g g , as where N s = 2 is the number of spin degrees of freedom, N c = 3 is the number of colours, and N f is the number of flavours. In Fig. 1, we show the ratio κ n T /η as a function of µ/T for N f = 2 and N f = 3. We observe that, while for small µ/T , this ratio is almost constant, it drops rapidly for larger µ/T , indicating that the conductivity of the QGP is small relative to the shear viscosity at low temperature and high density. Although the qualitative behaviour remains the same for N f = 2 and N f = 3, the drop in κ n T /η is more pronounced for N f = 3. Moreover, we note that at µ/T > 1 the ratio κ n T /η is almost N f independent.
A further interesting quantity is the ratio of thermal conductivity to shear viscosity. The heat flow is related to the dissipative charge current via the relation q µ = −( + P )n µ /n, and is given as [45] where κ q is the coefficient of thermal conductivity. Using the first-order relations in the relaxation-time approximation, we find In Fig. 2, we plot κ q /η scaled by the factor µ 2 /π 2 T versus µ/T , for a two and three flavor QGP. We observe a constant behaviour in the limit of small as well as large µ/T . Moreover, for large µ/T , we see that (κ q /η)µ 2 /π 2 T is independent of number of flavors. These limiting behaviours have interesting consequences.
In the limit of both small and large µ/T , Eq. (32) reduces to which is similar to the Wiedemann-Franz law [45,46]. The factor π 2 in the above equation is due to quantum statistics and it does not appear for a classical Boltzmann gas. In the limit of small µ/T , the constant C in Eq. (33) is C = (4g g + 7g q )/9g q . Thus, C = 37/27 and 95/81 for two and three flavor QGP, respectively; see Fig. 2.
On the other hand, for large µ/T we find C = 5/3, independent of the number of flavors, as shown in Fig. 2.
These values of C are comparable to C = 8/9 (here a factor of 1/9 indicates that the baryon chemical potential is three times the quark chemical potential) obtained in the calculations for strongly coupled conformal plasmas with finite chemical potential [45]. However, it should be noted that for a strongly coupled conformal plasma, the coefficient C depends on the number of space-time dimensions and is shown to be equal to 32/9, 8/9 and 2/9 for four, five and seven dimensions, respectively [47]. Nevertheless, it is intriguing that up to a constant of proportionality, the limiting behaviour of the ratio κ q /η, obtained in the relaxation-time approximation, is exactly the same as that derived in the case of a strongly coupled conformal fluid.

V. CONCLUSIONS AND OUTLOOK
In this paper we employed the iterative Chapman-Enskog method to derive the second-order dissipative hydrodynamical equations for a system of massless quarks and gluons. The bulk viscous pressure vanishes for such a system and therefore the dissipation is solely due to the shear stress tensor and the dissipative charge current. For the equilibrium distribution function, we considered quantum statistics with non vanishing quark chemical potential. We obtained novel, exact relations for the second-order transport coefficients corresponding to the dissipative charge current evolution.
Moreover, we demonstrated that the evolution equations for shear stress tensor and dissipative charge current can be decoupled. We also found that, for large values of the ratio of chemical potential to temperature, the charge conductivity is small relative to the shear viscosity. Finally, we showed that in the relaxation-time approximation, the limiting behaviour of the ratio of heat conductivity to shear viscosity is qualitatively similar to that of a conformal fluid in the strong coupling regime.
At this juncture, we would like to stress that the iterative Chapman-Enskog approach employed here to obtain the dissipative evolution equations from the Boltzmann equation in the relaxation-time approximation is compatible with the gradient expansion inherent in the formulation of dissipative hydrodynamics, as opposed to the moment method [32]. Looking forward, it would be interesting to determine the effect of the dissipative charge current in high-energy heavy-ion collisions, by implementing the dissipative equations derived here, in realistic hydrodynamic simulations. A further challenging problem would be to extend the current second-order formulation to third order. We leave these questions for future studies.