Equation of state at finite densities for QCD matter in nuclear collisions

TL;DR

This paper develops a QCD equation of state at finite densities incorporating conserved charges, crucial for modeling heavy ion collisions, by combining lattice QCD data with the hadron resonance gas model.

Contribution

It introduces a comprehensive QCD equation of state at finite chemical potentials considering multiple conserved charges, tailored for relativistic heavy ion collision applications.

Findings

Finite electric charge and strangeness chemical potentials significantly affect particle production.

Hydrodynamic simulations demonstrate the impact of these chemical potentials on collision outcomes.

The model aligns with conditions of strangeness neutrality and fixed charge-to-baryon ratios in nuclei.

Abstract

We construct the QCD equation of state at finite chemical potentials including net baryon, electric charge, and strangeness, based on the conserved charge susceptibilities determined from lattice QCD simulations and the equation of state of the hadron resonance gas model. For the application to relativistic heavy ion collisions we consider the situation of strangeness neutrality and matter with a fixed electric charge-to-baryon ratio, resembling that of heavy nuclei. The importance of finite electric charge and strangeness chemical potentials for particle production in heavy ion collisions is demonstrated using hydrodynamic simulations.

Tables

Table 1. Table 1: Ratios of protons to nucleons Z / A 𝑍 𝐴 Z/A for the nuclei used in the collider experiments at RHIC and LHC.

Nucleus	$Z/A$
${}_1{}^{}{}_{}{}^1$H	1.000
${}_1{}^{}{}_{}{}^2$H	0.500
${}_2{}^{}{}_{}{}^3$He	0.667
${}_{13}{}^{}{}_{}{}^{27}$Al	0.481
${}_{29}{}^{}{}_{}{}^{63}$Cu	0.460
${}_{40}{}^{}{}_{}{}^{96}$Zr	0.417
${}_{44}{}^{}{}_{}{}^{96}$Ru	0.458
${}_{54}{}^{}{}_{}{}^{127}$Xe	0.425
${}_{79}{}^{}{}_{}{}^{197}$Au	0.401
${}_{82}{}^{}{}_{}{}^{208}$Pb	0.394
${}_{92}{}^{}{}_{}{}^{238}$U	0.387

Table 2. Table 2: The list of parameters used for the parametrization of the second-order susceptibilities.

Susceptibility	$h_0$	$h_1$ (GeV2)	$T_1$ (GeV)	$\Delta T_1$ (GeV)	$T_2$ (GeV)	$\Delta T_2$ (GeV)	$h_2$ (GeV-n)	$n$
${\chi }_2^B$	$3.37\times {10}^{-1}$	$9.65\times {10}^{-3}$	$1.73\times {10}^{-1}$	$2.13\times {10}^{-2}$	$1.69\times {10}^{-1}$	$1.57\times {10}^{-2}$	$3.42\times {10}^5$	$8$
${\chi }_2^Q$	$6.71\times {10}^{-1}$	$6.45\times {10}^{-3}$	$1.40\times {10}^{-1}$	$2.87\times {10}^{-2}$	$1.00\times {10}^{-1}$	$1.00\times {10}^{-3}$	$1.00\times {10}^4$	$6$
${\chi }_2^S$	$1.02\times {10}^0$	$1.55\times {10}^{-2}$	$1.69\times {10}^{-1}$	$3.35\times {10}^{-2}$	$1.52\times {10}^{-1}$	$3.17\times {10}^{-2}$	$1.24\times {10}^6$	$8$
${\chi }_{1,1}^{B,Q}$	$1.97\times {10}^{-4}$	$-2.09\times {10}^0$	$1.34\times {10}^{-1}$	$7.28\times {10}^{-2}$	$1.59\times {10}^{-1}$	$2.80\times {10}^{-2}$	$9.00\times {10}^4$	$8$
${\chi }_{1,1}^{B,S}$	$-3.38\times {10}^{-1}$	$1.28\times {10}^{-2}$	$1.64\times {10}^{-1}$	$3.16\times {10}^{-2}$	$9.99\times {10}^{-2}$	$2.74\times {10}^{-2}$	$-3.42\times {10}^5$	$8$
${\chi }_{1,1}^{Q,S}$	$3.39\times {10}^{-1}$	$1.46\times {10}^{-2}$	$1.79\times {10}^{-1}$	$3.33\times {10}^{-2}$	$1.60\times {10}^{-1}$	$2.90\times {10}^{-2}$	$4.34\times {10}^5$	$8$

Table 3. Table 3: The list of parameters used for the parametrization of the fourth-order susceptibilities.

