The QCD Equation of State to $\mathcal{O}(\mu_B^6)$ from Lattice QCD

TL;DR

None

Contribution

None

Abstract

We calculated the QCD equation of state using Taylor expansions that include contributions from up to sixth order in the baryon, strangeness and electric charge chemical potentials. Calculations have been performed with the Highly Improved Staggered Quark action in the temperature range $T\in [135~{\rm MeV}, 330~{\rm MeV}]$ using up to four different sets of lattice cut-offs corresponding to lattices of size $N_\sigma^3\times N_\tau$ with aspect ratio $N_\sigma/N_\tau=4$ and $N_\tau =6-16$. The strange quark mass is tuned to its physical value and we use two strange to light quark mass ratios $m_s/m_l=20$ and $27$, which in the continuum limit correspond to a pion mass of about $160$ MeV and $140$ MeV espectively. Sixth-order results for Taylor expansion coefficients are used to estimate truncation errors of the fourth-order expansion. We show that truncation errors are small for baryon chemical potentials less then twice the temperature ($\mu_B\le 2T$). The fourth-order equation of state thus is suitable for the modeling of dense matter created in heavy ion collisions with center-of-mass energies down to $\sqrt{s_{NN}}\sim 12$ GeV. We provide a parametrization of basic thermodynamic quantities that can be readily used in hydrodynamic simulation codes. The results on up to sixth order expansion coefficients of bulk thermodynamics are used for the calculation of lines of constant pressure, energy and entropy densities in the $T$-$\mu_B$ plane and are compared with the crossover line for the QCD chiral transition as well as with experimental results on freeze-out parameters in heavy ion collisions. These coefficients also provide estimates for the location of a possible critical point. We argue that results on sixth order expansion coefficients disfavor the existence of a critical point in the QCD phase diagram for $\mu_B/T\le 2$ and $T/T_c(\mu_B=0) > 0.9$.

Figures

Tables

Table 1. Table 1: Parameters used in the ansatz given in Eq.  37 for the interpolation of the expansion coefficients of the net baryon-number density of (2+1)-flavor QCD with vanishing net strangeness and a fixed ratio of electric charge and net baryon-number density, n Q / n B = 0.4 subscript 𝑛 𝑄 subscript 𝑛 𝐵 0.4 n_{Q}/n_{B}=0.4 . These interpolations have been determined for the temperature interval T ∈ [ 130 ​ MeV , 280 ​ MeV ] 𝑇 130 MeV 280 MeV T\in[130\leavevmode\nobreak\ {\rm MeV},280\leavevmode\nobreak\ {\rm MeV}] . Also given are parameters needed for the interpolation of the expansion coefficients for the electric charge chemical potential (Eq.  38 ) and the coefficients for the parametrization of the pressure at μ B = 0 subscript 𝜇 𝐵 0 \mu_{B}=0 given in Eq.  39 .

$k$	$N_{k,0n}^B$	$N_{k,1n}^B$	$N_{k,2n}^B$	$N_{k,3n}^B$	$N_{k,4n}^B$	$N_{k,1d}^B$	$N_{k,2d}^B$	$N_{k,3d}^B$	$N_{k,4d}^B$
1	0.302182	-0.929305	1.230560	-0.798724	0.204722	-2.011836	1.190147	0.003869	-0.076244
3	0.000446650	0.00983742	-0.0315076	0.0323632	-0.0107642	-1.327047	0.0472047	0.0	0.323696
5	0.0000104211	-0.000327321	0.00122751	-0.00158725	0.000672708	-1.467875	-0.264770	0.796010	-0.044968
$k$	$q_{k,0n}$	$q_{k,1n}$	$q_{k,2n}$	$q_{k,3n}$	$q_{k,4n}$	$q_{k,1d}$	$q_{k,2d}$	$q_{k,3d}$	$q_{k,4d}$
1	-0.114472	-0.631833	2.102001	-2.165174	0.739905	16.565265	-35.328733	19.940335	0.384797
3	0.0505332	-0.312052	0.700958	-0.662171	0.219351	-23.224117	82.688725	-89.160400	31.381036
5	0.0000842	-0.0005250	0.00113467	-0.00103897	0.00034414	-2.095094	0.987940	0.146830	-0.0210650
	$p_{0n}$	$p_{1n}$	$p_{2n}$	$p_{3n}$	$p_{4n}$	$p_{1d}$	$p_{2d}$	$p_{3d}$	$p_{4d}$
0	0.00556035	128.702341	-293.064074	228.763685	-58.084225	12.713331	0.0	-31.330957	26.524394

Table 2. Table 2: Pressure, energy density and entropy density, characterizing lines of constant physics which correspond to the conditions met for μ B = 0 subscript 𝜇 𝐵 0 \mu_{B}=0 at T 0 = 145 subscript 𝑇 0 145 T_{0}=145  MeV, 155 155 155  MeV and 165 165 165  MeV. Columns 2-4 give results in appropriate units of temperature, while columns 5-7 give the same results expressed in units of G ​ e ​ V 𝐺 𝑒 𝑉 GeV and f ​ m 𝑓 𝑚 fm .

	at ${\mu }_B=0$			on LCP
$T_0\left[\mathrm{MeV}\right]$	$p/T_0^4$	$ϵ/T_0^4$	$s/T_0^3$	$p\left[\mathrm{GeV}/{\mathrm{fm}}^3\right]$	$ϵ\left[\mathrm{GeV}/{\mathrm{fm}}^3\right]$	$s\left[{\mathrm{fm}}^{-3}\right]$
145	0.586(80)	3.52(47)	4.11(53)	0.0337(46)	0.203(27)	1.63(21)
155	0.726(95)	4.61(55)	5.34(63)	0.0546(71)	0.346(41)	2.59(30)
165	0.898(110)	5.76(59)	6.66(69)	0.0868(106)	0.556(57)	3.90(40)

Table 3. Table 3: Ensemble parameters for calculations with light to strange quark mass ratio m l / m s = 1 / 20 subscript 𝑚 𝑙 subscript 𝑚 𝑠 1 20 m_{l}/m_{s}=1/20 on lattices of size N σ 3 ​ N τ superscript subscript 𝑁 𝜎 3 subscript 𝑁 𝜏 N_{\sigma}^{3}N_{\tau} with N τ = 6 , 8 subscript 𝑁 𝜏 6 8 N_{\tau}=6,\ 8 and N σ = 4 ​ N τ subscript 𝑁 𝜎 4 subscript 𝑁 𝜏 N_{\sigma}=4N_{\tau} . Columns 4 and 8 give the number of gauge field configurations, separated by 10 RHMC steps, that contributed to the analysis of up to sixth order generalized susceptibilities χ i ​ j ​ k B ​ Q ​ S superscript subscript 𝜒 𝑖 𝑗 𝑘 𝐵 𝑄 𝑆 \chi_{ijk}^{BQS} .

$N_{\tau }=6$				$N_{\tau }=8$
$\beta$	$m_l$	T[MeV]	#conf.	$\beta$	$m_l$	T[MeV]	#conf.
6.245	0.00415	179.52	14521	6.515	0.00302	178.36	16933
6.341	0.00370	198.61	3745	6.550	0.00291	184.84	15853
6.423	0.00335	216.33	1481	6.575	0.00282	189.58	11853
6.515	0.00302	237.81	1408	6.608	0.00271	196.01	16760
6.664	0.00257	276.43	1364	6.664	0.00257	207.32	8358
				6.800	0.00224	237.07	5816
				6.950	0.00193	273.88	9550
				7.150	0.00160	330.23	9184

Table 4. Table 4: Same as Table  3 but for the light to strange quark mass ratio m l / m s = 1 / 27 subscript 𝑚 𝑙 subscript 𝑚 𝑠 1 27 m_{l}/m_{s}=1/27 and including results for N τ = 12 subscript 𝑁 𝜏 12 N_{\tau}=12 .

$N_{\tau }=6$				$N_{\tau }=8$				$N_{\tau }=12$
$\beta$	$m_l$	T[MeV]	#conf.	$\beta$	$m_l$	T[MeV]	#conf.	$\beta$	$m_l$	T[MeV]	#conf.
5.980	0.00435	135.29	81200	6.245	0.00307	134.64	180320	6.640	0.00196	134.94	5834
6.010	0.00416	139.71	120790	6.285	0.00293	140.45	172110	6.680	0.00187	140.44	5833
6.045	0.00397	145.05	120770	6.315	0.00281	144.95	138150	6.712	0.00181	144.97	13846
6.080	0.00387	150.59	79390	6.354	0.00270	151.00	107510	6.754	0.00173	151.10	14200
6.120	0.00359	157.17	66180	6.390	0.00257	156.78	135730	6.794	0.00167	157.13	15476
6.150	0.00345	162.28	79660	6.423	0.00248	162.25	115850	6.825	0.00161	161.94	16772
6.170	0.00336	165.98	49760	6.445	0.00241	165.98	120270	6.850	0.00157	165.91	19542
6.200	0.00324	171.15	122700	6.474	0.00234	171.02	139980	6.880	0.00153	170.77	21220
6.225	0.00314	175.76	122730	6.500	0.00228	175.64	133070	6.910	0.00148	175.76	12303

Equations

$$P(T,{\mu }_B)$$

$$ϵ(T,{\mu }_B)$$

$$(\Delta (P/T^4){)}_{2n}\equiv \frac{(P_B(T,{\mu }_B)-P_B(T,0){)}_{2n}}{T^4}$$

$$\frac{n_B}{T^3}=\frac{P_B(T,0)}{T^4}\sinh {\hat{\mu }}_B=\sum _{k=1}^{\infty }\frac{{\chi }_{2k}^{B,HRG}(T)}{(2k-1)!}{\hat{\mu }}_B^{2k-1}\simeq \frac{P_B(T,0)}{T^4}\sum _{k=1}^{\infty }\frac{1}{(2k-1)!}{\hat{\mu }}_B^{2k-1}.$$

$$\frac{P_{ideal}}{T^4}=\frac{8{\pi }^2}{45}+\sum _{f=u,d,s}\left[\frac{7{\pi }^2}{60}+\frac{1}{2}{\left(\frac{{\mu }_f}{T}\right)}^2+\frac{1}{4{\pi }^2}{\left(\frac{{\mu }_f}{T}\right)}^4\right],$$

$$\frac{P}{T^4}=\frac{1}{V{T}^3}\ln \mathcal{Z}(T,V,{\hat{\mu }}_u,{\hat{\mu }}_d,{\hat{\mu }}_s)=\sum _{i,j,k=0}^{\infty }\frac{{\chi }_{ijk}^{BQS}}{i!j!k!}{\hat{\mu }}_B^i{\hat{\mu }}_Q^j{\hat{\mu }}_S^k,$$

