Flow harmonics v_{n} at finite density

# Flow harmonics vn at finite density

Yoshitaka Hatta Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan    Akihiko Monnai RIKEN BNL Research Center, Brookhaven National Laboratory, Upton 11973 NY, USA    Bo-Wen Xiao Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics, Central China Normal University, Wuhan 430079, China
July 13, 2019
###### Abstract

We investigate the Gubser solution of viscous hydrodynamics at finite density and analytically compute the flow harmonics . We explicitly show how and their viscous corrections depend on the chemical potential. The difference in between particles and antiparticles is also analytically computed and shown to be proportional to various chemical potentials and the viscosity. Excellent agreement is obtained between the results and the available experimental data from the SPS, RHIC and the LHC.

###### pacs:
47.75.+f, 12.38.Mh, 11.25.Hf
preprint: YITP-15-39

## I Introduction

Relativistic hydrodynamics is a general theoretical framework to describe the collective dynamics of high-energy systems near local thermal equilibrium. Its first application to hadron physics dates back to Landau’s attempt to describe multi-particle production in hadron-hadron collisions Landau (). It has become a topic of great interest since the discovery of the quark-gluon plasma (QGP) as a nearly-perfect fluid in the “Little Bangs” at BNL Relativistic Heavy Ion Collider (RHIC) Adcox:2004mh (); Adams:2005dq (); Back:2004je (); Arsene:2004fa () and CERN Large Hadron collider (LHC) Aamodt:2010pa (); ATLAS:2011ah (); Chatrchyan:2012wg (). This is supported by the observations that the azimuthal momentum anisotropy of hadronic distribution Ollitrault:1992bk (); Poskanzer:1998yz (), characterized by flow harmonics , are found to reflect the geometrical anisotropy of the overlapping region of two colliding nuclei, and that they are in good quantitative agreement with theoretical estimations. Nowadays the viscous hydrodynamic modeling is considered as one of the most powerful tools to quantify the QGP medium near the crossover phase transition Schenke:2010rr ().

The recent Beam Energy Scan (BES) experiments at RHIC pose us intriguing challenges to study the properties of the medium at finite density and to explore the QCD phase diagram to find signs of a critical point star (). Conserved charges such as net baryon number, strangeness and isospin would play important roles in the collisions with lower energies, as the differences between particle and antiparticle yields are clearly seen Adamczyk:2013gv (); Adamczyk:2013gw (). Historically, it had long been speculated based on several idealized calculations that the strong coupling limit is achieved only at highest energies of RHIC experiments. On the other hand, recent improvements in off-equilibrium hydrodynamic modeling motivates us to reexamine the validity of hydrodynamics in exploring the dense quark matter created at mid-low energies, especially since the differential elliptic flow is found to remain large in Phase I of the BES experiments. The applicability of hydrodynamic models is closely related to the origin of fluidity, about which little is known, and thus its verification would be a very important step towards a full understanding of the hot QCD medium.

So far many hydrodynamic analyses have been performed numerically because it is generally quite nontrivial to solve the partial differential equations involved. Analytical solutions of relativistic hydrodynamics, on the other hand, can be obtained with certain symmetry conditions and they are very instructive in understanding the essence of heavy-ion dynamics. The boost-invariant Bjorken flow Bjorken:1982qr () is one such classic example. More recently, Gubser found an exact boost-invariant solution of the Navier-Stokes equation which has a nontrivial dependence on the transverse coordinate Gubser:2010ze (). The latter solution has the advantage that one can add azimuthally anisotropic perturbations Gubser:2010ui (); Hatta:2014upa () and analytically compute the corresponding flow harmonics including the viscosity effects Hatta:2014upa (); Hatta:2014jva () (see, also, Csanad:2008af (); Peschanski:2009tg (); Csanad:2014dpa ()).