Susceptibility	$h_3$	$h_4$ (GeV)	$h_5$ (GeV2)	$T_3$ (GeV)	$\Delta T_3$ (GeV)	$h_6$	${\chi }_2$
${\chi }_4^B$	$1.45\times {10}^{-2}$	$2.49\times {10}^{-3}$	$0$	$1.62\times {10}^{-1}$	$2.27\times {10}^{-2}$	$1.00\times {10}^0$	${\chi }_2^B$
${\chi }_4^Q$	$1.35\times {10}^{-1}$	$0$	$0$	$1.61\times {10}^{-1}$	$1.74\times {10}^{-2}$	$1.25\times {10}^0$	${\chi }_2^Q$
${\chi }_4^S$	$6.36\times {10}^{-1}$	$-1.12\times {10}^{-1}$	$2.09\times {10}^{-2}$	$1.65\times {10}^{-1}$	$1.93\times {10}^{-2}$	$8.85\times {10}^{-1}$	${\chi }_2^S$
${\chi }_{3,1}^{B,Q}$	$0$	$0$	$0$	$1.63\times {10}^{-1}$	$1.16\times {10}^{-2}$	$9.96\times {10}^{-1}$	${\chi }_{1,1}^{B,Q}$
${\chi }_{2,2}^{B,Q}$	$4.42\times {10}^{-2}$	$1.31\times {10}^{-3}$	$-4.79\times {10}^{-4}$	$1.59\times {10}^{-1}$	$1.42\times {10}^{-2}$	$7.95\times {10}^{-1}$	${\chi }_2^B$
${\chi }_{1,3}^{B,Q}$	$4.25\times {10}^{-2}$	$4.54\times {10}^{-3}$	$-1.91\times {10}^{-3}$	$1.58\times {10}^{-1}$	$1.70\times {10}^{-2}$	$8.79\times {10}^{-1}$	${\chi }_2^B$
${\chi }_{3,1}^{B,S}$	$-2.87\times {10}^{-2}$	$7.93\times {10}^{-3}$	$-1.90\times {10}^{-3}$	$1.62\times {10}^{-1}$	$2.18\times {10}^{-2}$	$6.60\times {10}^{-1}$	${\chi }_{1,1}^{B,S}$
${\chi }_{2,2}^{B,S}$	$7.87\times {10}^{-2}$	$-1.35\times {10}^{-2}$	$2.60\times {10}^{-3}$	$1.68\times {10}^{-1}$	$2.46\times {10}^{-2}$	$-8.80\times {10}^{-1}$	${\chi }_{1,1}^{B,S}$
${\chi }_{1,3}^{B,S}$	$-2.04\times {10}^{-1}$	$1.85\times {10}^{-3}$	$-7.88\times {10}^{-4}$	$1.62\times {10}^{-1}$	$1.97\times {10}^{-2}$	$9.85\times {10}^{-1}$	${\chi }_{1,1}^{B,S}$
${\chi }_{3,1}^{Q,S}$	$2.31\times {10}^{-2}$	$-9.73\times {10}^{-4}$	$3.42\times {10}^{-4}$	$1.60\times {10}^{-1}$	$3.06\times {10}^{-2}$	$1.08\times {10}^0$	${\chi }_2^B$
${\chi }_{2,2}^{Q,S}$	$6.88\times {10}^{-2}$	$-2.24\times {10}^{-3}$	$9.64\times {10}^{-4}$	$1.63\times {10}^{-1}$	$2.60\times {10}^{-2}$	$1.12\times {10}^0$	${\chi }_2^B$
${\chi }_{1,3}^{Q,S}$	$2.02\times {10}^{-1}$	$1.04\times {10}^{-3}$	$-6.41\times {10}^{-4}$	$1.80\times {10}^{-1}$	$3.08\times {10}^{-2}$	$1.16\times {10}^0$	${\chi }_2^B$
${\chi }_{2,1,1}^{B,Q,S}$	$2.24\times {10}^{-2}$	$9.45\times {10}^{-5}$	$-2.33\times {10}^{-5}$	$1.62\times {10}^{-1}$	$1.30\times {10}^{-2}$	$5.81\times {10}^{-2}$	${\chi }_2^B$
${\chi }_{1,2,1}^{B,Q,S}$	$-2.30\times {10}^{-2}$	$1.00\times {10}^{-3}$	$-4.84\times {10}^{-4}$	$1.54\times {10}^{-1}$	$1.51\times {10}^{-2}$	$-1.39\times {10}^{-1}$	${\chi }_2^B$
${\chi }_{1,1,2}^{B,Q,S}$	$-6.72\times {10}^{-2}$	$-6.89\times {10}^{-4}$	$3.00\times {10}^{-4}$	$1.63\times {10}^{-1}$	$1.66\times {10}^{-2}$	$-1.07\times {10}^{-1}$	${\chi }_2^B$

Table 4. Table 4: The list of parameters used for the parametrization of the sixth-order susceptibilities.

Susceptibility	$h_7$	$h_8$	$h_9$	$h_{10}$	$T_4$ (GeV)	$\Delta T_4$ (GeV)	$T_5$ (GeV)
${\chi }_6^B$	$7.54\times {10}^{-2}$	$2.70\times {10}^{-2}$	$-1.64\times {10}^{-2}$	-	$1.27\times {10}^{-1}$	$1.73\times {10}^{-2}$	$1.57\times {10}^{-1}$
${\chi }_{5,1}^{B,Q}$	$2.59\times {10}^{-2}$	$1.39\times {10}^{-2}$	$1.81\times {10}^{-2}$	$8.73\times {10}^{-4}$	$1.21\times {10}^{-1}$	$1.12\times {10}^{-2}$	$1.52\times {10}^{-1}$
${\chi }_{5,1}^{B,S}$	$-5.52\times {10}^{-2}$	$4.38\times {10}^{-3}$	$-6.94\times {10}^{-3}$	-	$1.00\times {10}^{-1}$	$7.50\times {10}^{-3}$	$1.52\times {10}^{-1}$
	$\Delta T_5$ (GeV)	$T_6$ (GeV)	$\Delta T_6$ (GeV)	$T_7$ (GeV)	$\Delta T_7$ (GeV)	$T_8$ (GeV)	$\Delta T_8$ (GeV)
${\chi }_6^B$	$1.09\times {10}^{-2}$	$2.17\times {10}^{-1}$	$5.12\times {10}^{-2}$	$2.63\times {10}^{-1}$	$1.43\times {10}^{-2}$	-	-
${\chi }_{5,1}^{B,Q}$	$1.11\times {10}^{-2}$	$1.64\times {10}^{-1}$	$7.24\times {10}^{-3}$	$1.96\times {10}^{-1}$	$2.58\times {10}^{-2}$	$2.49\times {10}^{-1}$	$1.55\times {10}^{-2}$
${\chi }_{5,1}^{B,S}$	$1.20\times {10}^{-2}$	$1.34\times {10}^{-1}$	$1.07\times {10}^{-2}$	$1.72\times {10}^{-1}$	$1.13\times {10}^{-2}$	$2.02\times {10}^{-1}$	$1.81\times {10}^{-2}$

Equations

$$\frac{P}{T^4}$$

$${\chi }_{l,m,n}^{B,Q,S}={\frac{{\partial }^l{\partial }^m{\partial }^{n}P(T,{\mu }_B,{\mu }_Q,{\mu }_S)/T^4}{\partial ({\mu }_B/T{)}^l\partial ({\mu }_Q/T{)}^m\partial ({\mu }_S/T{)}^n}|}_{{\mu }_{B,Q,S}=0}.$$

$$P$$

$$\frac{P}{T^4}$$

$${\tilde{P}}_{\mathrm{lat}}(T,{\mu }_J)=P_{\mathrm{lat}}(T_s,{\mu }_J)\times \frac{T^4}{T_s^4},$$

$$\frac{{\partial }^{2}P}{\partial T^2}$$

$$\frac{{\partial }^{2}P}{\partial {\mu }_J^2}$$

$$\frac{{\partial }^{2}P}{\partial T^2}$$

$$\frac{{\partial }^{2}P}{\partial {\mu }_J^2}$$

$$s$$

$$n_J$$

$$e$$

$$c_s^2$$

$$f(T,{\mu }_B)=\tanh [(T-T_c({\mu }_B))/\Delta T_c],$$