$${\mu }_u$$

$${\mu }_d$$

$${\mu }_s$$

$${\chi }_{ijk}^{BQS}\equiv {\chi }_{ijk}^{BQS}(T)={\frac{\partial P(T,\hat{\mu })/T^4}{\partial {\hat{\mu }}_B^i\partial {\hat{\mu }}_Q^j\partial {\hat{\mu }}_S^k}|}_{\hat{\mu }=0}.$$

$$\frac{n_X}{T^3}=\frac{\partial P/T^4}{\partial {\hat{\mu }}_X},X=B,Q,S.$$

$$\Delta (T,{\hat{\mu }}_B,{\hat{\mu }}_Q,{\hat{\mu }}_S)\equiv \frac{ϵ-3P}{T^4}=T\frac{\partial P/T^4}{\partial T}=\sum _{i,j,k=0}^{\infty }\frac{{\Xi }_{ijk}^{BQS}}{i!j!k!}{\hat{\mu }}_B^i{\hat{\mu }}_Q^j{\hat{\mu }}_S^k,$$

$${\Xi }_{ijk}^{BQS}(T)=T\frac{\mathrm{d}{\chi }_{ijk}^{BQS}(T)}{\mathrm{d}T}.$$

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

$$\frac{s}{T^3}$$

$${\hat{\mu }}_Q(T,{\mu }_B)$$

$${\hat{\mu }}_S(T,{\mu }_B)$$

$$\frac{P(T,{\mu }_B)}{T^4}-\frac{P(T,0)}{T^4}$$

$$\frac{n_X}{T^3}$$

$${\hat{\mu }}_B\frac{\mathrm{d}P/T^4}{\mathrm{d}{\hat{\mu }}_B}={\hat{\mu }}_B\frac{n_B}{T^3}+{\hat{\mu }}_B\frac{\mathrm{d}{\hat{\mu }}_Q}{\mathrm{d}{\hat{\mu }}_B}\frac{n_Q}{T^3}+{\hat{\mu }}_B\frac{\mathrm{d}{\hat{\mu }}_S}{\mathrm{d}{\hat{\mu }}_B}\frac{n_S}{T^3},$$

$$P_{2n}=\frac{1}{2n}\left(N_{2n-1}^B+\sum _{k=1}^n(2k-1)\left(s_{2k-1}{N}_{2n-2k+1}^S+q_{2k-1}{N}_{2n-2k+1}^Q\right)\right).$$

$$T\frac{\mathrm{d}P/T^4}{\mathrm{d}T}$$

$$\Delta (T,{\hat{\mu }}_B)=\frac{ϵ-3P}{T^4}={\left(\frac{ϵ-3P}{T^4}\right)}_{{\hat{\mu }}_B=0}+\sum _{n=1}^{\infty }(T{P}_{2n}^{\prime }(T)-h_{2n}(T)){\hat{\mu }}_B^{2n},$$

$$h_{2n}=\sum _{k=1}^n\left(s_{2k-1}^{\prime }{N}_{2n-2k+1}^S+q_{2k-1}^{\prime }{N}_{2n-2k+1}^Q\right).$$

$$t_{2n}=\sum _{k=1}^n\left(s_{2k-1}{N}_{2n-2k+1}^S+q_{2k-1}{N}_{2n-2k+1}^Q\right).$$

$$\frac{ϵ(T,{\mu }_B)}{T^4}-\frac{ϵ(T,0)}{T^4}$$

$$\frac{s(T,{\mu }_B)}{T^3}-\frac{s(T,0)}{T^3}$$

$$\mathrm{Tr}{M}_f^{-1}{M}_f^{\prime }{M}_f^{-1}{M}_f^{\prime }....M_f^{-1}{M}_f^{\prime },$$

$$P_{2n}=\frac{1}{2n}{N}_{2n-1}^B=\frac{1}{(2n)!}{\chi }_{2n}^B.$$

$$\frac{P(T,{\mu }_B)-P(T,0)}{T^4}$$

$$\frac{n_B}{T^3}$$

$$\frac{n_X}{n_B}=\frac{{\chi }_{11}^{BX}+\frac{1}{6}{\chi }_{31}^{BX}{\hat{\mu }}_B^2+\frac{1}{120}{\chi }_{51}^{BX}{\hat{\mu }}_B^4}{{\chi }_2^B+\frac{1}{6}{\chi }_4^B{\hat{\mu }}_B^2+\frac{1}{120}{\chi }_6^B{\hat{\mu }}_B^4},X=Q,S.$$

$$\Delta \left(ϵ/T^4\right)$$

$$\Delta \left(s/T^3\right)$$

Full text

March 18, 2024

I Introduction

The temperature and density dependence of bulk thermodynamic quantities, commonly summarized as the equation of state (EoS), provide the most basic characterization of equilibrium properties of strong-interaction matter. Its analysis within the framework of lattice regularized Quantum Chromodynamics (QCD) has been refined ever since the early calculations performed in pure $SU(N)$ gauge theories Engels:1980ty . Quite recently, the continuum extrapolated results for the EoS of QCD with physical light and strange quark masses have been calculated Borsanyi:2013bia ; Bazavov:2014pvz . Bulk thermodynamic observables such as pressure ($P$), energy density ($\epsilon$) and entropy density ($s$) as well as second order quantities such as the specific heat ($C_{V}$) and velocity of sound ($c_{s}$) have now been obtained at vanishing chemical potentials for the three quark flavors $(\mu_{u},\mu_{d},\mu_{s})$. In accordance with the analysis of the chiral transition temperature, $T_{c}\simeq(154\pm 9)$ MeV Bazavov:2011nk , bulk thermodynamic observables change smoothly in the transition region. At low temperature they are found to be in quite good agreement with hadron resonance gas (HRG) model calculations, although some systematic deviations have been observed, which may be attributed to the existence of additional resonances which are not taken into account in HRG model calculations based on well established resonances listed in the particle data tables Majumder:2010ik ; Bazavov:2014xya .

The EoS at vanishing chemical potentials does already provide important input into the modelling of the hydrodynamic evolution of hot and dense matter created in heavy ion collisions. While this is appropriate for the thermal conditions met in these collisions at the LHC and the highest RHIC beam energies, knowledge of the EoS at non-vanishing baryon ($\mu_{B}$), strangeness ($\mu_{S}$) and electric charge ($\mu_{Q}$) chemical potentials is indispensable for the hydrodynamic modelling of the conditions met in the beam energy scan (BES) at RHIC. Due to the well-known sign problem for lattice QCD formulations at non-zero chemical potential a direct calculation of the EoS at non-zero $(\mu_{B},\ \mu_{Q},\ \mu_{S})$ is unfortunately not yet possible. At least for small values of the chemical potentials this can be circumvented by using a Taylor expansion of the thermodynamic potential Gavai:2001fr ; Allton:2002zi . In this way some results for EoS at non-zero baryon chemical potential have been obtained on coarse lattices Allton:2002zi ; Gavai:2003mf ; Allton:2003vx . These calculations have even been extended to sixth order in the baryon chemical potential Allton:2005gk ; Ejiri:2005uv . First continuum extrapolated results for the EoS using second order Taylor expansion coefficients have been obtained within the stout discretization scheme for staggered fermions Borsanyi:2012cr and simulations at imaginary chemical potential have been used to arrive at a sixth order result for the QCD EoS Gunther:2016vcp and up to eighth order for some generalized susceptibilities DElia:2016jqh through analytic continuation.

Results for higher order expansion coefficients are clearly needed if one wants to cover the range of chemical potentials, 0\leq\mu_{B}/T\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\ 3 that is expected to be explored with the BES at RHIC by varying the beam energies in the range $7.7\leavevmode\nobreak\ {\rm GeV}\leq\sqrt{s_{NN}}\leq 200\leavevmode\nobreak\ {\rm GeV}$. Of course, the Taylor expansions will break down, should the elusive critical point in the QCD phase diagram CEP ; Stephanov turn out to be present in this range of baryon chemical potentials. The convergence of the series thus needs to be monitored carefully.

This paper is organized as follows. In the next section we briefly discuss Taylor series for a HRG model in Boltzmann approximation. This helps to argue for the significance of sixth order Taylor expansions. In Section  III, we present the basic framework of Taylor series expansions, introduce expansions in the presence of global constraints and discuss some details of our calculations and the ensembles used. In Section IV we discuss the $6^{th}$ order Taylor expansion of QCD thermodynamics in the simplified case of vanishing strangeness and electric charge chemical potentials. Section V is devoted to the corresponding discussion of strangeness neutral systems $n_{S}=0$ with fixed net electric charge ($n_{Q}$) to net baryon-number ($n_{B}$) ratio, which is of relevance for the description of hot and dense matter formed in heavy ion collisions where typically $n_{Q}/n_{B}\simeq 0.4$. We discuss the relevance of a non-vanishing electric charge chemical potential by considering electric charge neutral ($n_{Q}/n_{B}=0$) as well as isospin symmetric ($n_{Q}/n_{B}=1/2$) systems. At the end of this section we present a parametrization of the equation of state that can easily be used as input for the modeling of the thermal conditions met in heavy ion collisions. In Section VI we present results on lines of constant pressure, energy density and entropy density and compare their dependence on $\mu_{B}$ with empirical results for the freeze-out conditions observed in heavy ion collisions. We comment on the radius of convergence of the Taylor series for the pressure and resulting constraints for the location of a possible critical point in Section VII. Finally we present our conclusions in Section  VIII. Details on (A) the statistics and simulation parameters, (B) explicit expressions for the expansions of electric charge and baryon number chemical potentials, and (C) explicit expressions for the expansion parameters of the lines of constant physics are given in three Appendices A-C.

II Taylor expansions and the low and high temperature limits of strong interaction matter

The main aim of this work is to supply a EoS of strong-interaction matter using up to sixth order Taylor expansions for bulk thermodynamic observables. As we will see later at present results on sixth order expansion coefficients in the Taylor series will mainly help to constrain truncation errors in the fourth order expansion rather than providing accurate results on the sixth order contribution to thermodynamic quantities. We will argue that our analysis provides reliable results for the EoS for baryon chemical potentials up to $\mu_{B}/T\simeq 2$ at temperatures below $T\simeq 160$ MeV and for an even larger range in $\mu_{B}/T$ at higher temperatures.

Before turning to a discussion of lattice QCD results on the EoS, it may be useful to analyze truncation effects in the hadron resonance gas (HRG) model, which seems to provide a good approximation for thermodynamics in the low temperature, hadronic regime. For simplicity let us consider the case of vanishing electric charge and strangeness chemical potentials, $\mu_{Q}=\mu_{S}=0$. At temperatures close to the transition temperature $T_{c}\simeq 154$ MeV and for baryon chemical potentials less than a few times the transition temperature, the baryon sector of a HRG is well described in the Boltzmann approximation. In a HRG model calculation based on non-interacting hadrons the pressure may then be written as