In this study, we investigate at finite density by analytically solving the viscous hydrodynamic equations coupled with conserved currents assuming conformal and boost-invariant symmetries. Aside from the fact that the solution itself is new and of theoretical importance, it gives us a theoretical guidance about the behavior of over a wide range of the beam energy for which there are not many numerical simulations Kolb:2000sd (); Teaney:2001av (); Petersen:2009vx (); Shen:2012vn (); Karpenko:2015xea () and the previous knowledge obtained through the precision analyses in the RHIC-LHC energy regime, such as the value of the shear viscosity , are no longer fully applicable. We discuss extensively the nature of flow in the presence of currents and estimate the beam energy (or chemical potential) dependence of . The difference in between particles and antiparticles is also analytically computed. The results are compared with the experimental data from SPS, RHIC and the LHC Alt:2003ab (); Adamczyk:2013gv (); Adamczyk:2013gw (); CMS:2013bza (). We see that they are in qualitative agreement, which suggests that a reasonable description of the low-energy experimental data might be possible within a hydrodynamic framework.

The paper is organized as follows. The basic setup of relativistic hydrodynamics is outlined in Section II. We then present analytical formulas of the flow harmonics in the ideal and viscous cases in Sections III and IV, respectively. Phenomenological inputs for our model are summarized in Section V. Using these formulas and input parameters, we compare our results with the experimental data in Section VI. Section VII is devoted to summary and conclusions.

## Ii Hydrodynamic equations

### ii.1 Setup

We shall consider hydrodynamics of a conformal theory. The system is characterized by the local temperature and a set of local chemical potentials where the subscript labels various conserved charges of the theory. The flow velocity is denoted by with the normalization . The energy-momentum tensor in the Navier-Stokes approximation takes the form

 Tμν=4ε3uμuν+ε3gμν−2ησμν, (1)

where is the shear tensor and is the shear viscosity. In (1), the conformal equation of state between the energy density and the pressure has been used. The conserved current associated with the chemical potential can be written as

 Jμi=niuμ−κi(uμuν+gμν)∂ν(μiT), (2)

where is the charge density and is the charge conductivity. The hydrodynamic equations consist of the conservation equations for and

 ∇μTμν = 0,∇μJμi=0, (3)

where is the covariant derivative.

Since there is no intrinsic mass scale in a conformal theory, the energy density and the charge densities can be generically written as

 ε=T4f(μ1T,μ2T,⋯),ni=μiT2gi(μ1T,μ2T,⋯). (4)

With a view to applying to heavy-ion collisions, we shall focus on the following representative situation. We assume that there is the leading current (‘baryon number current’) and the corresponding chemical potential is treated to all orders. In addition, there is one subleading current (‘isospin number current’) whose chemical potential is small and treated only to linear order. We take to be ‘orthogonal’ to , in that is invariant under a sign flip (i.e., cross terms like are absent). With these assumptions, we can parameterize

 ε=T4f(μT),n=μT2g(μT),~n=~μT2~g(μT). (5)

The last equation may be written as where is the susceptibility.

### ii.2 Gubser flow

We shall solve the hydrodynamic equations (3) for a given flow velocity

 uτ=cosh[tanh−12τx⊥L2+τ2+x2⊥],u⊥=sinh[tanh−12τx⊥L2+τ2+x2⊥], (6)

and . The parameter is the characteristic length scale of the system. In heavy-ion collisions, it is roughly the transverse size of the colliding nuclei. Eq. (6) is called Gubser flow Gubser:2010ze (); Gubser:2010ui () expressed in the coordinate system

 ds2=−dτ2+τ2dζ2+dx2⊥+x2⊥dϕ2, (7)

where is the proper time, is the spacetime rapidity and is the transverse coordinate. The condition means that the flow is boost invariant along the beam () direction.

Gubser flow takes a very simple form in a cleverly chosen coordinate system which is related to the Minkowski coordinates via a Weyl rescaling of the metric.

 d^s2=ds2τ2=−dρ2+cosh2ρ(dΘ2+sin2Θdϕ2)+dζ2, (8)

where

 sinhρ=−L2−τ2+x2⊥2Lτ,tanΘ=2Lx⊥L2+τ2−x2⊥. (9)

In this coordinate system, the flow velocity is simply . In addition to the boost invariance, the flow respects the symmetry with respect to the ‘polar’ angles . Variables in this coordinate system will be denoted with a ‘hat’, e.g., , .

## Iii Inviscid case

In this section, we solve the hydrodynamic equations (3) in the ideal case . We then deform the solution in the azimuthal direction and compute flow harmonics .