$$c_s^2$$

$$\frac{P}{T^4}$$

$${\mu }_u$$

$${\mu }_d$$

$${\mu }_s$$

$${\chi }_2^B$$

$${\chi }_{1,1}^{B,Q}$$

$${\chi }_4^B$$

$${\chi }_{3,1}^{B,S}$$

$${\chi }_{3,1}^{B,Q}$$

$${\chi }_{3,1}^{Q,S}$$

$${\chi }_{2,1,1}^{B,Q,S}$$

$$\left(\begin{array}{c}{n}_B\\ n_Q\\ n_S\end{array}\right)=T^2\left(\begin{array}{ccc}{\chi }_2^B& {\chi }_{1,1}^{B,Q}& {\chi }_{1,1}^{B,S}\\ {\chi }_{1,1}^{B,Q}& {\chi }_2^Q& {\chi }_{1,1}^{Q,S}\\ {\chi }_{1,1}^{B,S}& {\chi }_{1,1}^{Q,S}& {\chi }_2^S\end{array}\right)\left(\begin{array}{c}{\mu }_B\\ {\mu }_Q\\ {\mu }_S\end{array}\right).$$

$${\mu }_B$$

$${\mu }_Q$$

$${\mu }_S$$

$${\chi }_2$$

$$\displaystyle g_{i}^{\pm}(T_{i},\Delta T_{i})=\frac{1}{2}\bigg{[}1\pm\tanh\bigg{(}\frac{T-T_{i}}{\Delta T_{i}}\bigg{)}\bigg{]}.$$

$${\chi }_4$$

$${\chi }_6^B$$

$${\chi }_{5,1}^{B,Q}$$

$${\chi }_{5,1}^{B,S}$$

Full text

I Introduction

The nearly perfect fluidity of the quantum chromo dynamic (QCD) matter discovered in heavy ion collisions at the BNL Relativistic Heavy Ion Collider (RHIC) and CERN Large Hadron Collider (LHC) has provided us with rare opportunities to experimentally explore the nuclear equation of state, which encodes the fundamental thermodynamic properties of nuclear matter. The system created in these collisions is expected to be deconfined from hadronic matter to the quark-gluon plasma (QGP) above approximately two trillion degrees Kelvin.

The theoretical study of the non-perturbative QCD equation of state dates back to the MIT bag model Chodos et al. (1974a, b) where hadrons are assumed to be located in a bag embedded in the QCD vacuum characterized by the bag constant. This prescription introduces confinement phenomenologically. Also, several model approaches, such as the potential model De Rujula et al. (1975) and the Nambu-Jona-Lasinio model Nambu and Jona-Lasinio (1961a, b), have been proposed to understand the thermodynamic properties of QCD. A more comprehensive picture became available with the advent of first-principle calculations based on lattice QCD techniques. The quark-hadron phase transition is found to be a crossover by (2+1)-flavor lattice QCD calculations, in contrast to the SU(3) pure gauge case where a first-order phase transition is predicted Brown et al. (1990); Ali Khan et al. (2000); Aoki et al. (2006). Recent lattice QCD simulations with a physical pion mass have been able to provide us with the realistic equation of state at finite temperatures and vanishing chemical potentials Borsanyi et al. (2014); Bazavov et al. (2014).

At finite densities the equation of state is not well known, owing to the fermion sign problem of the first principle method (for a review see de Forcrand (2009)). Several techniques have been developed in lattice QCD, including the Taylor expansion method Gavai and Gupta (2001); Allton et al. (2002), the imaginary chemical potential method de Forcrand and Philipsen (2002); D’Elia and Lombardo (2003); Guenther et al. (2017), Lefschetz thimble decomposition Pham (1983); Witten (2011), and the complex Langevin method Parisi (1983); Klauder and Petersen (1985); Ambjorn and Yang (1985), but so far no complete calculations are available at larger chemical potentials. There can be non-trivial structures in the $\mu_{B}$-$T$ QCD phase diagram other than the QGP phase Fukushima and Hatsuda (2011); it is conjectured that there is a critical point based on the chiral model that predicts that the quark-hadron crossover becomes a first-order transition at some finite temperature and chemical potential Asakawa and Yazaki (1989). Beam energy scan (BES) programs are being performed to explore finite-density QCD matter and determine its detailed phase structure at RHIC and the CERN Super Proton Synchrotron (SPS), and are planned at several facilities including the GSI Facility for Antiproton and Ion Research (FAIR), JINR Nuclotron-based Ion Collider fAсility (NICA) and JAEA/KEK Japan Proton Accelerator Research Complex (J-PARC).

From the viewpoint of hydrodynamic modeling of relativistic nuclear collisions, the equation of state is needed to close the set of equations of motion, by characterizing the thermodynamic properties of the system. The equation of state at vanishing density obtained from lattice QCD calculations has been employed in comparisons of hydrodynamic simulations with experimental data from heavy-ion collisions Pratt et al. (2015); Sangaline and Pratt (2016); Bernhard et al. (2016); Pang et al. (2018); Monnai and Ollitrault (2017); Paquet et al. (2017), where bulk observables are generally well reproduced. For quantitative predictions and analysis of the BES experimental data, an equation of state at finite chemical potentials is needed as input to hydrodynamic models. Several works have been devoted to the construction of such quantitative models of the finite-density equation of state, including Refs. Nonaka and Asakawa (2005); Bluhm et al. (2005, 2007); Steinheimer et al. (2011); Huovinen and Petreczky (2011); Hempel et al. (2013); Albright et al. (2014, 2015); Rougemont et al. (2017); Critelli et al. (2017); Vovchenko et al. (2018); Parotto et al. (2018); Vovchenko et al. (2019); Fu et al. (2018a, b); Motornenko et al. (2019); Plumberg et al. (2018).

In this work, we present a framework to construct a QCD equation of state model (neos) with multiple charges: net baryon ($B$), strangeness ($S$) and electric charge ($Q$) based on state-of-the-art lattice QCD Borsanyi et al. (2012); Bellwied et al. (2015); Borsanyi et al. (2018); Bazavov et al. (2012a); Ding et al. (2015); Bazavov et al. (2017) and hadron resonance gas results. A version of the equation of state, which only has baryon chemical potential, has previously been introduced and used in Refs. Denicol et al. (2016); Monnai and Schenke (2016); Shen et al. (2017); Shen and Schenke (2018); Denicol et al. (2018); Shen and Schenke (2019); Gale et al. (2019). We numerically calculate the equation of state with conditions on the conserved charges close to those in relativistic heavy-ion collisions. This analysis is expected to be relevant in mid to low energy heavy-ion collisions and at forward rapidity where the fragments of shattered nuclei are relatively abundant Li and Shen (2018); Monnai (2012). The presence and interplay of the chemical potentials are expected to play an important role in for example the hadron chemistry or particle abundances (see e.g. Andronic et al. (2006)), and the separation of hadron and anti-hadron flow observables Karpenko et al. (2015); Hatta et al. (2015, 2016).