[TABLE]

where we introduced the notation $\hat{\mu}_{B}\equiv\mu_{B}/T$ and $P_{M}(T)$ ($P_{B}(T,\hat{\mu}_{B})$) denote the meson (baryon) contributions to the pressure. A similar relation holds for the energy density,

[TABLE]

with $\epsilon_{M/B}\equiv T^{2}\left(\partial(P_{M/B}/T)/\partial T\right)_{\hat{\mu}_{B}}$. The $\mu_{B}$-dependent contribution thus is simple and can easily be represented by a Taylor series. Truncating this expansion at $(2n)$-th order we obtain

[TABLE]

where in the last equality we made use of the fact that in HRG models constructed from non-interacting, point-like hadrons, all expansion coefficients are identical when using a Boltzmann approximation for the baryon sector, i.e. all baryon number susceptibilities are identical, $\chi_{2k}^{B,HRG}=P_{B}(T,0)$. The ratios of these susceptibilities are unity, $\chi_{2k}^{B,HRG}/\chi_{2(k-1)}^{B,HRG}=\chi_{2k}^{B,HRG}/\chi_{2}^{B,HRG}=1$. Similarly one finds for the net baryon-number density,

[TABLE]

Higher order corrections are thus more important in the net baryon-number density than in the expansions of the pressure or energy density. For instance, the contribution to $\mu_{B}n_{B}/T^{4}$ at ${\cal O}(\hat{\mu}_{B}^{2k})$ is a factor $2k$ larger than the corresponding ${\cal O}(\hat{\mu}_{B}^{2k})$ expansion coefficient of the pressure.

In Fig. 1 we show results from a Taylor series expansion of the $\mu_{B}$-dependent part of the pressure in a HRG model truncated after leading order (LO), next-to-leading order (NLO) and next-to-next-to-leading (NNLO) order. These truncated expansions are compared to the exact result, i.e. $\left(\Delta P\right)_{\infty}(T)=P_{B}(T)(\cosh(\hat{\mu}_{B})-1)$. The insertion shows the deviation of the $n$-th order truncated Taylor series ($\left(\Delta P\right)_{n}(T)$) from the exact result ($\left(\Delta P\right)_{\infty}(T)$). As can be seen already the fourth order Taylor series provides a good approximation for the pressure (and energy as well as entropy density) of a HRG for all $\mu_{B}\leq 2T$. At $\mu_{B}=2T$ the fourth order Taylor series for the $\mu_{B}$-dependent contribution to the pressure deviates by less than 5% from the exact result. These deviations are, of course, even smaller in the total pressure which in the temperature range of interest is dominated by the meson contribution. Even at $T=170$ MeV, which certainly is already above the range of applicability of HRG models, the baryonic contribution to the pressure (energy density) amounts only to about 20% (30%). A 5% truncation error in the $\mu_{B}$-dependent contribution to the pressure or energy density thus amounts to less than a 2% effect in the total pressure or energy density. Similar estimates hold for the more general case of non-vanishing $\mu_{Q}$ and $\mu_{S}$.

Of course, the good convergence properties of the Taylor series for the pressure in HRG models also reflect that the radius of convergence of this series is infinite. If there exists a critical point in the QCD phase diagram one cannot expect to find that the Taylor series is that well behaved. Still the HRG result provides a benchmark also for the QCD case. If the radius of convergence of the Taylor series for the QCD pressure is finite and, in particular, smaller than $\mu_{B}\simeq 3T$, one should find large deviations in the generalized susceptibilities from the corresponding HRG results. Ratios of susceptibilities have to grow asymptotically like, $\chi_{2k}^{B,QCD}/\chi_{2(k-1)}^{B,QCD}\sim k^{2}$ in order to yield a finite radius of convergence for a Taylor expansion. We will come back to a discussion of this asymptotic behavior after having discussed our sixth order calculation of Taylor expansion coefficients.

Let us briefly mention also the high temperature limit. At large values of the temperature, the pressure approaches that of a massless ideal gas of quarks and gluons. In this limit the pressure is just a second order polynomial in $\hat{\mu}^{2}_{f}$,

[TABLE]

In this limit a fourth order Taylor expansion thus provides the exact results for the basic bulk thermodynamic observables. This also is correct in leading order perturbation theory, i.e. at ${\cal O}(g^{2})$ Vuorinen:2003fs .

III Outline of the Calculation

III.1 Taylor series in baryon number, electric charge and strangeness

chemical potentials

Our goal is the calculation of Taylor expansion coefficients for basic bulk thermodynamic observables of strong-interaction matter in terms of chemical potentials $\mu_{X}$ for conserved charges ($X=B,\ Q,\ S)$. We start with the expansion of the pressure, $P$, in terms of the dimensionless ratios $\hat{\mu}_{X}\equiv\mu_{X}/T$, which are the logarithms of fugacities,

[TABLE]

with $\chi_{000}^{BQS}\equiv P(T,0)/T^{4}$. The chemical potentials for conserved charges are related to the quark chemical potentials $(\mu_{u},\mu_{d},\mu_{s})$,

[TABLE]

The expansion coefficients $\chi_{ijk}^{BQS}$, i.e. the so-called generalized susceptibilities, can be calculated at vanishing chemical potential222We often suppress the argument ($T$) of the generalized susceptibilities. We also suppress superscripts and subscripts of $\chi_{ijk}^{BQS}$ whenever one of the subscripts vanishes, e.g. $\chi_{i0k}^{BQS}\equiv\chi_{ik}^{BS}$.,

[TABLE]

From Eq. 6 it is straightforward to obtain the Taylor series for the number densities,

[TABLE]

This only requires knowledge of the expansion coefficients entering the series for $P/T^{4}$. The energy ($\epsilon$) and entropy ($s$) densities, on the other hand, also require derivatives of the generalized susceptibilities with respect to temperature, which are the expansion coefficients of the trace anomaly,

[TABLE]

with $i+j+k$ even and

[TABLE]

With this one finds for the Taylor expansions of the energy and entropy densities,

[TABLE]

III.2 Constrained series expansions

In our calculations we generated all generalized susceptibilities up to $6^{th}$ order, which are needed to set up the general Taylor series in terms of the three conserved charge chemical potentials as discussed in the previous subsection. In the following we will, however, consider only thermodynamic systems, in which the electric charge and strangeness chemical potentials are fixed by additional constraints and become functions of the baryon chemical potential and temperature. We only consider constraints that can be fulfilled order by order in the Taylor series expansion. That is, for the construction of the $6^{th}$ order Taylor series of the pressure in terms of $\hat{\mu}_{B}$ we need to know the expansion of $\hat{\mu}_{Q}(T,\mu_{B})$ and $\hat{\mu}_{S}(T,\mu_{B})$ up to fifth order in $\hat{\mu}_{B}$,

[TABLE]

The above parametrization includes the cases of vanishing electric charge and strangeness chemical potentials, $\mu_{Q}=\mu_{S}=0$, which we are going to discuss in the next section as well as the strangeness neutral case with fixed electric charge to baryon-number ratio, which we will analyze in Section V.

Implementing the constraints specified in Eq. 14 in the Taylor series for the pressure and net conserved-charge number densities one obtains series in terms of the baryon chemical potential only,

[TABLE]

Using

[TABLE]

and the series expansions of $\hat{\mu}_{Q}$ and $\hat{\mu}_{S}$ given in Eq. 14 one easily finds the relation between the expansion coefficients for the pressure and number densities,

[TABLE]

When imposing constraints on the electric charge and strangeness chemical potentials, these generally become temperature dependent functions as indicated in Eq. 14. The temperature derivative of $P/T^{4}$ at fixed $\hat{\mu}_{B}$ in the constraint case and the partial derivative of $P/T^{4}$ at fixed $(\hat{\mu}_{B},\hat{\mu}_{Q},\hat{\mu}_{S})$, which defines the trace anomaly $\Delta$ (Eq. 10), thus are related through

[TABLE]

where the (total) temperature derivative ${\rm d}/{\rm d}T$ is taken at fixed $\hat{\mu}_{B}$ and $\hat{\mu}_{X}^{\prime}={\rm d}\hat{\mu}_{X}/{\rm d}T$. With this we obtain the Taylor series for the trace anomaly,

[TABLE]

with

[TABLE]

We also introduce

[TABLE]

With this the Taylor series expansion of the energy and entropy densities for constraint cases, in which $\hat{\mu}_{Q}$ and $\hat{\mu}_{S}$ satisfy Eq. 14, becomes

[TABLE]

with $\epsilon_{2n}(T)=3P_{2n}(T)+TP_{2n}^{\prime}(T)-h_{2n}(T)$ and $\sigma_{2k}(t)=4P_{2n}(T)+TP_{2n}^{\prime}(T)-N_{2n-1}^{B}(T)-h_{2n}(T)-t_{2n}(T)$.

III.3 Numerical calculation of generalized susceptibilities up to

${\cal O}(\mu^{6})$

The generalized susceptibilities $\chi_{ijk}^{BQS}$ have been calculated on gauge field configurations generated for (2+1)-flavor QCD using the Highly Improved Staggered Quark (HISQ) action Follana:2006rc and the tree-level improved Symanzik gauge action.

All calculations are performed using a strange quark mass $m_{s}$ tuned to its physical value. We performed calculations with two different light to strange quark mass ratios, $m_{l}/m_{s}=1/27$ and $1/20$. The former corresponds to a pseudo-scalar Goldstone mass, which in the continuum limit yields a pion mass $m_{\pi}\simeq 140$ MeV, the latter leads to a pion mass $m_{\pi}\simeq 160$ MeV. These parameters are fixed using the line of constant physics determined by HotQCD from the $f_{K}$ scale. Using $f_{K}=155.7(9)/\sqrt{2}$ MeV allows to determine the lattice spacing $a(\beta)$ at a given value of the gauge coupling $\beta$ and the corresponding set of quark masses ($m_{l},m_{s})$, which in turn fixes the temperature on a lattice with temporal extent $N_{\tau}$, i.e. $T=(N_{\tau}a)^{-1}$. More details on the scale determination are given in Bazavov:2011nk .