### iii.1 Isotropic ideal solution

The isotropic solution (i.e., independent of ) has been obtained already in Gubser:2010ze (); Gubser:2010ui () in the presence of a current . Assuming that all the quantities depend only on , we can readily solve the hydrodynamic equations for and in the coordinates (8). We then perform the Weyl transformation back to the Minkowski space , to get

 ε0=T40f(μ0T0) ∝ 1τ4(coshρ)8/3, (10)
 n0=μ0T20g(μ0T0)∝1τ3cosh2ρ, ~n0=~μ0T20~g(μ0T0)∝1τ3cosh2ρ. (11)

These equations can be solved for and . It is consistent to look for the solution where and have the same -dependence such that the ratios , are independent of . We find

 T0=Cτ(coshρ)2/3,μ0=αCτ(coshρ)2/3,~μ0=~αCτ(coshρ)2/3, (12)

and therefore,

 ε0=f(α)C4τ4(coshρ)8/3,n0=αg(α)C3τ3cosh2ρ,~n0=~α~g(α)C3τ3cosh2ρ. (13)

The parameter is related to the particle multiplicity to be extracted from the experimental data. For a massless particle species (‘pion’), the relation is Gubser:2010ze (); Hatta:2014jva ()

 dNidY≈gi4C3π, (14)

where is the momentum rapidity and is the degeneracy factor.

### iii.2 Anisotropic ideal solution

We now perturb the solution anisotropically to introduce the dependence. In doing so, we shall focus on the early time regime (or , see (9)). As observed in Hatta:2014jva (), in this regime the perturbed solution is fully under analytical control including the viscous case to be discussed in the next section.

Following Gubser:2010ui (), we consider the following deformation of the isotropic solution

 ^ε0→^ε = ^ε0(1−ϵnAδ)4, ^n0→^n = ^n0(1−ϵnAδ′)3, (^uρ,^uΘ,^uϕ,^uζ)=(1,0,0,0) → (1,−ϵnνs^gΘΘ∂ΘA,−ϵnνs^gϕϕ∂ϕA,0), (15)

where

 A≡(2Lx⊥L2+x2⊥)ncosnϕ, (16)

is proportional to the spherical harmonics in the early time regime . Note that we preserve boost invariance in this paper, but the case was also considered in Gubser:2010ui (). is the eccentricity111In a conformal theory, the definition of eccentricity requires some care. We use Gubser:2010ui (); Hatta:2014jva () (17) . which we assume to be small and keep only linear terms in . and have to be determined by solving the hydrodynamic equations linearized around the isotropic solution. Plugging (15) into (3), we find the following equation for

 ∂ρδ′=νs3cosh2ρn(n+1). (18)

This turns out to be exactly the same as the equation satisfied by Gubser:2010ui (). Therefore, in the ideal case we have , which means that and are rescaled by the same factor , and . The ratios and are thus unchanged. At early times , the right hand side of (18) is negligible and we can set Hatta:2014jva ().

### iii.3 vn at finite μ

In order to compute flow harmonics , we use the Cooper-Frye formula Cooper:1974mv ()

 (2π)3dNidYpTdpTdϕp=gi∫Σ(−pμdσμ)(exp(u⋅p+kμiT)+δf)∝1+2vn(pT)cosnϕp, (19)

where we assumed the Boltzmann distribution and is the deviation from the equilibrium distribution. The use of the Boltzmann distribution (rather than the Fermi/Bose distributions) may be justified for the purpose of computing the integrated Hatta:2014jva (). generically represents a set of chemical potentials for net baryon number, isospin and strangeness. We assign for particles with positive/negative quantum numbers mentioned above, and for neutral particles with respect to the corresponding quantum number. In principle, since we are assuming conformal symmetry, the formula (19) should be used only for massless particles, or particles that can be approximately treated as massless (i.e., pions). However, for the sake of discussion in Section VI.3, we shall later introduce massive particles and compute their in the ‘probe approximation’, namely, by neglecting their backreaction to the flow velocity. Since we add in particles in the final state that do not exist in the fluid, the total energy is not conserved at freezeout. But the fraction of the change is exponentially suppressed by the mass and will be neglected.