The multi-dimensional phase diagram has been studied Toublan and Kogut (2005); Xu et al. (2011); Kamikado et al. (2013); Ueda et al. (2013); Barducci et al. (2004); Nishida (2004); Son and Stephanov (2001) and some conjecture non-trivial phase structures. In this study, we consider a crossover equation of state as a baseline for the application to relativistic nuclear collisions.

The paper is organized as follows. In Sec. II, the construction of finite-density equations of state based on the Taylor expansion method for lattice QCD and the hadron resonance gas is presented. The numerical evaluation of the hybrid equation of state is performed in Sec. III, where the strangeness neutrality condition and the fixed charge-to-baryon ratio of nuclei are taken into account. In Sec. IV, particle ratios are estimated in hydrodynamic simulations assuming different conditions on the charge content of the system. Sec. V presents conclusions and discussions. Natural units $c=\hbar=k_{B}=1$ and the Minkowski metric $g^{\mu\nu}=\mathrm{diag}(+,-,-,-)$ are used.

II The equation of state

Based on the Taylor expansion method Gavai and Gupta (2001); Allton et al. (2002) we employ lattice QCD results of the conserved charge susceptibilities to construct the equation of state in the QGP phase. In the hadronic phase we use the equation of state of a hadron resonance gas, because the Taylor expansion method is not reliable at low temperatures. The use of a non-interacting resonance gas model is partly motivated by the good agreement between thermodynamic quantities at vanishing chemical potential, including susceptibilities, from lattice QCD and the hadron resonance gas. Also, the Cooper-Frye prescription Cooper and Frye (1974) of kinetic freeze-out requires that the hydrodynamic equation of state precisely matches that of the kinetic theory description of the hadron resonance gas on the freeze-out hypersurface for the successful conservation of energy-momentum and all charges. If at low temperatures the lattice result was used instead of the hadron resonance gas model, the truncation of the Taylor expansion at finite order would lead to an underestimation of the pressure in the hadronic phase, because higher order susceptibilities can be large for the hadron resonance gas in the Boltzmann limit, e.g., $\chi_{2n}^{B}/\chi_{2}^{B}=1$.

II.1 Lattice QCD equation of state

We consider the Taylor expansion method of lattice QCD as mentioned earlier. For the three-flavor QCD system, the expansion of the pressure around the vanishing density limit reads

[TABLE]

where $P$ is the pressure, $P_{0}$ is the pressure at vanishing chemical potentials, $T$ is the temperature, and $\mu_{B,Q,S}$ are the chemical potentials of baryon number, electric charge, and strangeness, respectively. $\chi_{l,m,n}^{B,Q,S}$ is the ($l$+$m$+$n$)-th order susceptibility defined at vanishing chemical potentials:

[TABLE]

The number $l$+$m$+$n$ should be even, owing to the matter-antimatter symmetry. The expansion is valid only when the fugacity $\mu_{B,Q,S}/T$ is sufficiently small. The lattice QCD results are parametrically extrapolated to high temperatures under the condition that they do not violate the Stefan-Boltzmann limits. See Appendix A.

II.2 Hadron resonance gas equation of state

The hydrostatic pressure of the hadron resonance gas can be written as

[TABLE]

where $i$ is the index for particle species, $g_{i}$ is the degeneracy, $m_{i}$ the particle’s mass, and $K_{2}(x)$ is the modified Bessel function of the second kind. The index $k$ describes the expansion of quantum distributions around the classical ones. It is generally sufficient to take into account the contributions of $k\leq 3$ for pions, $k\leq 2$ for kaons and $k=1$ for the heavier particles. The upper signs are for fermions and the lower signs for bosons. The hadronic chemical potential is $\mu_{i}$ = $B_{i}\mu_{B}+Q_{i}\mu_{Q}+S_{i}\mu_{S}$ where $B_{i}$, $Q_{i}$, and $S_{i}$ are the quantum numbers of net baryon, electric charge, and strangeness, respectively.

II.3 Hybrid equation of state

The complete nuclear equation of state is constructed by connecting the pressure of the lattice QCD equation of state to that of the hadron resonance gas model Huovinen and Petreczky (2010)

[TABLE]

where $J=\{B,Q,S\}$. Here $f(T,\mu_{J})$ is an arbitrary function for the connection of the two functions which satisfies $f\to 1$ when $T\gg T_{c}$ and $f\to 0$ when $T\ll T_{c}$, where $T_{c}$ is the connecting temperature. In this work we choose $f$ to be a hyperbolic tangent, defined in (14). A temperature shift $T_{s}(T,\mu_{J})$ is introduced phenomenologically to preserve the monotonicity conditions of thermodynamic variables at larger chemical potentials. Since $T_{s}$ is generally a function of $T$ and $\mu_{J}$, one can define

[TABLE]

which is the shifted QGP equation of state. Here $\tilde{P}_{\mathrm{lat}}$ should reduce to $P_{\mathrm{lat}}$ at small chemical potentials. While in this work we will not use a temperature shift, the shifting temperature is used in similar constructions of the equation of state in previous works  Denicol et al. (2016); Monnai and Schenke (2016); Shen et al. (2017); Shen and Schenke (2018); Denicol et al. (2018); Shen and Schenke (2019); Gale et al. (2019).

We require that the thermodynamic variables monotonously increase as functions of $T$ and $\mu_{J}$, respectively, as

[TABLE]

Those conditions may be trivially satisfied for the hadron resonance gas or lattice QCD equation of state when $\mu_{J}/T$ is not large, but the connection procedure can make it non-trivial. The conditions can be expressed as

[TABLE]

[TABLE]

Assuming that the thermodynamic quantities on the lattice QCD side are larger than those on the hadron resonance gas side, the conditions are still not trivially satisfied as the second order derivatives of $f$ can be negative. These conditions will be handled numerically in the next section.

The thermodynamic variables, the entropy density $s$, the conserved charge densities $n_{J}$, the energy density $e$, and the sound velocity $c_{s}$ are given as