All calculations have been performed on lattices of size $N_{\sigma}^{3}N_{\tau}$ with an aspect ratio $N_{\sigma}/N_{\tau}=4$. We perform calculations in the temperature interval $T\in[135\leavevmode\nobreak\ {\rm MeV},330\leavevmode\nobreak\ {\rm MeV}]$ using lattices with temporal extent $N_{\tau}=6,\ 8,\ 12$ and $16$, which corresponds to four different values of the lattice spacings at fixed temperature. At temperatures $T\leq 175$ MeV all calculations have been performed with the lighter, physical quark mass ratio $m_{l}/m_{s}=1/27$. In the high temperature region quark mass effects are small and we based our calculations on existing data sets for $m_{l}/m_{s}=1/20$, which have previously been generated by the HotQCD collaboration and used for the calculation of second order susceptibilities Bazavov:2012jq . These data sets have been extended for the calculation of higher order susceptibilities. Gauge field configurations are stored after every $10^{th}$ molecular dynamics trajectory of unit length.

All calculations of $4^{th}$ and $6^{th}$ order expansion coefficients have been performed on lattices with temporal extent $N_{\tau}=6$ and 8. In these cases we gathered a large amount of statistics. At low temperatures we have generated up to $1.2$ million trajectories for $N_{\tau}=6$ and up to $1.8$ million trajectories for $N_{\tau}=8$. At high temperature less than a tenth of this statistics turned out to be sufficient. The $2^{nd}$ order expansion coefficients have been calculated on lattices with four different temporal extends, $N_{\tau}=6,\ 8,\ 12,\ 16$. At fixed temperature this corresponds to four different values of the lattice cut-off, which we used to extract continuum extrapolated results for the second order expansion coefficients. We also extrapolated results for the higher order expansion coefficients to the continuum limit. However, having at hand results from only two lattice spacings for these expansion coefficients we consider these extrapolations as estimates of the results in the continuum limit.

On each configuration the traces of all operators needed to construct up to sixth order Taylor expansion coefficients have been calculated stochastically. For the calculation of $2^{nd}$ and $4^{th}$ order expansion coefficients we follow the standard approach of introducing a non-zero chemical potential in the QCD Lagrangian as an exponential prefactor for time-like gauge field variables Hasenfratz:1983ba , i.e. the chemical potential $\mu_{f}$ for quark flavor $f$ is introduced through a factor ${\rm e}^{\mu_{f}a}$ (${\rm e}^{-\mu_{f}a}$) on time-like links directed in the forward (backward) direction. This insures that all observables calculated are free of ultra-violet divergences. For the calculation of all $6^{th}$ order expansion coefficients we use the so-called linear-$\mu$ approach Gavai:2011uk ; Gavai:2014lia . This becomes possible as no ultra-violet divergences appear in $6^{th}$ order cumulants and above. In the linear-$\mu$ formulation the number of operators that contribute to cumulants is drastically reduced and their structure is simplified. All operators appearing in the exponential formulation, that involve second or higher order derivatives of the fermion matrix Allton:2005gk , vanish. The remaining operators are identical in both formulations. One thus only has to calculate traces of observables that are of the form,

[TABLE]

where $M_{f}$ is the staggered fermion matrix for light ($f=l$) or strange ($f=s$) quarks, respectively, and $M_{f}^{\prime}$ denotes its derivative with respect to the flavor chemical potential $\hat{\mu}_{f}$. The final error on these traces depends on the noise due to the use of stochastic estimators for the inversion of the fermion matrices $M_{f}$, as well as on the gauge noise resulting from a finite set of gauge configurations that get analyzed. We analyzed the signal to noise ratio for all traces of operators that we calculate and identified the operator $D_{1}=M_{f}^{-1}M^{\prime}_{f}$, as being particularly sensitive to the stochastic noise contribution. This operator has been measured using 2000 random noise vectors. For the calculation of traces of all other operators we used 500 random noise vectors. We checked that this suffices to reduce the stochastic noise well below the gauge noise. The simulation parameters and the statistics accumulated in this calculation are summarized in the tables of Appendix A.

All fits and continuum extrapolations shown in the following are based on spline interpolations with coefficients that are allowed to depend quadratically on the inverse temporal lattice size. Our fitting ansatz and the strategy followed to arrive at continuum extrapolated results are described in detail in Ref. Bazavov:2014pvz . For the current analysis we found it sufficient to use spline interpolations with quartic polynomials and 3 knots whose location is allowed to vary in the fit range.

IV Equation of state for $\mu_{Q}=\mu_{S}=0$

Let us first discuss the Taylor expansion for bulk thermodynamic observables in the case of vanishing electric charge and strangeness chemical potentials. This greatly simplifies the discussion and yet incorporates all the features of the more general case. Also the discussion of truncation errors presented in this section carries over to the more general situation.

IV.1 Pressure and net baryon-number density

For $\mu_{Q}=\mu_{S}=0$ the Taylor expansion coefficients $P_{2n}$ and $N_{2n-1}^{B}$, introduced in Eqs. 15 and 16, are simply related by

[TABLE]

The series for the pressure and net baryon-number density simplify to,

[TABLE]

In Eqs. 26 and 27 we have factored out the leading order (LO) $\mu_{B}$-dependent part in the series for the pressure as well as the net baryon-number density. This helps to develop a feeling for the importance of higher order contributions and, in particular, the approach to the HRG limit at low temperatures. Note that all ratios $\chi_{2n}^{B}/\chi_{2}^{B}$ are unity in a HRG and, in the infinite temperature, ideal quark gas limit, $\chi_{4}^{B}/\chi_{2}^{B}=2/(3\pi^{2})\simeq 0.068$ is the only non-vanishing higher order expansion coefficient. From Eqs. 26 and 27 it is evident that contributions from higher order expansion coefficients become more important in the number density than in the pressure. Relative to the LO result, the contributions of the NLO and NNLO expansion coefficients for $n_{B}/T^{3}$ are a factor two and three larger respectively than for the corresponding expansion coefficients in the pressure series.

We show the leading order coefficient $\chi_{2}^{B}(T)$ in Fig. 2 and the NLO ($\chi_{4}^{B}$) and NNLO ($\chi_{6}^{B}$) coefficients divided by $\chi_{2}^{B}(T)$ in Fig. 3. The left hand part of Fig. 2 shows the leading order contribution $\chi_{2}^{B}$ in the entire temperature interval used in the current analysis. For the LO expansion coefficients we also used data from simulations on $48^{3}\times 12$ lattices. Here we used existing data for $m_{l}/m_{s}=1/20$ Bazavov:2014pvz and generated new ensembles for $m_{l}/m_{s}=1/27$ at nine temperature values below $T=175$ MeV. Furthermore, we used data on $64^{3}\times 16$ lattices at a corresponding set of low temperature values. These data are taken from an ongoing calculation of higher order susceptibilities performed by the HotQCD Collaboration333We thank the HotQCD Collaboration for providing access to the second order quark number susceptibilities.. This allowed us to update the continuum extrapolation for $\chi_{2}^{B}$ given in Bazavov:2012jq . The new continuum extrapolation shown in Fig. 2 is consistent with our earlier results, but has significantly smaller errors in the low temperature region. In the right hand part of this figure we compare the continuum extrapolated lattice QCD data for $\chi_{2}^{B}$ with HRG model calculations. It is obvious that the continuum extrapolated QCD results overshoot results obtained from a conventional, non-interacting HRG model calculations with resonances taken from the particle data tables (PDG-HRG) and treated as point-like excitations. We therefore compare the QCD results also with a HRG model that includes additional strange baryons,which are not listed in the PDG but are predicted in quark models and lattice QCD calculations. We successfully used such an extended HRG model (QM-HRG) in previous calculations Majumder:2010ik ; Bazavov:2014xya . As can be seen in Fig. 2 (left), continuum extrapolated results for $\chi_{2}^{B}$ agree well with QM-HRG calculations.

As can be seen in the left hand part of Fig. 3, the ratio $\chi_{4}^{B}/\chi_{2}^{B}$ approaches unity with decreasing temperature, but is small at high temperatures where the leading order correction is large. The relative contribution of the NLO correction thus is largest in the hadronic phase, where $\chi_{4}^{B}/\chi_{2}^{B}\simeq 1$. For temperatures T\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}155 MeV we find $\chi_{4}^{B}/\chi_{2}^{B}\leq 0.8$. The relative contribution of the NLO correction to the $\mu_{B}$-dependent part of the pressure (number density) in the crossover region and below thus is about 8% (16%) at $\mu_{B}/T=1$ and rises to about 33% (66%) at $\mu_{B}/T=2$. At temperatures larger than $180$ MeV the relative contribution of the NLO correction to pressure and number density at $\mu_{B}/T=2$ is less than 8% and 16%, respectively.

The relative contribution of the ${\cal O}(\hat{\mu}_{B}^{6})$ correction, $\chi_{6}^{B}/\chi_{2}^{B}$, is shown in the right hand part of Fig. 3. The ideal gas limit for this ratio vanishes. Obviously the ratio is already small for all temperatures $T>180$ MeV, i.e. $\chi_{6}^{B}/\chi_{2}^{B}\leq 0.5$. Consequently, for $\hat{\mu}_{B}=2$ the correction to the leading order result is less than 2.2% for the $\mu_{B}$-dependent part of the pressure and less than 7% for the net baryon-number density. At lower temperatures the statistical errors on current results for $\chi_{6}^{B}/\chi_{2}^{B}$ are still large. However, a crude estimate for the magnitude of this ratio at all temperatures larger than $130$ MeV suggests, $\left|\chi_{6}^{B}/\chi_{2}^{B}\right|\leq 3$. In the low temperature, hadronic regime and for $\hat{\mu}_{B}=2$ the ${\cal O}(\hat{\mu}_{B}^{6})$ corrections to the $\mu_{B}$-dependent part of the pressure can be about 13%. However, in the total pressure, which also receives large contributions from the meson sector, this will result only in an error of less than 3%. In the calculation of the net baryon-number density, on the other hand, the current uncertainty on ${\cal O}(\hat{\mu}_{B}^{6})$ expansion coefficients results in errors of about 40% at temperatures below $T\simeq 155$ MeV. In fact, as discussed already in section II, higher order corrections are larger in the Taylor expansion of the number density. From Eq. 25 it follows for the ratio of NLO and LO expansion coefficients, $N_{5}^{B}/N_{1}^{B}=3P_{6}/P_{2}$. Clearly better statistics is needed in the low temperature range to control higher order corrections to $n_{B}/T^{3}$.

In Fig. 4 we show results for the $\mu_{B}$-dependent part of the pressure (left) and the net baryon-number density (right) calculated from Taylor series up to and including LO, NLO and NNLO contributions, respectively. This suggests that up to $\mu_{B}\simeq 2T$ results for the pressure at low temperature are well described by a Taylor series truncated at NNLO, while at higher temperature NNLO corrections are small even at $\mu_{B}\simeq 3T$. This also is the case for $n_{B}/T^{3}$, although the NNLO correction is large at low temperatures and, at present, does not allow for a detailed quantitative analysis of the baryon-number density in this temperature range.