The integral in (19) is taken over the hypersurface of constant energy density where the kinetic freezeout occurs. In the ideal case, constant means constant since is a constant. Let us write the condition of constant energy density as

 ε(τ,x⊥,ϕ)=T4f(α)≡C4B4(2L)4f(α)=εc. (20)

Typically, is of the order of the critical energy density of the QCD phase transition. We take in this paper. Following Hatta:2014jva (), we assume that the condition (20) is reached within the early time regime where we can use the approximate solution (15). The parameter in (20) is then related to the (position-dependent) freezeout time as

 τf(x⊥,ϕ)=(2L)5B3(L2+x2⊥)2(1−3ϵn(2Lx⊥L2+x2⊥)ncosnϕ). (21)

For consistency with our early freezeout scenario, we must have .

Under these assumptions, the integral (19) can be performed analytically and the integrated is obtained from the formula

 vn=∫dpTvn(pT)dNdYdpT∫dpTdNdYdpT. (22)

In the ideal case , does not depend on since the factor cancels in the ratio (22). The result is Hatta:2014jva ()222See (61) of Hatta:2014jva (). We have corrected a mistake by a factor of 2 in the overall normalization. The same comment applies to (48) below.

 vnϵn=964Γ(3n)Γ(4n)(128B3)nΓ2(n2)n2(3n+2)2(n−1)2(4n+1)∼B−3n∝(f3/4ε3/4cL3dNdY)n. (23)

This determines the dependence of . Quite generally, is an increasing function . On the other hand, is a decreasing function of . We shall see that, in heavy-ion collisions, the latter dependence is stronger, and as a result (23) is a decreasing function of , or equivalently, an increasing function of the collision energy . Incidentally, we note that the directed flow vanishes, consistently with our assumption of boost-invariance.

## Iv Viscous case

We now turn to the viscous case . Although the system is out of equilibrium, from the Landau matching condition we can define the local and using the same relations as in equilibrium

 ε=T4f(μT),n=μT2g(μT),~n=~μT2~g(μT), (24)

but now cannot be a constant.

### iv.1 Isotropic viscous solution

First consider the isotropic case . Although in (2) is not a constant anymore, it depends only on (see below). Then we still have so that

 n∝1τ3cosh2ρ, (25)

is the same as in the ideal case Gubser:2010ui (). However, the solution of the Navier-Stokes equation has an extra -dependence proportional to the shear viscosity . In the case of vanishing chemical potentials, this -dependence can be obtained exactly Gubser:2010ze ()

 εNS=1τ4fC4(coshρ)8/3[1+^η9f1/4Csinh3ρ 2F1(32,76,52;−sinh2ρ)]4, (26)

where is independent of .

However, at finite density, will depend on , and this makes it difficult to find an exact solution. Related to this, can now depend on both and , and this relation can be model-dependent. We can get around this problem by assuming that is small. Specifically, we rescale by the entropy density

 η=¯η(μT)s, (27)

as is often done in hydrodynamic simulations. We then regard as a small parameter () and keep only terms linear in . In this approximation, we may replace and in (27) by their equilibrium values at , namely, and

 s≈1T0(ε0+p0−μ0n0)=C3τ3cosh2ρ(43f(α)−α2g(α))≡C3τ3cosh2ρh(α). (28)

We then find the solution valid to

 εNS = T4f(μT)≈f(α)C4τ4(coshρ)8/3[1+4h(α)¯η(α)9f(α)Csinh3ρ 2F1(32,76,52;−sinh2ρ)] (29) ≈f(α)C4τ4(coshρ)8/3[1−2h(α)¯η(α)f(α)C(e−ρ2)2/3], nNS = μT2g(μT)=αg(α)C3τ3cosh2ρ, (30)

where in the second line of (29) we focus on the early-time regime where is negative and large.333The viscous Gubser solution is known to become unphysical (the temperature becomes negative) as Gubser:2010ze (). Physically, this corresponds to very early times and/or very large values of . Our results are not sensitive to these regions. We can simply choose the initial time of the evolution to be small, but not too small. Besides, all the -integrals to be performed below are fully convergent at .