[TABLE]

respectively, using the standard thermodynamic relations. Here $\{J,K\}=B,Q,S$ and $J\neq K$.

II.4 Multiple charges in nuclear collisions

A standard nucleus is made of protons and neutrons so the averaged density of strangeness is zero, which may be expressed as $n_{S}(T,\mu_{B},\mu_{Q},\mu_{S})=0$. However, neglecting electric charge for the moment, the conventional choice of $\mu_{S}=0$ leads to $n_{S}\neq 0$, because the strangeness density is dependent on $\mu_{B}$ as the strange quark carries both net baryon number and strangeness. Thus, in the presence of a finite net-baryon number, $\mu_{S}$ should generally be non-vanishing, so $n_{S}=0$ can be fulfilled. The condition, of course, can in principle be locally broken in the presence of geometrical fluctuations or diffusion processes. The equation of state with the strangeness neutrality condition can be expressed in terms of $T$ and $\mu_{B}$ because $\mu_{S}=\mu_{S}(T,\mu_{B})$.

The electric charge density $n_{Q}(T,\mu_{B},\mu_{Q},\mu_{S})$ is non-vanishing in nuclei as $n_{Q}=(Z/A)n_{B}$, where $Z$ is the proton number and $A$ is the nucleon number. The list of $Z/A$ ratios of the nuclei used in collider experiments is shown in Table 1. The typical ratio for heavy nuclei such as Au or Pb is $Z/A\approx 0.4$. The precise $n_{Q}$ dependence is expected to become more important when comparing collisions of isobar systems.

III Numerical construction

In this section, we numerically construct the hybrid QCD equation of state at finite densities (4). On the lattice QCD side, we make use of one of the latest (2+1)-flavor calculations of the equation of state at vanishing chemical potentials Bazavov et al. (2014) and the diagonal and off-diagonal susceptibilities up to the fourth order Bazavov et al. (2012a); Ding et al. (2015); Bazavov et al. (2017); Sharma . In addition, we introduce some of the most relevant sixth order susceptibilities, to allow for a proper matching of all quantities, as discussed later. The specific functional forms of parametrization are summarized in Appendix B.

All hadron resonances from the particle data group Tanabashi et al. (2018) with $u$, $d$, and $s$ constituent quark components and masses smaller than 2 GeV are taken into account in the resonance gas model. We use

[TABLE]

where $T_{c}(\mu_{B})=0.16\ \mathrm{GeV}-0.4\times(0.139\ \mathrm{GeV}^{-1}\mu_{B}^{2}+0.053\ \mathrm{GeV}^{-3}\mu_{B}^{4})$ motivated by and modified from the chemical freeze-out curve Cleymans et al. (2006). The dependence of the connecting temperature on strangeness and electric charge chemical potentials is assumed to be weak and neglected for the moment. The connecting width is $\Delta T_{c}=0.1T_{c}(0)$. Here, we choose not to perform a temperature shift and use $T_{s}=T$. It should be noted that this is not a unique choice of the parameters, but the thermodynamic monotonicity conditions (6) and (7) leave a rather narrow window for the possible parameter values. The smooth matching leads to an equation of state with a crossover transition. Implementation of a QCD critical point and the first-order phase transition is also possible for different choices of $f$. It will be discussed elsewhere as the location of the critical point and the critical behavior near it are currently not well known. For a possible approach to include a critical point see Critelli et al. (2017); Parotto et al. (2018); Plumberg et al. (2018).

The sixth-order susceptibilities should be relevant near and below the crossover transition. The term involving $\chi_{6}^{B}$ is naively expected to give the largest contribution to the pressure and the net baryon number because of the hierarchy in the chemical potentials $\mu_{B}>\mu_{S}>\mu_{Q}$ in nuclear collisions. The strangeness density and the electric charge density are not directly sensitive to $\chi_{6}^{B}$, because they are derivatives with respect to $\mu_{S}$ or $\mu_{Q}$, respectively, implying that the terms involving $\chi_{5,1}^{B,S}$ and $\chi_{5,1}^{B,Q}$ will be the important ones for them. We introduce those three susceptibilities in a phenomenological approach so that $n_{B}$, $\mu_{S}$, and $\mu_{Q}$ are smooth functions of $T$ and $\mu_{B}$, and that the results of the hadron resonance gas model are preserved below $T_{c}$, because of the relatively large uncertainties in the current lattice calculations of higher order susceptibilities. We find that the effects of the sixth-order susceptibilities are limited to the large chemical potential regions near the crossover transition.

For the strangeness and electric charges, we consider three cases: (i) $\mu_{S}=\mu_{Q}=0$, (ii) $n_{S}=0$ and $\mu_{Q}=0$, and (iii) $n_{S}=0$ and $n_{Q}=0.4n_{B}$. They are reffered to as neos B, neos BS and neos BQS, respectively. The first is the commonly used scenario in which one assumes that the net baryon chemical potential is the only non-vanishing one in the system. The second imposes the strangeness neutrality condition but neglects the electric charge chemical potential. The third is the most realistic scenario for the collision of heavy nuclei where $Z/A\sim 0.4$ Bazavov et al. (2017). It is also straightforward to calculate the equation of state as functions of $\mu_{B}$, $\mu_{Q}$, and $\mu_{S}$ for more general systems.

III.1 Vanishing strangeness and electric charge chemical potentials $\mu_{S}=\mu_{Q}=0$

First, the case where $\mu_{S}=\mu_{Q}=0$ is investigated. The pressure of the resulting equation of state is plotted in Fig. 1 (a). One can see the monotonous increase of $P$ as a function of $T$ or $\mu_{B}$. The equation of state reduces to that of lattice QCD at $\mu_{B}=0$ at the vanishing density limit. The constant entropy density over net baryon density lines, $s/n_{B}$ = 420, 144, 51, and 30, are plotted to illustrate the most relevant regions for the BES programs. They correspond to Au+Au collisions at $\sqrt{s_{NN}}=$ 200, 62.4, 19.6, and 14.5 AGeV, respectively Guenther et al. (2017). Note that the ratio is roughly constant during the time evolution in nuclear collisions when the entropy and the net baryon number are conserved, which is the case for the nearly-perfect fluid. The trajectory of $s/n_{B}$ is a straight line at higher temperatures where the system is relatively close to conformal, because $s/n_{B}\sim T/\mu_{B}$. It turns around near the crossover towards lower temperatures as pions begin to dominate over protons because of the mass difference, and large baryon chemical potential is required to have protons at lower temperatures for keeping the $s/n_{B}$ ratio fixed. In the limit $T\to 0$, the chemical potential approaches the proton mass.