It also is obvious that the Taylor series for the pressure and $n_{B}/T^{3}$ in the temperature range up to $T\simeq 180$ MeV are sensitive to the negative contributions of the $6^{th}$ order expansion coefficient. The occurrence of a dip in the sixth order expansion coefficient of the pressure has been expected to show up on the basis of general scaling arguments for higher order derivatives of the QCD pressure in the vicinity of the chiral phase transition Friman:2011pf . It may, however, also reflect the influence of a singularity on the imaginary chemical potential axis Bonati:2016pwz (Roberge-Weiss critical point Roberge ) on Taylor series of bulk thermodynamic observables in QCD. Even with improved statistics it thus is expected that the wiggles, that start to show up in the expansion of pressure and net baryon-number density above $\mu_{B}/T\simeq 2$ (see Fig. 4) and reflect the change of sign in the sixth order expansion coefficient, will persist. Getting the magnitude of the dip in $\chi_{6}^{B}/\chi_{2}^{B}$ at $T\simeq 160$ MeV under control in future calculations thus is of importance for the understanding of this non-perturbative regime of the QCD equation of state in the high temperature phase close to the transition region. This also indicates that higher order corrections need to be calculated in order to control the equation of state in this temperature regime.

IV.2 Net strangeness and net electric charge densities

For vanishing strangeness and electric charge chemical potentials the corresponding net strangeness ($n_{S}$) and net electric charge ($n_{Q}$) densities are nonetheless non-zero because the carriers of these quantum numbers also carry baryon number. The ratios of number densities are given by

[TABLE]

In a hadron resonance gas the ratios $n_{S}/n_{B}$ and $n_{Q}/n_{B}$ are independent of the baryon chemical potential and, irrespective of the value of $\hat{\mu}_{B}$, these ratios approach $-1$ and [math], respectively, in the $T\rightarrow\infty$ limit. One thus may expect that these ratios only show a mild dependence on $\hat{\mu}_{B}$, which indeed is apparent from the results of the NNLO expansions shown in Fig. 5.

For $\mu_{Q}=\mu_{S}=0$ non-vanishing electric charge and strangeness densities only arise due to a non-zero baryon-chemical potential. In the low temperature HRG phase $n_{Q}$ and $n_{S}$ thus only receive contributions from charged baryons or strange baryons, respectively. The ratios $n_{Q}/n_{B}$ and $n_{S}/n_{B}$ thus are sensitive to the particle content in a hadron resonance gas and a comparison with PDG-HRG and QM-HRG is particularly sensitive to the differences in the baryon content in these two models. It is apparent from Fig. 5 that at low temperatures the QM-HRG model provides a better description of the lattice QCD results than the PDG-HRG model.

IV.3 The energy and entropy densities

In order to calculate the energy and entropy densities, defined in Eqs. 23 and 24, we need to extract the temperature derivative of the expansion coefficients of the pressure. We use as a starting point the representation of the pressure given in Eq. 26 and calculate the temperature derivatives of $\chi_{n}^{B}$ from the splines used to fit this observable. With this we construct the expansion coefficients $\epsilon_{n}^{B}(T)$ and $\sigma_{n}^{B}$ defined in Eqs. 12 and 13,

[TABLE]

We show the LO and NLO expansion coefficients for energy and entropy densities together with the expansion coefficient for the pressure in Fig. 6. Because of Eq. 25 the expansion coefficients of the net baryon-number density are simply proportional to those of the pressure.

Clearly the temperature dependence of the expansion coefficients of the energy and entropy densities shows more structure than in the case of the pressure. Qualitatively this can be understood in terms of the pseudo-critical behavior of bulk thermodynamic observables. Once thermodynamic quantities are dominated by contributions from the singular part of the free energy, which is expected to happen in the transition region, they become functions of $(T-T_{c})+\kappa\hat{\mu}_{B}^{2}$. The temperature derivative of the expansion coefficient $P_{2}$, which gives $\epsilon_{2}$, thus will show properties similar to those of $P_{4}$. The LO correction $\epsilon_{2}^{B}/T^{4}$ has a mild peak, which results from the strongly peaked $T$-derivative of $\chi_{2}^{B}$ which is qualitatively similar to $\chi_{4}^{B}$, and the NLO correction is negative in a small temperature interval above $T_{c}$, which arises from the negative $T$-derivative of $\chi_{4}^{B}$ at high temperature, which resembles the negative part of $\chi_{6}^{B}$ at high temperature.

Although the temperature dependence of $\epsilon_{n}$ and $\sigma_{n}$ differs from that of the pressure coefficient, $P_{n}$, the conclusions drawn for the relative strength of the expansion coefficients are identical in all cases. As can be seen from the inset in Fig. 6 (right) the relative contribution of the NLO expansion coefficients never exceeds 10%. In particular, at temperatures larger than $180$ MeV the magnitude of the NLO expansion coefficients never exceeds 2% of the LO expansion coefficients. Again this leads to the conclusion that at $\mu_{B}/T=2$ and temperatures above $180$ MeV the NLO correction contributes less than 8% of the leading correction to $\mu_{B}$-dependent part of the energy and entropy densities. For T\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}155 MeV, however, the NLO contribution can rise to about 30%. A similar conclusion holds for the ${\cal O}(\hat{\mu}_{B}^{6})$ corrections, although it requires higher statistics to better quantify the magnitude of this contribution. In Fig. 7 we show results for the total pressure and total energy density. For $P/T^{4}$ and $\epsilon/T^{4}$ at $\mu_{B}=0$ we used the results obtained by the HotQCD Collaboration Bazavov:2014pvz and added to it the results from the ${\cal O}(\hat{\mu}_{B}^{6})$ expansions presented above. This figure also makes it clear that despite of the large error of higher order expansion coefficients, which we have discused above, the error on the total pressure and energy density still is dominated by errors on their values at $\mu_{B}=0$.

V Equation of state in strangeness neutral systems

V.1 Taylor expansion of pressure, baryon-number, energy

and entropy densities

We now want to discuss the equation of state for strangeness neutral systems with a fixed ratio of electric charge to baryon-number density, i.e. we impose the constraints Bazavov:2012vg

[TABLE]

These constraints can be realized through suitable choices of the electric charge and strangeness chemical potentials. This thus is a particular case of the constraint expansion discussed in Subsection III.2. The expansion coefficients $q_{n}$, $s_{n}$, $n=1,\ 3,\ 5$ needed to satisfy these constraints are given in Appendix B. For $r=0.4$ the constrained EoS obtained in this way is usually considered to be most appropriate for applications to heavy ion collisions. We will, however, in the following also comment on other choices of $r$, including the case of isospin symmetric systems ($r=1/2$) and electric charge neutral matter ($r=0$).

Using the constraints specified in Eq. 31 and the definition of the pressure in terms of generalized susceptibilities, $\chi_{ijk}^{BQS}$, the expansion coefficients $P_{2n}$ can easily be determined. Here it advantageous to use the relation between the Taylor expansion coefficients of the pressure, $P_{2n}$, and number densities, $N_{2n-1}^{X}$, given in Eq. 18, which simplifies considerably for strangeness neutral systems. It now involves only the net baryon-number density coefficients,

[TABLE]

Explicit expressions for all $N_{n-1}^{B}$ and $q_{n-1}$, for $n=2,\ 4,\ 6$, are given in Appendix B. The resulting expansion coefficients for the pressure are shown in Fig. 8. Also shown in the bottom-right panel of this figure is the ratio of the expansion coefficients for the net baryon-number density, $N_{n-1}^{B}$ and the appropriately rescaled expansion coefficients of the pressure, $nP_{n}$. In electric charge neutral systems, $r=0$ as well as in the isospin symmetric limit $r=1/2$, for which the expansion coefficients $q_{i}=0$ vanish for all $i$, this ratio is unity. In both cases the simple relation given in Eq. 25 holds. Also for other values of $r$ the contribution from terms proportional to $r$ are small. In Fig. 8 (bottom, right) we show the ratio $N_{2n-1}^{B}/nP_{n}$ for the case $r=0.4$ and $n=2$ and $4$, respectively. At ${\cal O}(\hat{\mu}_{B}^{2})$ differences between $N_{1}^{B}$ and $2P_{2}$ never exceed 2% and at ${\cal O}(\hat{\mu}_{B}^{4})$ the difference between $N_{3}^{B}$ and $4P_{4}$ varies between 3% at low temperature and -6% at high temperature. In the infinite temperature ideal gas limit the ratios become $N^{B}_{1}/2P_{2}=1.018$ and $N^{B}_{3}/4P_{4}=0.927$, respectively.

In general one finds that the dependence of bulk thermodynamic observables on the net electric charge to net baryon number-ratio is weak. The ${\cal O}(\hat{\mu}_{B}^{2})$ expansion coefficient of the pressure in strangeness neutral systems differs by at most 10% in electric charge neutral ($r=0$) and isospin symmetric systems ($r=1/2$), respectively. The expansion coefficient $P_{2}$ evaluated for different values of $r$ is shown in Fig. 9. For chemical potentials $\hat{\mu}\leq 2$ this amounts to differences less than 1.5% of the total pressure. On the other hand, strangeness neutral systems differ substantially from systems with vanishing strangeness chemical potential. In this case the ${\cal O}(\hat{\mu}_{B}^{2})$ expansion coefficients differ by almost 50% in the high temperature limit. For $T<150$ MeV this difference is only about 10% reflecting that the different treatment of the strangeness sector becomes less important for the thermodynamics at low temperature. This is also shown in Fig. 9.