Using (29) and (30), we can eliminate

 (μTg(μT))4/3f(μT)=(αg(α))4/3f(α)(1−2h(α)¯η(α)f(α)C(e−ρ2)2/3). (31)

Writing

 μT=α+δα(ρ), (32)

we find the deviation from constancy due to the viscosity

 δα(ρ)≈2h¯ηCf(e−ρ2)2/343α+4g′3g−f′f=γh¯ηCf(L2+x2⊥2Lτ)2/3, (33)

where

 γ(α)≡243α+4g′3g−f′f. (34)

Note that as . At the freezeout time , we have the relation

 μT=α+δα|freezeout≈α+γK(L2+x2⊥)2(2L)4. (35)

Finally, we can solve for and using (32). The result is

 T=Cτ(coshρ)2/3(1−δα3(1α+g′g)), μ=αCτ(coshρ)2/3(1+δα3(2α−g′g)). (36)

### iv.2 Anisotropic viscous solution

We now perturb the solution as in (15). First consider the current in (2). now depends not only on , but also on and . However, the dependence is of order . (See (43) below. Remember that for the ideal solution is constant even in the anisotropic case.) Therefore, if we neglect terms of order , we can approximate . Then (18) is still valid and we get

 n=μT2g(μT)≈αg(α)C3τ3cosh2ρ(1−ϵnA)3. (37)

As for the energy density, we find

 ε=T4f(μT)≈f(α)C4τ4(coshρ)8/3(1−2h(α)¯η(α)f(α)C(e−ρ2)2/3)(1−ϵnAδ)4, (38)

where Hatta:2014jva ()

 δ≈1+h(α)¯η(α)2f(α)C(e−ρ2)2/3. (39)

From the constant energy condition

 ε=C4B4(2L)4f(α)=εc, (40)

we can determine the freezeout surface in the viscous case Hatta:2014jva ()

 τf(x⊥,ϕ)=(2L)5B3(L2+x2⊥)2(1−3K(L2+x2⊥)22(2L)4−3ϵn(2Lx⊥L2+x2⊥)ncosnϕ), (41)

where the ‘Knudsen number’ is proportional to the shear viscosity

 K=h(α)¯η(α)B2f(α)C. (42)

(32) and (36) are also modified as where

 δα′=δα(1+ϵnA)=δα(1+ϵn(2Lx⊥L2+x2⊥)ncosnϕ), (43)

and

 T = μ = αCτ(coshρ)2/3[1−ϵnA+δα3(2α−g′g)]. (44)

### iv.3 vn at finite μ and η

The computation of is more complicated than the case. This is because does not mean , and therefore one cannot treat in the Boltzmann factor (19) as a constant when integrating over the hypersurface of constant energy. In order to cope with this, we write (IV.2) as

 T=Tc−T0f′4fδα, (45)

where

 Tc≡CB2L=T0⎛⎝1−ϵnA−h¯η2fC(L2+x2⊥2Lτ)2/3⎞⎠, (46)

is constant by virtue of (40). We then expand the Boltzmann factor as444In this subsection we set . The case with will be treated in the next subsection.

 exp(u⋅p+kμT) ≈ ekαeu⋅p/Tcexp[δα(u⋅pT2cT0f′4f+k(1+ϵnA))] (47) ≈ ekαeu⋅p/Tc{1+δα(u⋅pTcf′4f+k)(1+ϵnA)},

where we approximated in the term.

The first term in (47), proportional to unity, gives the same result as in Hatta:2014jva ()555For simplicity, here we ignore the contribution from the nonequilibrium part in (19). This has been computed in Hatta:2014jva () for a particular choice of . However, its -dependence is strongly affected by the choice of which is not unique Dusling:2009df (). Moreover, even the overall sign of this contribution is sensitive to the -cutoff.

 vnϵn=964Γ(3n)Γ(4n)(128B3)nΓ2(n2){n2(3n+2)2(n−1)2(4n+1)−3n3(n−1)K16(3n−1)(3n2+3n+2)}. (48)

Note that for (see, however, Staig:2011wj ()). The second term in (47) leads to a new order contribution to . To compute it, we borrow some results from Hatta:2014jva (). First, the perturbed flow velocity on the freezeout surface has the following components in the coordinates (7)

 u⊥=u0⊥+δu⊥ϵncosnϕ,uϕ=δuϕϵnsinnϕ, (49)