As discussed earlier, the condition $\mu_{S}=0$, which is often assumed in nuclear collision analyses, leads to a non-vanishing strangeness density $n_{S}$. The value of $-n_{S}/T^{3}$ is shown in Fig. 1 (b). Positive baryon chemical potential leads to negative strangeness density because the $s$ quark has a negative strangeness chemical potential. The high temperature behavior can be understood as $n_{S}/T^{3}\sim\chi_{1,1}^{B,S}\mu_{B}/T\sim-\mu_{B}/3T$ (21). The strangeness density is relatively small at lower temperatures because kaons, the lightest strange hadrons, have net baryon number zero.

III.2 Strangeness neutrality $n_{S}=0$ and vanishing electric charge chemical potential $\mu_{Q}=0$

We next study the case where $n_{S}=0$ and $\mu_{Q}=0$. The pressure is shown as a function of $T$ and $\mu_{B}$ in Fig. 2 (a). The $n_{Q}/n_{B}$ ratio is arbitrary, and approaches $n_{Q}/n_{B}\sim 0.5$ in the parton gas limit (29). One can see that the equation of state is modified at larger baryon chemical potentials compared with that of the $\mu_{S}=0$ case. The constant $s/n_{B}$ lines are also shifted to larger $\mu_{B}$ (Fig. 2 (b)) because the strangeness neutrality implies $\mu_{S}\sim\mu_{B}/3$ at high temperatures. For clarity we show the projections of the constant $s/n_{B}$ lines onto the $\mu_{B}$-$T$ plane in Fig. 3. Here one can see that $\mu_{B}$ has to be about 3/2 times larger for a given $n_{B}$ when $\mu_{S}\neq 0$. The gray area in the figure shows a conjectured region $\mu_{B}/T>3$ where the Taylor expansion method of lattice QCD is not well-defined. A larger value of $\mu_{B}$ will result in a larger thermodynamic force $\nabla^{\mu}(\mu_{B}/T)$ for the net baryon diffusion current Denicol et al. (2018). It will have an important effect on constraining the net baryon diffusion constant in future phenomenological studies. While initial strangeness fluctuations and a strangeness diffusion current can break strangeness neutrality locally, this should not diminish the effect of enhanced $\mu_{B}$ on the net baryon diffusion, since it is a sub-leading effect and the strangeness is still globally conserved at zero.

Importantly, $\mu_{s}$ is non-zero at freeze-out, which will affect results on particle-antiparticle ratios of strange hadrons in hydrodynamic models, as is the case in thermal models Andronic et al. (2006). The potentially large effect of the strangeness neutrality condition is also discussed in Ref. Bazavov et al. (2012b).

It should be noted that while we have now imposed more realistic conditions compared to the previous case, the thermodynamic properties of the QCD system itself remain the same, we merely look at different slices of the multi-dimensional equation of state.

III.3 Strangeness neutrality $n_{S}=0$ and fixed electric charge-to-baryon ratio $n_{Q}=0.4n_{B}$

Finally, we investigate the case where $n_{S}=0$ and $n_{Q}=0.4n_{B}$. This is the setup most relevant to Au+Au and Pb+Pb collisions. The dimensionless pressure $P/T^{4}$ is plotted in Fig. 4 (a). The difference from the previous case is small in this setup but should be meaningful for correctly understanding particle-antiparticle ratios of charged particles.

The electric charge chemical potential shown in Fig. 4 (b) is negative, owing to the interplay of multiple conserved charges. Since the number of neutrons is larger than that of protons in heavy nuclei, $d$ quarks are slightly more abundant than $u$ quarks in the QGP phase and $\pi^{-}$ more abundant than $\pi^{+}$ in the hadronic phase. While the overall system is positively charged, a negative electric chemical potential is needed for describing this situation. $\mu_{Q}$ becomes positive for the system of 3He since $Z/A>1/2$. This would have to be taken into account for the collisions involving such nuclei.

It should be noted that $\mu_{Q}$ is small and is rather sensitive to the fine structure of the equation of state, including higher-order susceptibilities, at large chemical potentials. This implies that improvement in the lattice QCD calculations, including higher order susceptibilities, will be important in quantitative analyses.

III.4 Discussion

We have constructed the nuclear equation of state under several different conditions. We now study the differences between the different scenarios in more detail.

The sound velocities of the equation of state under the three different conditions along two constant $s/n_{B}$ lines are plotted in Fig. 5. One can see that finite-density effects are visible comparing the sound velocities of $s/n_{B}=420$ and $30$. Around the crossover temperature, the EoS becomes soft and $c_{s}^{2}$ has a minimum. The location of the minimum shifts towards lower temperatures as the net baryon density increases. Also, the sound velocity becomes larger in the QGP phase and smaller in the hadronic phase at larger chemical potentials. This is because the net baryon contribution in $c_{s}^{2}$ (13) is positive for the former phase and negative for the latter phase. At higher temperatures, it starts to approach the Stefan-Boltzmann limit $c_{s}^{2}=1/3$. For the three presented equations of state, $c_{s}^{2}$ reaches 94.8 % of the value the Stefan-Boltzmann limit at $T=0.6$ GeV and 97.2 % at $T=0.8$ GeV for $s/n_{B}=420$.

Comparing neos B to neos BS, the strangeness neutrality condition is found to slightly enhance the sound velocity in the QGP phase. It should be noted that if one neglected the derivatives involving $n_{S}$ in the calculation of $c_{s}^{2}$ (13) for neos B, the sound velocity would be further underestimated than our current result, because

[TABLE]

when $\mu_{S}=0$, i.e., $n_{S}\neq 0$, which again highlights the importance of adequate treatment of the multiple conserved charges. The effects of the fixed charge to baryon ratio on the sound velocity is almost negligible. Since the effect of the electric chemical potential is not large, the difference in the sound velocity is also not large when one neglects the derivatives involving $\mu_{Q}$ and $n_{Q}$ in the definition (13).