Compared to the leading ${\cal O}(\hat{\mu}_{B}^{2})$ contributions to bulk thermodynamic observables the ${\cal O}(\hat{\mu}_{B}^{4})$ and ${\cal O}(\hat{\mu}_{B}^{6})$ corrections are smaller in the strangeness neutral case than in the case $\mu_{Q}=\mu_{S}=0$, which we have discussed in the previous section. This is evident from Fig. 10, where we show the ratios $12P_{4}/P_{2}$ and $360P_{6}/P_{2}$. These combinations are unity in a HRG with $\mu_{S}=\mu_{Q}=0$ but smaller than unity in the strangeness neutral case. Higher order corrections in Taylor series for strangeness neutral systems thus are of less importance than in the case $\mu_{S}=0$. This also means that the errors, which are large on e.g. sixth order expansion coefficients, are of less importance for the overall error budget of Taylor expansions in strangeness neutral systems. This is indeed reflected in the $\mu_{B}$-dependence of $(P(T,\mu_{B})-P(T,0))/T^{4}$ and $n_{B}(T,\mu_{B})/T^{3}$ shown in the upper panels of Fig. 11 for the case $r=0.4$. As can be seen in these two figures, at low temperatures the $\mu_{B}$-dependent part of the pressure as well as the net baryon-number density agree quite well with HRG model calculations that describe the thermodynamics of a gas of non-interacting, point-like hadron resonances. This agreement, however, gets worse at larger values of $\mu_{B}$. Not unexpectedly, at higher temperatures deviations from HRG model calculations become large already at small values of $\mu_{B}$. This is apparent from the lower two panels of Fig. 11, where we show the ratio of the $\mu_{B}$-dependent part of the pressure and the corresponding HRG model result (left) and the net baryon-number density divided by the corresponding HRG model result (right). In the HRG model calculation $(P(T,\mu_{B})-P(T,0))/T^{4}$ as well as $n_{B}(T,\mu_{B})/T^{3}$ only depend on the baryon sector of the hadron spectrum. The results shown in Fig. 11 thus strongly suggest that HRG model calculations using resonance spectra in model calculations for non-interacting, point-like hadron gases may be appropriate (within $\sim 10$% accuracy) to describe the physics in the crossover region of strongly interacting matter at vanishing or small values of the baryon chemical potential, but fail444It has been pointed out that the point-like particle approximation is appropriate in the meson sector but not in the baryon sector at high density. Introducing a non-zero size of hadron resonances Hagedorn:1980kb ; Dixit:1980zj may, for some observables, improve the comparison with QCD thermodynamics Andronic:2012ut ; Vovchenko:2016rkn . However, it seems that the introduction of several additional parameters will be needed to achieve overall good agreement with the many observables calculated now in QCD in the temperature range of interest, i.e. in the crossover region from a hadron gas to strongly interacting quark-gluon matter. to do so at large $\mu_{B}/T$ and/or T\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}160 MeV. At $T=165$ MeV QCD and HRG model results for the net baryon-number density differ by 40% at $\mu_{B}/T=2$. This has consequences for the determination of freeze-out conditions in heavy ion collisions. We will come back to this discussion in Section VI.

The $\mu_{B}$-dependent contributions to the energy and entropy densities have been defined in Eqs. 23 and 24. In strangeness neutral systems the expansion coefficients simplify considerably,

[TABLE]

Results for the ${\cal O}(\mu_{B}^{2})$ and ${\cal O}(\mu_{B}^{4})$ expansion coefficients are shown in Fig. 12 together with the corresponding expansion coefficients for the pressure and net baryon-number density. Results for the total energy density as well as the total pressure for $\mu_{B}/T=0$ and $2$ are shown in Fig. 13. As discussed in the previous section also here it is evident that current errors on the total pressure and energy density are dominated by errors on these observables at $\mu_{B}=0$.

In Fig. 13 we also show results for the total pressure obtained within the stout discretization scheme. The result for $\hat{\mu}_{B}=0$ is taken from Borsanyi:2013bia . The $\hat{\mu}_{B}$-dependent contribution is based on calculations with an imaginary chemical potential Gunther:2016vcp . These results have been analytically continued to real values of $\hat{\mu}_{B}$ using a $6^{th}$ order polynomial in $\hat{\mu}_{B}$. As can be seen the total pressure agrees quite well with the results obtained with a sixth order Taylor expansion, although the results obtained the analytic continuation within the stout discretization scheme tend to stay systematically below the central values obtained from the analysis of Taylor series expansions in the HISQ discretization scheme.

V.2 Parametrization of the equation of state

At $\mu_{B}=0$ the HotQCD Collaboration presented a parametrization of the pressure, obtained as interpolating curves for the continuum extrapolated fit, that also provided an adequate description of all the other basic thermodynamic quantities, i.e. the energy and entropy densities as well as the specific heat and the velocity of sound Bazavov:2014pvz . Here we want to extend this parametrization to the case $\hat{\mu}_{B}>0$. Similar to what has been done at $\mu_{B}=0$ it turns out that a ratio of fourth order polynomials in the inverse temperature is flexible enough to describe the temperature dependence of all required Taylor expansion coefficients in the temperature range $T\in[130\leavevmode\nobreak\ {\rm MeV},280\leavevmode\nobreak\ {\rm MeV}]$. We use such an ansatz for the three expansion coefficients of the net baryon-number density ($N_{1}^{B},\ N_{3}^{B},\ N_{5}^{B}$) and the three electric charge chemical potentials ($q_{1},\ q_{3},\ q_{5}$). This suffices to calculate all thermodynamic observables in strangeness neutral systems.

We use a ratio of fourth order polynomials in $1/T$ as an ansatz for the expansion coefficients of the net baryon-number density,

[TABLE]

Here $\bar{t}=T_{c}/T$ and the QCD transition temperature $T_{c}=154$ MeV is used as a convenient normalization. Similarly we define the parametrization of the expansion coefficients for the electric charge chemical potential,

[TABLE]

The parameters for these interpolating curves are summarized in Table 1.

The expansion coefficients of the pressure are then obtained by using Eqs. 32-34. The resulting interpolating curves for $P_{k}$ are shown as darker curves in Fig. 8. All other interpolating curves shown as darker curves in other figures have been obtained by using the above interpolations. In particular, interpolating curves for the energy and entropy densities are obtained by using Eqs. 35 and 36 and calculating analytically temperature derivatives of the parametrizations of $P_{n}$ and $q_{n}$ given in Eqs. 37 and 38. The resulting interpolating curves for the second and fourth order Taylor expansion coefficients are shown in Fig. 12.

We also used a ratio of fourth order polynomials to interpolate results for the pressure at $\mu_{B}=0$. We write the pressure as

[TABLE]

The coefficients $p_{in}$ and $p_{id}$ are also given in Table 1.

VI Lines of constant physics to ${\cal O}(\mu_{B}^{4})$

We want to use here the Taylor series for bulk thermodynamic observables, i.e. the pressure, energy and entropy densities, to discuss contour lines in the $T$-$\mu_{B}$ plane on which these observables stay constant. It has been argued quite successfully that the thermal conditions at the time of chemical freeze-out in heavy ion collisions can be characterized by lines in the $T$-$\mu_{B}$ plane on which certain thermodynamic observables or ratios thereof stay constant Cleymans:1999st ; Cleymans:2005xv , although the freeze-out mechanism in the rapidly expanding fireball created in a heavy ion collision is of dynamical origin and will in detail be more complicated (see for instance Tomasik:2002qt ). While lines of constant physics (LCPs) involving total baryon-number densities, as used in Cleymans:1999st ; Cleymans:2005xv , are not appropriate for calculations within the framework of quantum field theories, other criteria like lines of constant entropy density in units of $T^{3}$ Cleymans:2004hj or constant pressure Rafelski:2009jr ; Petran:2013qla ; Rafelski:2015hta have been suggested to characterize freeze-out parameters $(T_{f},\mu_{B}^{f})$ corresponding to heavy ion collisions at different values of the beam energy ($\sqrt{s_{NN}}$). Generally such criteria have been established by comparing experimental data with model calculations based on some version of a HRG model. We will determine here LCPs from the lattice QCD calculations of pressure, energy and entropy densities and confront them with freeze-out parameters that have been obtained by comparing particle yields, measured at different values of $\sqrt{s_{NN}}$, to HRG model calculations.

We consider an observable $f(T,\mu_{B})$, i.e. the pressure, energy density or entropy density which are even functions of $\mu_{B}$. We parametrize a ’line of constant $f$’ by,

[TABLE]

In order to determine the expansion coefficients $\kappa_{2}^{f}$ and $\kappa_{4}^{f}$ we need to expand the function $f(T,\mu_{B})$ up to $4^{th}$ order in $\mu_{B}$ and up to second order in $T$ around some point $(T_{0},0)$,

[TABLE]

Note that we expand here in terms of $\mu_{B}$ rather than in $\hat{\mu}_{B}\equiv\mu_{B}/T$. Replacing the temperature $T$ in Eq. VI by the ansatz for a line of constant $f$, Eq. 40, and keeping terms up to ${\cal O}(\mu_{B}^{4})$ gives

[TABLE]

We then can determine $\kappa_{2}^{f}$ and $\kappa_{4}^{f}$ by demanding that the expansion coefficients at ${\cal O}(\mu_{B}^{2})$ and ${\cal O}(\mu_{B}^{4})$ vanish, i.e.

[TABLE]

As we will deal with observables that are given as a Taylor series in $\hat{\mu}_{B}$ at fixed $T$, i.e. $f(T,\mu_{B})=\sum_{k=0}^{\infty}f_{2k}(T)\hat{\mu}_{B}^{2k}$, the derivatives with respect to $\mu_{B}$ appearing in Eqs. 42 and 43 can be replaced by suitable Taylor expansion coefficients of $f(T,\mu_{B})$,

[TABLE]

We will in the following work out detailed expressions for the quadratic correction coefficient, $\kappa_{2}^{f}$, for lines of constant pressure ($f\equiv P$), energy density ($f\equiv\epsilon$) and entropy density ($f\equiv s$) in strangeness neutral systems with electric charge to net baryon-number ratio $r=0.4$. Details for the quartic coefficient, $\kappa_{4}^{f}$, are given in Appendix C.

• pressure $f\equiv P$:

The function $f(T,\mu_{B})$ is given by $P=T^{4}\sum_{n}P_{n}(\mu_{B}/T)^{n}$, with $P_{0}=P(T,0)/T^{4}$ denoting the pressure in units of $T^{4}$ at vanishing baryon chemical potential and $P_{n}(T)$, $n>0$, denoting the expansion coefficients of $P(T,\mu_{B})/T^{4}$ as introduced in Eq. 15. In the denominator of Eq. 44 we use the thermodynamic relation between pressure and entropy density $s=(\partial P/\partial T)_{\mu_{B}}$. The numerator is given by $f_{2}(T)=T^{4}P_{2}(T)$. This gives

[TABLE]

where $s/T^{3}$ is evaluated at $\hat{\mu}_{B}=0$.

• energy density $f\equiv\epsilon$:

The function $f(T,\mu_{B})$ is given by $\epsilon=T^{4}\sum_{n}\epsilon_{n}(\mu_{B}/T)^{n}$, with $\epsilon_{0}=\epsilon(T,0)/T^{4}$ denoting the energy density in units of $T^{4}$ at vanishing baryon chemical potential . In the denominator of Eq. 44 we use the thermodynamic relation between energy density and specific heat $C_{V}=(\partial\epsilon/\partial T)_{\mu_{B}}$ In the numerator we have $f_{2}(T)=T^{4}\epsilon_{2}(T)$. This gives

[TABLE]

where $C_{V}/T^{3}$ is evaluated at $\hat{\mu}_{B}=0$.

• entropy density $f\equiv s$:

The function $f(T,\mu_{B})$ is given by $s=(\epsilon+P-\mu_{B}n_{B}-\mu_{Q}n_{Q})/T=(\epsilon+P-\mu_{B}n_{B}(1+r\mu_{Q}/\mu_{B}))/T$. As $n_{B}$ is of ${\cal O}(\mu_{B})$ we need for the ratio of electric charge and strangeness chemical potentials only the leading order relation $\mu_{Q}/\mu_{B}=q_{1}$ defined in Eq. 14. In the denominator we use,

[TABLE]

In the numerator we have $f_{2}(T)=T^{3}(\epsilon_{2}+P_{2}-N_{1}^{B}(1+rq_{1}))$. With this we get,

[TABLE]

where we have used Eq. 32 to replace $N_{1}^{B}$ in favor of $P_{2}$.

We note that $\kappa_{2}^{\epsilon}>\kappa_{2}^{s}$, i.e. with increasing $\mu_{B}$ the entropy density decreases on lines of constant energy density.

The second order coefficients for the lines of constant physics thus can directly be calculated using the continuum extrapolated results for the pressure and energy density obtained at vanishing chemical potential in Bazavov:2014pvz and the leading order expansion coefficient of the pressure shown in Fig. 10. Similarly we obtain the quartic coefficients from the fourth order expansion of the pressure using the relations given in Appendix C. We show results for $\kappa_{2}^{f}$ and $\kappa_{4}^{f}$ in Fig. 14.

In the interval around $T_{c}$, i.e. $T\in[145\leavevmode\nobreak\ {\rm MeV},165\leavevmode\nobreak\ {\rm MeV}]$ we find,

[TABLE]

Apparently, at ${\cal O}(\mu_{B}^{2})$, lines of constant pressure and constant energy or entropy densities agree quite well and they also agree, within currently large errors, with the curvature of the transition line in (2+1)-flavor QCD. The coefficient of the quartic correction for the contour lines turns out to be about two orders of magnitude smaller than the leading order coefficients. This, of course, reflects the small contribution of the NLO corrections to the $\mu_{B}$-dependent part of pressure and energy density. For all fourth order coefficients we find $|\kappa_{4}^{f}|\leq 0.00024$ in the temperature interval around $T_{c}$. For $\mu_{B}/T\leq 2$ the contribution arising from $\kappa_{4}^{f}$ only leads to modifications of $T_{f}(\mu_{B})$ that stays within the error band arising from the uncertainty in $\kappa_{2}^{f}$.

The resulting lines of constant physics in the $T$-$\mu_{B}$ plane are shown in Fig. 15 (left) for three values of the temperature, $T=145$ MeV, $155$ MeV and $165$ MeV. These correspond to constant energy densities $\epsilon=0.203(27)$ GeV/fm3, $0.346(41)$ GeV/fm3 and $0.556(57)$ GeV/fm3, which roughly correspond to the energy density of cold nuclear matter, a hard sphere gas of nucleons at dense packing and the interior of a nucleus, respectively. Values of other bulk thermodynamic observables characterizing these LCPs are summarized in Table 2. The corresponding net baryon-number densities on these LCPs are shown in Fig. 15 (right). It is apparent from Fig. 15 (left) that LCPs for constant pressure, energy or entropy density agree well with each other up to baryon chemical potentials $\mu_{B}/T=2$, where the difference in temperature on different LCPs is at most $2$ MeV. We also note that the temperature on a LCP varies by about $7$ MeV or, equivalently, 5% between $\hat{\mu}_{B}=0$ and $\hat{\mu}_{B}=2$. Thus on a line of constant pressure, the entropy in units of $T^{3}$ changes by about 15%. I.e. constant $P$ or constant $s/T^{3}$, which both have been suggested as phenomenological descriptions for freeze-out conditions in heavy ion collisions, can not hold simultaneously, although a change of 15% of one of these observables may phenomenologically not be of much relevance. We also stress that at large values of $\hat{\mu}_{B}$ the comparison of experimental data with HRG model calculations, e.g. the use of single particle Boltzmann distributions used to extract freeze-out temperatures and chemical potentials, becomes questionable. As shown in Fig. 11 net baryon-number densities extracted from HRG and QCD calculations differ substantially at $\mu_{B}/T\simeq 2$.

Also shown in Fig. 15 (left) are results on freeze-out parameters and hadronization temperatures extracted from particle yields measured in heavy ion experiments Das:2014qca ; Floris:2014pta ; Becattini:2016xct by comparing data with model calculations based on the hadron resonance gas models. The region $\mu_{B}/T\leq 2$ corresponds to beam energies $\sqrt{s_{NN}}\geq 11.4$ GeV in the RHIC beam energy scan. Obviously, the freeze-out parameters extracted from the beam energy scan data Das:2014qca do not follow any of the LCPs. However, the discrepancy between the freeze-out parameters determined at the LHC Floris:2014pta and the highest beam energy at RHIC Das:2014qca suggests that also these determinations are not consistent among each other.

Finally we note that the lines of constant physics discussed above compare also well with the crossover line for the QCD transition. At non-zero values of the baryon chemical potential the change of the (pseudo)-critical temperature has been determined, using various approaches at real Kaczmarek:2011zz ; Endrodi:2011gv and imaginary Bonati:2015bha ; Bellwied:2015rza ; Cea:2015cya values of the chemical potential. To leading order one obtains,

[TABLE]

with $\kappa^{c}_{2}$ ranging from $0.0066(7)$ Kaczmarek:2011zz ; Endrodi:2011gv to $0.0135(20)$ Bonati:2015bha , $0.0149(21)$ Bellwied:2015rza and $0.020(4)$ Cea:2015cya . Lines that cover this spread in curvature parameters are also shown in Fig. 15 (left) for $T_{c}(0)=155$ MeV. While a small curvature for the crossover line would suggest that the crossover transition happens under more or less identical bulk thermodynamic conditions a large curvature obviously would indicate that the crossover transition happens already at significantly smaller values of pressure and energy density as $\mu_{B}/T$ increases.

VII Radius of convergence and the critical point

As discussed in the previous sections we generally find that the Taylor series for all basic thermodynamic quantities converge well for values of baryon chemical potentials $\mu_{B}\leq 2T$. Even in the low temperature regime the relative contribution of higher order expansion coefficients are generally smaller than in corresponding HRG model calculations. This, of course, also has consequences for our current understanding of the location of a possible critical point in the QCD phase diagram.

The results on the expansion coefficients of the Taylor series for e.g. the pressure can be cast into estimates for the location of a possible critical point in the QCD phase diagram. In general the radius of convergence can be obtained from ratios of subsequent expansion coefficients in the Taylor series for the pressure. Equally well one may use one of the derivatives of the pressure series. As one has to rely on estimates of the radius of convergence that generally are based on a rather short series, it may indeed be of advantage to use as a starting point the series for the net baryon-number susceptibility Gavai:2004sd , which diverges at the critical point, but still contains information from all expansion coefficients of the pressure series. The radius of convergence of this series is identical to that of the pressure. Model calculations also suggest that the estimators obtained from the susceptibility series converge faster to the true radius of convergence Karsch:2011yq . For $\mu_{Q}=\mu_{S}=0$ the expansion coefficients of the Taylor series for the net baryon-number susceptibility are again simply related to that of the pressure,

[TABLE]

From this one obtains estimators for the radius of convergence of the pressure and susceptibility series,

[TABLE]

Both estimators converge to the true radius of convergence in the limit $n\rightarrow\infty$. In order for this to correspond to a singularity at real values of $\hat{\mu}_{B}$, all expansion coefficients should asymptotically stay positive.

Obviously, the estimators $r_{2n}^{P}$ and $r_{2n}^{\chi}$ are proportional to each other, $r_{2n}^{P}=\sqrt{(2n+2)(2n+1)/[2n(2n-1)]}r_{2n}^{\chi}$. The difference between these to estimators may be taken as a systematic error for any estimate of the radius of convergence obtained from a truncated Taylor series. In the hadron resonance gas limit one finds for estimators involving sixth order cumulants, $r_{4}^{P}=1.58r_{4}^{\chi}$. In the following we restrict our discussion to an analysis of $r_{2n}^{\chi}$, which at finite $n$ leads to the smaller estimator for the radius of convergence. This seems to be appropriate in the present situation where we only can construct two independent estimators from ratios of three distinct susceptibilities. We thus may hope to identify regions in the QCD phase diagram at small values of $\hat{\mu}_{B}$ which are unlikely locations for a possible critical point.

An immediate consequence of the definitions given in Eq. 53 is that the ratios of generalized susceptibilities need to grow asymptotically like $|\chi_{n+2}^{B}/\chi_{n}^{B}|\sim n^{2}$ in order to arrive in the limit $n\rightarrow\infty$ at a finite value for the radius of convergence. At least for large values of $n$ one thus needs to find large deviations from the hadron resonance gas results $|\chi_{n+2}^{B}/\chi_{n}^{B}|^{HRG}=1$. As is obvious from the results presented in the previous sections, in particular from Fig. 3, the analysis of up to sixth order Taylor expansion coefficients does not provide any hints for such large deviations. The ratio $\chi_{4}^{B}/\chi_{2}^{B}$ turns out to be less than unity in the entire temperature range explored so far, i.e. for $T\geq 135$ MeV or $T/T_{c}>0.87(6)$. Below the crossover temperature, $T\sim 155$ MeV, the sixth order expansion coefficients also are consistent with HRG model results. They still have large errors. However, using the upper value of the error for $\chi_{6}^{B}/\chi_{4}^{B}$ provides a lower limit for the value of the estimator $r_{4}^{\chi}$. For temperatures in the interval $135\leavevmode\nobreak\ {\rm MeV}\leq T\leq 155\leavevmode\nobreak\ {\rm MeV}$ (or equivalently $0.87(5)\leq T/T_{c}\leq 1$) we currently obtain a lower limit on $r_{4}^{\chi}$ from the estimate $\chi_{6}^{B}/\chi_{4}^{B}\simeq\chi_{6}^{B}/\chi_{2}^{B}<3$. This converts into the bound $r_{4}^{\chi}\geq 2$, which is consistent with our observation that the Taylor series of all thermodynamic observables discussed in the previous sections is well behaved up to $\mu_{B}=2T$. A more detailed analysis, using the current errors on $\chi_{6}^{B}/\chi_{4}^{B}$ at five temperature values below and in the crossover region of the transition at $\mu_{B}=0$, is shown in Fig. 16. This shows that the bound arising from $r_{4}^{\chi}$ is actually more stringent at temperatures closer to $T_{c}$, where $\chi_{6}$ starts to become small and eventually tends to become negative.

