Lattice-based equation of state at finite baryon number, electric charge and strangeness chemical potentials

We construct an equation of state for Quantum Chromodynamics (QCD) at finite temperature and chemical potentials for baryon number $B$, electric charge $Q$ and strangeness $S$. We use the Taylor expansion method, up to the fourth power for the chemical potentials. This requires the knowledge of all diagonal and non-diagonal $BQS$ correlators up to fourth order: these results recently became available from lattice QCD simulations, albeit only at a finite lattice spacing $N_t=12$. We smoothly merge these results to the Hadron Resonance Gas (HRG) model, to be able to reach temperatures as low as 30 MeV; in the high temperature regime, we impose a smooth approach to the Stefan-Boltzmann limit. We provide a parameterization for each one of these $BQS$ correlators as functions of the temperature. We then calculate pressure, energy density, entropy density, baryonic, strangeness, electric charge densities and compare the two cases of strangeness neutrality and $\mu_S=\mu_Q=0$. We also calculate the isentropic trajectories and compare them in the two cases. Our equation of state can be readily used as an input of hydrodynamical simulations of matter created at the Relativistic Heavy Ion Collider (RHIC).


INTRODUCTION
Relativistic heavy ion collisions have successfully recreated the Quark Gluon Plasma (QGP) in the laboratory at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the Large Hadron Collider (LHC) at CERN. At low baryon densities, the transition from the hadron gas phase where quarks and gluons are confined within hadrons into a deconfined state where quark and gluons are the main degrees of freedom is a smooth cross-over [1][2][3]. At larger baryon densities, the phase transition is expected to become stronger, eventually turning into first-order. If this is the case, there has to be a critical point on the QCD phase diagram [4][5][6][7][8]. The search for the QCD critical point is the focus of the second Beam Energy Scan (BES II) at RHIC, running in 2019 and 2020.
The Quark Gluon Plasma acts as a nearly perfect fluid and as such can be well-described by event-by-event relativistic viscous hydrodynamical models. The hydrodynamical description of the fireball has proved to be very successful in describing the experimental data [9][10][11][12][13][14][15][16][17][18][19]. In order to close the hydrodynamical equations, an Equation of State (EoS) is required, which is based on first principle Lattice QCD calculations. Recently, a Bayesian analysis [20] has provided an important validation of the lattice QCD equation of state. This framework, based on a comparison of data from the LHC to theoretical models, has applied state-of-the-art statistical techniques to the combined analysis of a large number of observables while varying the model parameters. The posterior distribution over possible equations of states turned out to be consistent with results from lattice QCD simulations. Additionally, the correct description of the QCD equation of state is needed because differences in the equation can affect the extraction of transport coefficients [17]. Thus, a lattice-based QCD equation of state is a fundamental ingredient in the description of the state of matter created in a heavy-ion collision. The precise lattice QCD results for several thermodynamic quantities can thus be used in support of the heavy ion experimental program [21].
The EoS of QCD at zero baryonic density is known with high precision from first principles since a few years [22][23][24]. The calculation of the equation of state at finite chemical potential is hindered by the sign problem. Nevertheless, the thermodynamic quantities can be expanded as a Taylor series in powers of µ B /T , for which the coefficients χ n can be simulated on the lattice at µ B = 0. From these Taylor coefficients a variety of Lattice QCD based equations of state have been reconstructed [25][26][27] and later used within relativistic hydrodynamics [25,[28][29][30][31].
However, baryon number is not the only conserved charge in a heavy ion collision: strangeness and electric charge are also relevant quantum numbers. In fact, many questions remain regarding a possible separate freeze-out temperature for strange hadrons [32][33][34][35] and separations of electric charge due to a possible chiral magnetic effect [36], so many interesting questions need to be answered, that go beyond just baryon charge conservation. At the LHC, where the baryonic chemical potential µ B is basically vanishing, the chemical potentials for strangeness µ S and electric charge µ Q are also zero. At RHIC how-ever, as the baryonic density increases, the other two chemical potentials have finite values as well. Until now, the equation of state of QCD has only been extrapolated to finite µ B , either by keeping µ S = µ Q = 0, or along a specific trajectory in the four-dimensional parameter space, namely imposing that the strangeness density n S = 0 and that the electric charge density n Q = 0.4 n B to match the experimental situation.
After the early results for χ 2 , χ 4 and χ 6 [37], a continuum extrapolation for χ 2 was published in Ref. [38]; in Ref. [39] χ 4 was shown, but only at finite lattice spacing. The continuum limit for χ 6 was published for the first time in [40] in the case of strangeness neutrality, and later in [41] for both cases. In [42], a first determination of χ 8 at two values of the temperature and N t = 8 was presented. Finally, in Ref. [43] a determination of χ 8 was presented for the first time as a function of the temperature, at N t = 12, keeping µ S = µ Q = 0. Recently, the effect of introducing a critical point in the equation of state of QCD has also been tested [26].
However, a Taylor expansion of the equation of state, along a direction which satisfies the strangenessneutrality condition is not enough for the hydrodynamics approach, since the fluid cells have local fluctuations in strangeness density. Additionally, there is a complicated interplay between transport coefficients when B, Q, S are considered [44] that cannot be neglected at large baryon densities. For these reasons, an EoS fully expanded as a Taylor series in powers of µ B /T, µ S /T, µ Q /T is needed as an input of hydrodynamic simulations of the matter created at RHIC. In order to perform such an expansion, all of the diagonal and non-diagonal susceptibilities of these three conserved charges are needed from lattice QCD up to the chosen power. In this work, we perform the Taylor expansion of to total power four in the chemical potentials. These results recently became available [43], on N t = 12 lattices.
Alternative approaches to the Taylor series expansion have been suggested in [45,46] and [47,48], which have been shown to match well to lattice QCD data for the Fourier harmonics [49] at imaginary chemical potential. These Fourier harmonics appear to be important to distinguish baryon interactions within a hadron resonance gas (see also [50]), specifically for the thermodynamic regime above T > 150 MeV. We note that here we use lattice QCD data entirely in this regime (our hadron resonance gas model is only to constrain low temperatures below T 135 MeV where no lattice QCD results are available). However, due to the Taylor expansion our ap-proach is limited to chemical potentials µ B (2 − 2.5) T . To fully reproduce the Fourier harmonics we would need to reach µ B πT , for which higher order coefficients in the Taylor series would need to be included.
In this manuscript, we construct an equation of state for QCD at finite T , µ B , µ S , µ Q . We build the pressure as a Taylor series of the three chemical potentials, with coefficients taken from lattice simulations [43]. At low temperatures, we perform a smooth merging between the lattice and the Hadron Resonance Gas model results [51] and ensure continuity of higher order derivatives. At high temperatures, we impose a smooth approach to the Stefan-Boltzmann limit. We parameterize each one of these coefficients as a ratio of polynomials. From this we obtain the pressure and can then calculate all other quantities from thermodynamic relationships 1 .

METHODOLOGY AND RESULTS
The Taylor series of the pressure in terms of the three conserved charge chemical potentials can be written as We limit our calculation to i + j + k ≤ 4. The coefficients have recently been published from lattice QCD simulations on 48 3 × 12 lattices [43] in the temperature range (135 MeV)< T < (220 MeV). Since this is not enough to cover the hydrodynamical evolution of the system, we smoothly merge each coefficient at low temperature with the Hadron Resonance Gas model result, while at high temperature we calculate the Stefan-Boltzman limit for each one of them and assume that their value at T = 800 MeV is ∼ 10% away from the respective Stefan-Boltzmann limit. To simplify the notation, whenever i, j, k are zero, we only write the non-zero indices and only the corresponding conserved charges: for example, χ BQS 200 becomes χ B 2 , χ BQS 301 becomes χ BS 31 and so on. In order to provide a smooth pressure which can be easily derived to obtain the other thermodynamic quantities, we parameterize each coefficient by means of a ratio of up-to-ninth order polynomials in the inverse temperature: Only χ B 2 requires a different parameterization: In both equations above, t = T /154 MeV, t = T /200 MeV [53]. The values of the parameters for each coefficient are given in the appendix, together with the respective Stefan-Boltzmann limits. Figures 1 and 2 show all of the Taylor expansion coefficients as functions of the temperature. The black dots are the HRG model results, the red triangles correspond to the lattice QCD results and the thick blue line indicates the Stefan-Boltzmann limit.
Making use of this parameterization, we construct the pressure from Eq. (1). The other thermodynamic quantities are then derived from the pressure as follows: Everywhere in the above equation, i = j is intended.
In Fig. 3 we show the dependence of the normalized pressure, entropy density, energy density, baryonic, strangeness and electric charge densities on the temperature, along lines of constant µ B /T = 0.5, 1, 2, both with n S = 0, n Q = 0.4 n B (solid black lines), and in the case µ S = µ Q = 0 (dashed red lines). We find that the thermodynamic quantities that are less sensitive to the chemical composition of the system do not show large discrepancies between the two scenarios, for all three values of µ B /T . On the other hand, when realistic conditions on the global chemical composition of the system are imposed, the baryon density is largely affected, and substantially decreased; the opposite effect is visible for the electric charge density, which is heavily enhanced.
Finally, we compare i) the isentropic trajectories, ii) the temperature dependence of the speed of sound along lines of constant µ B /T and iii) the behavior of the speed of sound along parametrized chemical freeze-out lines between these two cases. The isentropic trajectories are shown in Fig. 4 for selected values of s/n B , which correspond to the indicated collision energies [40]. In the upper panel of Fig. 5 we show the speed of sound as a function of the temperature along lines with µ B /T = 0.5, 1, 2; the different colors correspond to different values of µ B /T . In the lower panel of Fig. 5 we show the behavior of the speed of sound along two parametrized chemical freeze-out lines. These two freeze-out lines are shifted from the one presented in [54], and have the form: Fig. 4 and in Fig. 5, the solid lines correspond to n S = 0, n Q = 0.4 n B while the dashed lines to µ S = µ Q = 0.
Since the EoS constructed in this work is a Taylor expansion carried out from lattice-QCD-calculated expansion coefficients, it is important to have an idea of the range of the validity of such expansion. It has been shown from lattice QCD simulations that the Taylor expansion of the Equation of State up to O(µ 4 B ) converges for µ B /T 2 − 2.5 [41], and the same can be said for our EoS. This roughly corresponds to a collision energy of √ s 10 GeV [54]. In order to have a better idea of where a possible breakdown of its validity occurs, we show in Fig. 6 the behavior of the electric chemical potential in the case with strangeness neutrality, along lines of constant µ B /T = 0.5 − 3. We see that a non-monotonic behavior appears around and above µ B /T ∼ 2.5. This is in line with the expectation that the convergence of the Taylor series is guaranteed in the regime µ B /T 2.5. We note again that with the Taylor expansion approach used here, we do not expect to fully incorporate the constraints from imaginary µ B -and thus reproduce the Fourier harmonics from [49] -since for them the coverage of the region µ B /T ≤ π would be required. Applying the constraints from imaginary µ B can be done in the near future to further improve our modeling of the QCD EoS, possibly concurrently with the inclusion of new continuum extrapolated lattice results.