where

 u⊥0=2x⊥(2L)5B3(L2+x2⊥)3, δu⊥=3(2L)5B3(L2+x2⊥)4(2Lx⊥L2+x2⊥)n−1L(n(L2−x2⊥)−4x2⊥), δuϕ=−3n2(2L)5B3(L2+x2⊥)2(2Lx⊥L2+x2⊥)n. (50)

(The viscosity can be neglected here.) The exponential factor in the Boltzmann distribution reads

 p⋅uTc=1Tc[−mTcosh(ζ−Y)+pTu⊥cos(ϕ−ϕp)−pTuϕx⊥sin(ϕ−ϕp)]≡U+ϵnδU, (51)

where is the transverse mass. The volume element of the constant energy hypersurface is

 −pμdσμ=x⊥τf(mTcosh(ζ−Y)−pTcos(ϕ−ϕp)∂τf∂x⊥+pTx⊥sin(ϕ−ϕp)∂τf∂ϕ)dζdx⊥dϕ,

where is given by (41) with the viscous term set to zero. Finally, we need the more precise version of (35)

 δα≈γK(L2+x2⊥)2(2L)4(1+2ϵnA). (53)

Armed with these formulas, let us decompose the contribution from the second term in (47) as

 (2π)3dNdYpTdpTdϕp ∼ ekα∫Σ(−pμdσμ)eU+ϵnδUδα((U+ϵnδU)f′4f+k)(1+ϵnA) (54) ≡ (δJ1+δJ2+δJ3)ϵncosnϕp,

corresponding to the three terms in (LABEL:sur). Consider first. To we have to evaluate

 δJ1∼ϵn2LγKekαB3mT∫dx⊥x⊥∫dζdϕcosh(ζ−Y)eUδU(f′4fU+f′4f+k). (55)

This can be efficiently evaluated using the trick introduced in Hatta:2014jva () (see Eq. (73) there). The -integral gives Bessel functions where

 z≡pTu⊥0Tc=2x⊥pT(2L)5TcB3(L2+x2⊥)3. (56)

This can be expanded as anticipating that the subsequent -integral is dominated by the region . We thus find

 δJ1 ≈ 2LγKekαB3mT∫∞0dx⊥x⊥4πzn2n(n−1)!1u⊥0(δu⊥−δuϕx⊥) (57) ×[f′4f(nK1(mT/Tc)−mT2T(K0(mT/Tc)+K2(mT/Tc)))+kK1(mT/Tc)] = 2LγKekαB39πL2mT(64pTTcB3)nn(n−1)Γ(3n)(3n−1)Γ(4n) ×[f′4f(nK1(mT/Tc)−mT2T(K0(mT/Tc)+K2(mT/Tc)))+kK1(mT/Tc)].

The correction to can be calculated from the formula (cf. (22))

 δv1n≡∫∞0dpTpTδJ1∫∞0dpTpTJ0ϵn2, (58)

where is the azimuthally symmetric part (cf. Eq. (45) of Hatta:2014jva ())

 J0=4πmTK1(mT/Tc)16L3B3ekα. (59)

In the massless case , the integral can be done exactly and we find

 δv1nϵn=9γK128(k−3f′4f)(128B3)nn2(n−1)Γ(3n)(3n−1)Γ(4n)Γ(n2+2)Γ(n2), (60)

and from (23),

 δv1nvidealn=γK4(k−3f′4f)n(n+2)(4n+1)(3n−1)(3n+2)2. (61)

It is important to emphasize that (61) is induced by the combined effect of the chemical potential and the viscosity. It vanishes when or because (cf. (34)). Compared with (48) which schematically reads , we notice that (61) is not enhanced by a factor of , hence subleading at large . However, it is the leading contribution to the difference in between particles () and antiparticles (). If is the baryon chemical potential, the protons have larger than the antiprotons. We shall study this effect in detail later.

In fact, for protons the approximation is not valid. Instead, we now assume and reevaluate . Note that when , is parametrically larger than , so it is enough to consider only .

When , the Bessel function is independent of the order

 Ki(mT/Tc)≈√π