We next focus on the interplay of the multiple conserved charges and plot an arbitrary constant pressure plane in the chemical potential $\mu_{B}$-$\mu_{S}$-$\mu_{Q}$ space at a constant temperature in the hadronic phase in Figure 6 (a). This quantifies the chemical potential dependences of this thermodynamic quantity. For demonstration, we choose $P/T^{4}=0.8$ and $T=0.14$ GeV. The intercepts of each axis, defined implicitly as $P(\mu_{B}^{\mathrm{int}},0,0)=P(0,\mu_{Q}^{\mathrm{int}},0)=P(0,0,\mu_{S}^{\mathrm{int}})$, are ordered as $\mu_{B}^{\mathrm{int}}>\mu_{S}^{\mathrm{int}}>\mu_{Q}^{\mathrm{int}}$, partly reflecting the fact that the lightest particles that carry those charges are protons, kaons, and pions in the hadronic phase, respectively (3).

Figure 6 (b) presents the same in the QGP phase. Here $P/T^{4}=2$ and $T=0.2$ GeV are considered. The ordering of the intercepts can be seen to be $\mu_{B}^{\mathrm{int}}>\mu_{Q}^{\mathrm{int}}>\mu_{S}^{\mathrm{int}}$ in the QGP phase. This is consistent with the naïve expectation that $\mu_{B}/3\sim 2\mu_{Q}/3\sim\mu_{S}$ in the massless parton gas limit (16)-(19). The intercept $\mu_{S}^{\mathrm{int}}$ is slightly larger owing to the fact that it is still close to the crossover transition and that the strange quarks have a non-negligible mass. $\mu_{B}$ takes a maximum value at some positive finite $\mu_{S}$ because $s$ quarks have positive net baryon number and negative strangeness (19), i.e., the leading-order off-diagonal susceptibility between the net baryon and strangeness is negative. This is not the case for the cross-coupling between the electric charge and the net baryon or strangeness.

Constant $s/n_{B}$ lines for the case that $n_{S}=0$ and $n_{Q}=0.4n_{B}$ are plotted in Fig. 7 to illustrate the typical range of the chemical potentials covered by heavy-ion collider experiments. The trajectories coincide at high temperatures because $s/n_{B}\sim T/\mu_{B}$ and the conditions on $n_{S}$ and $n_{Q}$ make $\mu_{S}$ and $\mu_{Q}$ roughy proportional to $\mu_{B}$. The trajectories slightly bend towards the larger strangeness chemical potential at large $\mu_{B}$ and small $\mu_{S}$ regions, which correspond to low temperatures below $T\sim 0.1$ GeV. This could be caused by the suppression of kaons compared with pions owing to the mass difference. The behavior can also be seen in Fig. 2 (b). The bending does not occur for the charge chemical potential because pions, the lightest hadrons, have electric charge. It is note-worthy that broader ranges may be explored in actual collider events since the system is geometrically fluctuating and large local variation of the entropy-to-conserved-charge ratios can occur.

We note that in the region where the lattice QCD contribution dominates, the validity of our parametrization is limited to the range where $\mu_{B}/T$ is sufficiently small. From a practical point of view, for the application to nuclear collisions, however, these regions are not expected to much affect the bulk physics, because most of the fluid elements do not go through the large $\mu_{B}/T$ regions near $T_{c}$. This can be seen in the constant $s/n_{B}$ lines shown in Figs. 1 through 4.

It would also be interesting to compare our results with the ones obtained by other lattice QCD approaches to finite density regions, such as the one from the imaginary chemical potential method Vovchenko et al. (2017).

IV Application to nuclear collisions

To study the effects of imposing strangeness neutrality and a realistic charge-to-baryon ratio on observables in heavy ion collisions, we perform hybrid model calculations of Pb+Pb collisions at center of mass energy $\sqrt{s}=17.3\,{\rm AGeV}$ involving viscous hydrodynamic simulations with the three different equations of state described above and a microscopic hadronic afterburner. A more detailed description of the hybrid model is given in Schenke et al. (2019). We compute particle yields and compare particle ratios to experimental data from the Super Proton Synchrotron (SPS) Afanasiev et al. (2002); Alt et al. (2005, 2006, 2008a, 2008b, 2008c) (compiled in NA (4)).

For the initial state, we use the dynamical model presented in Shen and Schenke (2018), which provides fluctuating distributions of net baryon and energy-momentum densities in three spatial dimensions. The 3+1D hydrodynamic simulation Music Schenke et al. (2010, 2011, 2012) is run here with zero bulk viscosity and a constant shear viscosity to entropy density ratio of $\eta/s=0.08$. We switch from hydrodynamics to the hadron cascade UrQMD Bass et al. (1998); Bleicher et al. (1999) at a switching energy density $e_{\rm sw}$, whose value we vary below.

In Fig. 8 we show the particle yields (top) and particle ratios (bottom) from these simulations using $e_{\rm sw}=0.26\,{\rm GeV/fm}^{3}$, and the neos equations of state with different constraints on strangeness and electric charge. One can see that imposing strangeness neutrality has a visible effect - mainly on the strange and anti-strange particle yields. This effect is amplified in the particle ratios. Yields of particles with positive strangeness are increased while those of particles with negative strangeness are decreased, which is due to the finite positive strangeness chemical potential present in neos BS (and neos BQS).

The agreement between the theoretical calculations and experimental data is improved for most particles with strangeness in neos BS. In the meson sector where $\mu_{B}$ is absent, the ratio of $K^{+}$ over $K^{-}$ gets enhanced by the strangeness neutrality condition and agreement with experimental data at SPS energy improves.

Protons, and to a lesser degree anti-protons, are modified, because in the presence of $\mu_{S}$, the baryon chemical potential $\mu_{B}$ also changes. The small change for pions is likely due to the modification of resonance abundances when going from neos B to neos BS.

Introducing the constraint on the electric charge by using neos BQS, we find only very mild modifications of the particle yields. The negative $\mu_{Q}$ leads to a slight increase of negative relative to positive charged particles, as can be best seen in the plot of the particle ratios as the difference between the points for neos BS and neos BQS. The introduction of $\mu_{Q}$ can explain at least qualitatively that $\pi^{-}$ are slightly more abundant than $\pi^{+}$. The ratio $\bar{\Omega}/\Omega$ behaves in the opposite way. Possibly changes of $\mu_{B}$ and $\mu_{S}$ when introducing $\mu_{Q}$ could contribute to this behavior. We note that the main effect of baryon–anti-baryon annihilation within UrQMD is the reduction of the anti-proton yield by approximately 35% at $\sqrt{s_{\mathrm{NN}}}=17.3$ GeV. Yields of $\bar{\Lambda}$ and $\bar{\Omega}$ are reduced by 25% and 20%, respectively.