CONCLUSIONS
In this manuscript, we constructed an equation of state for QCD at finite temperature and B, Q, S chemical potentials, based on a Taylor series up to fourth power in the chemical potentials. Our methodology is based on a smooth merging between the HRG model and lattice QCD results for each one of the Taylor expansion coefficients; for all coefficients except χ B 2 , the parameterization function is a ratio of up-to-ninth order polynomials. We provide all parameters in Tables I,II,III, so that our EoS can be readily used in the community. Furthermore, the code to generate the EoS and the tables for the thermodynamic quantities as functions of T, µ B , µ S , µ Q is available at the link mentioned in Ref. [55]. The Equation of State presented in this manuscript is important for the hydrodynamic description of the system created in heavy ion collisions at RHIC. There are numerous outstanding questions that remain to be understood at finite baryon densities that are influenced both by electric charge and strangeness. One recent surprise that arose from the first Beam Energy Scan was Λ polarization, that indicates that the Quark Gluon Plasma may be the most vortical fluid known to humanity [56]. However, considering that Λ's are simultaneously both strange particles and baryons, polarization studies should be done in hydrodynamic simulations that also consider all three conserved charges because of this interplay between strangeness and baryon number. As previously mentioned, this BQS equation of state can help shed light on the possible flavor hierarchy of freeze-out temperatures as well as the chiral magnetic effect. A variety of dynamical observables of conserved charges (e.g. kaon flow harmonics) have already been measured at the Beam Energy Scan I and many others are planned for the Beam Energy Scan II, which may help to further constrain the location of a possible critical point.
Finally, we point out that strange hadrons make up roughly 10% of all measured hadrons (assuming the kaon to pion ratio is a reasonable estimate for the ratio of all final state hadrons) and we can primarily only measure charged particles 2 . Thus, a BQS equation of state is required for a fully consistent description of the Quark Gluon Plasma at finite densities. Relativistic hydrodynamics in the presence of multiple conserved charges obtains cross terms that affect the transport coefficients [44,57,58]. Thus, it is misleading to extract transport   coefficients at finite baryon densities only considering finite baryon number and not also finite strangeness and electric charge. Furthermore, transport coefficients of different conserved charges have different characteristic temperatures, which further complicates the picture at large densities [59]. The consequences are still under development, but it is certain that a BQS equation of state is a vital first step to take into account any of these ef-fects.
At this point, our reconstructed BQS equation of state only consists of a cross-over transition. Unlike a previous work where an equation of state at finite µ B was coupled to the 3D Ising model in order to study criticality [26], such an endeavor with three conserved charges would be significantly more complicated. While the term "critical point" is used, there might actually be a critical line or even critical plane once one considers the full three dimensional space of µ B , µ S , and µ Q . Since there are large fluctuations in T , µ B , µ S , and µ Q throughout the evolution of a single event [60][61][62], certain elements of the fluid might pass through a critical region at an entirely different combination of T , µ B , µ S , and µ Q .

APPENDIX
We list the values of the parameters in Eq. (3) for each Taylor expansion coefficient in Table I. The Stefan-Boltzmann limit for the coefficients have the following values:  The " − " symbols in the table indicate that, for most of the coefficients, it is enough to consider a ratio of polynomials of order lower than nine.