These findings are consistent with recent results for susceptibility ratios obtained from calculations with an imaginary chemical potential DElia:2016jqh . Also in that case all susceptibility ratios are consistent with HRG model results. At present one thus cannot rule out that the radius of convergence may actually be infinite. Results for $r_{4}^{\chi}$ obtained in Ref. DElia:2016jqh lead to even larger estimators for the radius of convergence than our current lower bound. This is also shown in Fig. 16.

The observations and conclusions discussed above are in contrast to estimates for the location of a critical point obtained from a calculation based on a reweighting technique Fodor:2004nz as well as from Taylor series expansion in 2-flavor QCD Datta:2014zqa ; Datta:2016ukp . Both these calculations have been performed with unimproved staggered fermion discretization schemes and thus may suffer from large cut-off effects. Moreover, the latter calculation also suffers from large statistical errors on higher order susceptibilities. Results from Ref. Fodor:2004nz and Ref. Datta:2016ukp are also shown in Fig. 16.

We thus conclude from our current analysis that a critical point at chemical potentials smaller than $\mu_{B}=2T$ is strongly disfavored in the temperature range $135\leavevmode\nobreak\ {\rm MeV}\leq T\leq 155$ MeV and its location at higher values of temperature seems to be ruled out. Our results suggest that the radius of convergence in that temperature interval will turn out to be significantly larger than the current bound once the statistics on $6^{th}$ order cumulants gets improved and higher order cumulants become available.

VIII Conclusions

We have presented results on the equation of state of strong-interaction matter obtained from a sixth order Taylor-expansion of the pressure of (2+1)-flavor QCD with physical light and strange quark masses. We discussed expansions at vanishing strangeness chemical potential $\mu_{S}=0$ as well as for strangeness neutral systems $n_{S}=0$. We have discussed in detail the latter case for a fixed electric charge to net baryon-number ratio, $n_{Q}/n_{B}=0.4$, which is appropriate for situations met in heavy ion collisions. The results, however, can easily be extended to arbitrary ratios of $n_{Q}/n_{B}$. We find that the dependence of basic thermodynamic observables on $n_{Q}/n_{B}$ is small for $0\leq n_{Q}/n_{B}\leq 1/2$. This may be of interest for applications in heavy ion collisions where strong external magnetic fields and non-trivial topology in QCD can lead to charge asymmetries in different regions of phase space.

We have presented a parametrization of basic thermodynamic observables in terms of ratios of fourth order polynomials in the inverse temperature which is appropriate in the temperature range studied here, i.e. $T\in[130\leavevmode\nobreak\ {\rm MeV},330\leavevmode\nobreak\ {\rm MeV}]$.

We presented results for lines of constant pressure, energy and entropy density in the $T$-$\mu_{B}$ plane and showed that corrections of ${\cal O}(\hat{\mu}_{B}^{4})$ are negligible for $\hat{\mu}_{B}<2$. For all three observables the curvature term at ${\cal O}(\hat{\mu}_{B}^{2})$ is smaller than $\kappa^{\rm max}=0.012$. This suggest that, e.g. energy density and pressure, would drop on the crossover line for the chiral transition, if the corresponding curvature coefficient turns out to be larger than $\kappa^{\rm max}$.

The Taylor series for pressure and net baryon-number density as well as energy density and entropy density determined for $\mu_{S}=0$ as well as $n_{S}=0$ have expansion coefficients that are close to HRG model results at low temperature. In general ratios of subsequent expansion coefficients approach the corresponding HRG model values from below when lowering the temperature. As a consequence, in the entire temperature range explored so far, the expansions are ”better behaved” than the HRG model series, which have an infinite radius of convergence. Assuming that the current results obtained with expansion coefficients up to $6^{th}$ order are indicative for the behavior of higher order expansion coefficients and taking into account the current errors on $6^{th}$ order expansion coefficients we concluded that at temperatures $T>135$ MeV the presence of a critical point in the QCD phase diagram for $\mu_{B}\leq 2T$ is unlikely.

Acknowledgments

This work was supported in part through Contract No. DE-SC001270 with the U.S. Department of Energy, through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Nuclear Physics, the DOE funded BEST topical collaboration, the NERSC Exascale Application Program (NESAP), the grant 05P12PBCTA of the German Bundesministerium für Bildung und Forschung, the grant 56268409 of the German Academic Exchange Service (DAAD), grant 283286 of the European Union, the National Natural Science Foundation of China under grant numbers 11535012 and 11521064 and the Early Career Research Award of the Science and Engineering Research Board of the Government of India. Numerical calculations have been made possible through an INCITE grant of USQCD, ALCC grants in 2015 and 2016, and PRACE grants at CINECA, Italy, and the John von Neumann-Institute for Computing (NIC) in Germany. These grants provided access to resources on Titan at ORNL, BlueGene/Q at ALCF and NIC, Cori I and II at NERSC and Marconi at CINECA. Additional numerical calculations have been performed on USQCD GPU and KNL clusters at JLab and Fermilab, as well as GPU clusters at Bielefeld University, Paderborn University, and Indiana University. We furthermore acknowledge the support of NVIDIA through the CUDA Research Center at Bielefeld University.

Appendix A Details on simulation parameters and data sets

Our main data sets have been generated on lattices of size $N_{\sigma}^{3}\times N_{\tau}$, with $N_{\sigma}/N_{\tau}=4$ and $N_{\tau}=6,\ 8$ and $12$. We performed calculations with two different light to strange quark mass ratios, $m_{l}/m_{s}=1/20$ and $1/27$, respectively. The simulation parameters are summarized in Table 3 and Table 4.

Appendix B Constraints on chemical potential for strangeness neutral systems

with fixed electric charge to baryon-number ratio

We are interested in expansion coefficients for strangeness neutral systems in which the net electric-charge is proportional to the net baryon-number. I.e. we introduce the constraint given in Eq. 31. These constraints can be fulfilled order by order in the Taylor expansion of the number densities by choosing the expansion coefficients of the series for $\hat{\mu}_{Q}$ and $\hat{\mu}_{S}$, given in Eq. 14, appropriately, i.e. the coefficients $s_{n}$ and $q_{n}$ can be determined order by order. We start with the Taylor series for the number densities introduced in Eq. 16 and define the expansion coefficients as

[TABLE]

for $n=1,\ 3,\ 5$. At each order in the expansion we then have to solve a set of two linear equations, which always have the same structure. We find as solutions

[TABLE]

and

[TABLE]

At leading order one finds for the terms $m_{1}^{X}$,

[TABLE]

and the contributions to the next-to-leading order expansion terms, $m_{3}^{X}$, are given by

[TABLE]

Finally the contributions to the next-to-next-to-leading order expansion terms, $m_{5}^{X}$, are given by

[TABLE]

In (2+1)-flavor QCD calculations the light ($u,\ d$) quark masses are taken to be degenerate. A consequence of this degeneracy is that not all generalized susceptibilities $\chi_{ijk}^{BQS}$ that enter the above expressions are independent. In a given order $n\equiv 2l\equiv i+j+k$ this results in a set of relations among the expansion coefficients. In general, at order $n=2l$, there are $l(l+1)$ constraints, i.e. for $l=1$ this gives rise to two relations, Bazavov:2012vg

[TABLE]

for $l=2$ there are six constraints,

[TABLE]

and for $l=3$ there are twelve constraints,

[TABLE]

Using these constraints it is tedious, but straightforward, to show that in the isospin symmetric case, $r=1/2$, indeed all expansion coefficients for the electric charge chemical potential vanish, i.e. $\hat{\mu}_{Q}=0$ to all orders in $\mu_{B}$.

We show results for the LO expansion coefficients $s_{1}$ and $q_{1}$ and the ratios of the NLO and LO expansion coefficients, $s_{3}/s_{1}$ and $q_{3}/q_{1}$ in Fig. 17. As can be seen the NLO coefficients are already negligible for T\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}170 MeV. The absolute value of the NNLO expansion coefficients $s_{5}$ and $q_{5}$ never is larger than 1% of the corresponding LO coefficients.

In Fig. 17, we also show results from hadron resonance gas (HRG) model calculations. The black curves are the predictions of the usual HRG model which consists of all the resonances listed in the Particle Data Group Tables up to 2.5 GeV (PDG-HRG). The PDG-HRG results for $s_{1}$ are substantially smaller than the continuum extrapolated lattice QCD results. It has been argued in Bazavov:2014xya that this can be caused by contributions from additional, experimentally not yet observed, strange hadron resonances which are predicted in quark model calculations. A HRG model calculation based on such an extended resonance spectrum (QM-HRG) is also shown in Fig. 17. At finite values of the lattice cut-off we observe significant differences between lattice QCD calculations and both versions of the HRG models. This is in particular the case for the expansion coefficients of the electric charge chemical potentials. One thus may wonder whether these deviations can be understood in terms of taste violations in the staggered fermion formulation which result in a modification of the resonance spectrum and affect most strongly the light pseudo-scalar (pion) sector.

Appendix C The coefficient $\kappa_{4}^{f}$ of lines of constant physics at ${\cal O}(\mu_{B}^{4})$

We will present here results for the expansion coefficient $\kappa_{4}^{f}$ of lines of constant physics defined in Eq. (43),

[TABLE]

The coefficients $f_{2k}$ are defined by

[TABLE]

In particular, we will give explicit expressions for the case of constant pressure ($f\equiv P$), constant energy density ($f\equiv\epsilon$) and constant entropy density($f\equiv s$). For the pressure we had the earlier expression (Eq. (15))

[TABLE]

Comparing Eqs. (67) and (66) we have, $f_{0}=P(T,0)\equiv T^{4}P_{0}$, $f_{2}=T^{4}P_{2}$ and $f_{4}=T^{4}P_{4}$. Thus,

[TABLE]

Here $s$ and $C_{V}$ are the entropy density and specific heat per unit volume at vanishing chemical potential. Similarly,

[TABLE]

Putting everything together we get, for the pressure:

[TABLE]

where $\sigma_{2}$ denotes the ${\cal O}(\hat{\mu}_{B}^{2})$ expansion coefficient of the entropy density as introduced in Eq. 24.

Next we consider $\kappa_{4}^{\epsilon}$. Since the energy density is also of dimension four, we only need to replace $P_{2n}$ with $\epsilon_{2n}$ in the first line of Eq. (70). With this we obtain,

[TABLE]

Since $C_{V}\equiv\left(\partial\epsilon_{0}/\partial T\right)_{\mu_{B}}$, the above may be written as

[TABLE]

Finally we consider $\kappa_{4}^{s}$. Since the entropy density is of dimension three, Eqs. (68) become

[TABLE]

and therefore

[TABLE]

To zeroth order, the specific heat is also given by $C_{V}=\left(\partial(Ts)/\partial T\right)_{\mu_{B}}$. Thus,

[TABLE]