In Fig. 9 we study the effect of different switching energy densities on particle yields (top) and ratios (bottom) for neos BQS. Using a lower switching energy density means assuming that the system is fully chemically equilibrated down to lower temperatures. Anti-baryons are most sensitive to the switching energy density. The reason could be that lower $e_{\rm sw}$ means lower temperature at switching, which goes along with a larger baryon chemical potential (see Fig. 3). Lower temperature tends to decrease heavier particles’ yields, while the larger baryon chemical potential will lead to more baryons, weakening the effect of lower temperature, and to less anti-baryons, adding to the effect. We find that the experimental data on particle ratios prefers a switching temperature between $e_{\rm sw}=0.16$ and $0.36\,{\rm GeV/fm}^{3}$.

V Conclusions

By matching the Taylor expanded lattice QCD equation of state in the high-temperature region to that of a hadron resonance gas model at low temperature, we have constructed the QCD equation of state at finite net baryon, electric charge, and strangeness chemical potentials. We employ phenomenologically motivated sixth order susceptibilities to allow for a smooth matching that respects the thermodynamic monotonicity conditions. The equation of state is designed to be used in simulations of collisions of heavy nuclei such as Au or Pb in a wide range of collision energies explored by the beam energy scan programs.

The equation of state without strangeness chemical potential leads to the violation of the strangeness neutrality condition, which should hold in heavy ion collisions. A positive finite strangeness chemical potential is observed when the condition $n_{S}=0$ is kept. Also enforcing the approximate relation between electric charge and baryon number for heavy nuclei, $n_{Q}/n_{B}=0.4$, introduces a negative finite electric chemical potential. These constraints should be important for understanding the difference between particle yields and flow observables of particles and anti-particles within hydrodynamic models of heavy-ion collisions.

We have presented the effect of enforcing different constraints on strangeness and electric charge on the particle yields and ratios in $\sqrt{s}=17.3\,{\rm AGeV}$ Pb+Pb collisions. While strange and anti-strange particles are most affected by these constraints, modifications of non-strange particles were also observed, mostly driven by the finite $\mu_{S}$, less so the finite, negative, $\mu_{Q}$. This is understood by the fact that the introduction of $\mu_{S}$ (and $\mu_{Q}$) will also alter $\mu_{B}$.

These results are also important because they imply that we do not explore the $\mu_{B}$-$T$ plane in the BES experiments, but a certain slice in the $\mu_{B}$-$\mu_{Q}$-$\mu_{S}$-$T$ hyper-plane. This can affect the search of the QCD critical point because the traditional critical point at $\mu_{Q}=\mu_{S}=0$ may not be reached.

Further importance may arise with regard to isobar collisions. The equation of state can be different for different isobar collision systems, which should be taken into account for correctly understanding the background signals to the chiral magnetic effect. A similar discussion is applicable to small systems, where light nuclei such as proton, deuteron, or 3He tend to have larger $Z/A$ ratios, and the sign of the electric charge chemical potential can be flipped. It is also possible to perform event-by-event hydrodynamic analyses distinguishing protons and neutrons.

Future prospects for model improvements include introduction of the full sixth order susceptibilities from lattice QCD calculations. Although they are vanishing in the high temperature limit, they could play a non-trivial role near the crossover at larger chemical potentials.

Our equation of state model neos is publicly available neo .

Acknowledgements.

The authors thank Frithjof Karsch, Swagato Mukherjee, and Sayantan Sharma for useful discussion. AM is supported by JSPS KAKENHI Grant Number JP19K14722. BPS is supported under DOE Contract No. DE-SC0012704. CS is supported under DOE Contract No. DE-SC0013460. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This work is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration.

Appendix A The Stefan-Boltzmann limit

For the massless non-interacting system of $u$, $d$ and $s$, the analytical expression of the pressure reads

[TABLE]

where $N_{f}=3$ is the number of flavors. The speed of sound (13) is independent of the temperature and chemical potentials, $c_{s}^{2}=1/3$. It is note-worthy that the sixth and higher order susceptibilities are vanishing in this case. The susceptibilities of the conserved charges are given using Eqs. (1), (16) and

[TABLE]

The second-order susceptibilities in the Stefan-Boltzmann limit are

[TABLE]

and the fourth-order ones are

[TABLE]

They are used as anchors for the high temperature behaviors of the equation of state, where lattice QCD data are scarce, so that the basic thermodynamic features are not violated in the large $T$ limit.

One can analytically solve the linearized equations for charge densities $n_{B}$, $n_{Q}$, and $n_{S}$ obtained by keeping the second order diagonal and off-diagonal susceptibilities:

[TABLE]

The solutions are

[TABLE]

in the Stefan-Boltzmann limit. In the case where $n_{S}=0$ and $n_{Q}=0.4n_{B}$, those can be expressed as $\mu_{B}=4.6n_{B}/T^{2}$, $\mu_{Q}=-0.2n_{B}/T^{2}$, and $\mu_{S}=1.6n_{B}/T^{2}$. It is worth mentioning that the sign of $\mu_{Q}$ is rather sensitive to the proportionality constant between the net baryon and charge densities and turns positive at $n_{Q}=0.5n_{B}$.

Appendix B Parametrizations of susceptibilities

The parametrizations of the diagonal and off-diagonal susceptibilities at zero chemical potentials in the regime above $T_{c}$ are presented here. The second-order susceptibilities are parametrized as

[TABLE]

where

[TABLE]

The parameters are listed in Table 2.

The functional forms for the fourth-order susceptibilities are

[TABLE]

where the base $\chi_{2}$ is chosen for purely parametric purposes. The individual parameters can be found in Table 3. The Stefan-Boltzmann limits are used to regulate the high temperature behavior of the parametrizations. It is note-worthy that the $h_{3}$ values are typically not the exact Stefan-Boltzmann values because they are parameters for the fitting which is valid conservatively up to around 600 MeV, and the convergence of the fitting functions can be rather slow. Also, the lattice data itself does not approach the limit at around $3T_{c}$ in some cases. The fitting can be further improved when lattice QCD data become available for a wider temperature range.

The sixth-order susceptibilities used are

[TABLE]

where the parameter coefficients are listed in Table 4. Unlike the second- and fourth order ones, they are not based on lattice QCD but determined phenomenologically from the thermodynamic conditions as mentioned in the main text. It should thus be noted that they are effectively contaminated by the contributions of higher order susceptibilities and can be different from those obtained in lattice calculations.