Equation of state of a hot-and-dense quark gluon plasma: lattice simulations at real µ B vs. extrapolations

The equation of state of the quark gluon plasma is a key ingredient of heavy ion phenomenology. In addition to the traditional Taylor method, several novel approximation schemes have been proposed with the aim of calculating it at ﬁnite baryon density. In order to gain a pragmatic understanding of the limits of these schemes, we compare them to direct results at µ B > 0, using reweighting techniques free from an overlap problem. We use 2stout improved staggered fermions with 8 time-slices and cover the entire RHIC BES range in the baryochemical potential, up to µ B /T = 3.

Introduction The equation of state of strongly interacting matter under extreme conditions -high temperatures or baryon densities -plays a key role in many physical systems, such as the early Universe, heavy ion collisions and neutron stars.The most well established firstprinciples method to study the strongly coupled regime is lattice QCD [1], that maps the path integral formulation of QCD to a classical statistical-mechanical system, which can be simulated with Monte Carlo methods.Indeed, many properties of strongly interacting matter at zero baryon density have been elucidated using this method, such as the crossover nature of the transition [2], the value of the transition temperature [3,4] and the equation of state [5,6].Studies at finite baryon density are, however, hampered by the sign problem: the Boltzmann weights in the usual path integral representation become complex, making them unsuitable for importance sampling.Thus, most lattice results on the properties of hot-and-dense QCD matter rely on extrapolations from zero [7][8][9][10][11][12][13][14][15][16][17] or purely imaginary chemical potential [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33], two situations with no sign problem.
In spite of the difficulties, there has recently been considerable progress in the calculation of the equation of state of a hot-and-dense quark gluon plasma (QGP).First, Taylor coefficients of the pressure p in the baryochemical potential µ B have been calculated up to fourth order in the continuum limit [11,12,14] and up to eighth order at finite lattice spacing [17,26,31] -albeit with rather large uncertainties for the sixth and eighth order coefficients.Second, several resummation schemes have been proposed for the Taylor expansion [34][35][36][37][38], with a promise of better convergence properties.However, at the moment there is no theoretical understanding of the convergence properties of these schemes.It is therefore important for phenomenological applications to determine the region of validity of these techniques.This is the purpose of this work.To this end, we use some novel developments in simulation techniques for finite density QCD, concerning reweighting schemes where the pressure difference between the simulated and target theories can be calculated without encountering heavy tailed distributions -i.e., an overlap problem.Recently, studies of this type have become feasible for improved lattice actions with physical quark masses for at least two schemes: sign reweighting, and phase reweighting [39,40].These have the advantage of giving direct, reliable results at µ B > 0 -provided that the exponentially difficult sign problem in the reweighting is dealt with by sufficient statistics.By comparing the equation of state calculated using phase reweighting with that coming from the Taylor expansion and its resummations, we can quantify the systematic bias of the different truncations of the different schemes, giving unprecedented insight into the equation of state of the QGP.
We simulate 2stout improved staggered fermions with physical quark masses on lattices with 8 time slices -a discretization that is often used as the first or even the second point of continuum extrapolations of several thermodynamic quantities [2,5,24,30,[41][42][43][44][45][46].The cut-off effects on the equation of state are moderate (after applying a tree-level improvement [5]).We reimplement all discussed extrapolation schemes using this same setup.Thus, the differences in the predictions can only come from the systematics of the extrapolation.Since we are mostly interested in a comparison of various methods for a statistical physics system, the physical volume can be chosen freely.We use a fixed aspect ratio of LT = 2, with T the temperature and L the spatial size of the box.As we perform no finite volume scaling study, we must restrict ourselves to the study of the QGP equation of state, leaving that of the fate of the crossover transition to a later date.Observables The grand canonical partition function of lattice QCD is given schematically by: where S g (U ) is the gauge action, det M is the quark determinant, implicitly including all flavors, as well as staggered rooting, µ collectively denotes the chemical potentials of all quark flavors and U are the link variables.We work with the quark chemical potentials satisfying µ q ≡ µ u = µ d and µ s = 0. We concentrate on the pressure as a function of the temperature and dimensionless chemical potential μB ≡ µ B /T = 3µ q /T ≡ 3μ q , given by p(T, μB ) ≡ p(T, μB ) In particular, we study the pressure difference between zero and non-zero chemical potentials: We also compute the light quark density: (4) The integral of nL over μB is ∆p.We calculate the equation of state with several different methods.
Reweighting from µ B = 0 In what is arguably the simplest reweighting scheme, simulations are performed at µ B = 0 and ∆p is reconstructed via [47]: While Eq. ( 5) is exact with infinite statistics, the tails of the distribution of the weights det M(μ B ) det M(0) are heavy and so hard to sample correctly (overlap problem), and it is therefore hard to judge the reliability of results obtained with this method.It was proposed that the overlap problem can be mitigated by reweighting in more parameters [48][49][50].However, even with multi-parameter reweighting, the overlap problem remains the main bottleneck of the method (at least around the transition line, where this question was studied in Ref. [51]).We include results from (one-parameter) reweighting from µ B = 0 for completeness, and also because one of the resummation strategies we are going to test can be regarded as a truncation of this reweighting scheme.
Phase reweighting A way to avoid heavy-tailed distributions in the weights is to simulate a theory where these can only take values in a compact domain.Two examples are sign reweighting [39,40,52,53] and phase reweighting [40,54,55].Here we use the latter.Here in the simulated theory one replaces the quark determinant with its absolute value.The simulated ensemble is called the phase quenched ensemble, corresponding to a finite isospin chemical potential, i.e., µ u = −µ d .The pressure then schematically reads: It is also given by the integral of the isospin density: from which the pressure at finite µ B is obtained as: p(T, μB = 3μ q ) − pI (T, μq ) = 1 where e iθ = det M(μ B ) |det M(μ B )| is the complex phase factor of the fermion determinant and . . .P Q means taking an expectation value in the phase quenched theory.One finally arrives at: Alternatively, one can calculate nL directly with reweighting: from which the pressure difference can be calculated via integration.The two methods are not a priori guaranteed to give the same results if the observable in the numerator of Eq. (10), namely e iθ ∂ ∂ μB ln det M, has an overlap problem.In principle this could happen, as only the pressure difference p(T, μB = μq ) − p I (T, μI = μq ) ∝ e iθ P Q in Eq. ( 8) is guaranteed to be free of one -due to the compactness of the observable e iθ .
The phase diagram at finite isospin density was calculated in Ref. [45].The phase quenched ensemble has a pion condensed phase for low T and µ q > m π /2.In this region, the sign problem is expected to be severe.Here we avoid this issue since we concentrate on the equation of state of the QGP, and we find a mild sign problem: e iθ P Q = cos θ P Q never gets below 0.1 in any of our ensembles and is always more than 10σ away from zero.
Taylor expansion The pressure is expanded in powers of the baryochemical potential: p(T, μB ) = p 0 (T ) + p 2 (T ) μ2 B + p 4 (T ) μ4 B + . . . .(11) We calculate the Taylor expansion coefficients in two different ways: by using configurations generated at µ B = 0 to calculate them directly, and by using simulations at imaginary chemical potentials to obtain them from a fit.These are the two procedures used in the literature so far.For a recent example of the first method see, e.g., Ref. [17], for the second one Refs.[26,31].More details on the determination of the coefficients are given in the section on our numerical results.
Resummations based on shifting sigmoid functions A resummation of the Taylor expansion was introduced in Ref. [34], defined by the implicit equation This is motivated by the empirical observation of the existence of an approximate scaling variable T (1 + κμ 2 B ), coming from lattice studies at imaginary chemical potential.These show that certain observables collapse on single (sigmoid shaped) curves when plotted against this variable [33,34].Up to μB = 1.5 the existence of an 2. The direct results for nL at non-zero µB compared with different approximations (or extrapolations): the Taylor method to different orders up to p8 and the exponential resummations [35] to order N = 2, 4 and 6, calculated from the ensemble at µB = 0, as well as the shifting nL/n SBL L method calculated from imaginary chemical potential simulations and reweighting from µB = 0.
approximate scaling variable has also been confirmed directly at a real µ B by reweighting from the sign quenched ensemble [40].We call this the nL (T,μ B ) μB shifting method.This is a systematically improvable expansion.Its validity is not predicated on the existence of this approximate scaling variable, only its fast convergence: i.e., if an approximate scaling variable exists, we expect the method to converge fast.
A disadvantage of the previous scheme is that it was designed to work only near the crossover.To make it more suitable for high T , it was later refined [38] by introducing a Stefan-Boltzmann correction: nL (T, μB ) nSBL where we introduced the Stefan-Boltzmann limit: We call this method the nL (T,μ B ) nSBL L shifting method.Both shifting methods can be implemented by fitting the κ n or λ n coefficients to data at imaginary µ B .
Exponential resummation A different resummation is based on truncating the reweighting from µ B = 0 [35], by approximating det M(μB) det M(0) (top).The direct Taylor data from µB = 0 simulations has smaller errors than the fit to imaginary µB data.This is mainly due to the small volume in our study.For larger volumes, the signal-to-noise ratio of the direct p6 and p8 would be considerably larger.A spline interpolation of the direct results is included to lead the eye.
lattice data for the Taylor coefficients.The procedure introduced in Ref. [35] had the disadvantage of using a biased estimator for the expectation value in Eq. ( 15), by simply exponentiating the stochastic estimators of the D n .This bias was studied in Ref. [36].We remove this bias by calculating the D n exactly for each configuration, using the reduced matrix formalism [56].Note that unlike the Taylor expansion or the resummations based on shifting sigmoids, this scheme is not defined in terms of thermodynamic quantities, but by an approximation of the integrand of the path integral.Thus, the convergence properties of the scheme are not guaranteed to be governed by the analytic properties of the free energy.
Lattice setup and numerical results We used a treelevel Symanzik improved gauge action, and two steps of stout smearing [57] with a smearing parameter ρ = 0.15 on the gauge links entering the quark determinant, with physical quark masses, using the kaon decay constant f K for scale setting (see Ref. [58] for details).We apply a tree-level improvement factor to the pressure, that makes the N τ = 8 2stout results already close to the continuum [5].This corresponds to division by 1.28 at infinite volume, which is the correction we apply for nL and ∆p.
We simulate the phase quenched ensemble for μ2 B = 1.5, 3, 4.5, 6, 7.5 and 9.These simulations are performed without an explicit symmetry breaking term, commonly used in simulations at a finite isospin density [45,59].Instead, we follow the method of Ref. [40].The resulting ensemble corresponds to an isospin chemical potential µ I = µ u = −µ d .We use these ensembles to reweight to a baryochemical potential with µ u = µ d = µ B /3. Determinant ratios are calculated using the reduced matrix formalism [56,60].We performed the calculation of ∆p in the two inequivalent ways: the one with no possible overlap problem described by Eq. ( 9), and by integration of Eq. (10).We obtain compatible results for ∆p with the two methods, as can be seen in Fig. 1 (left) for T = 160 MeV.
We also simulate at µ B = 0 for four temperatures: T /MeV = 140, 150, 160 and 170.This is to perform reweighting from µ B = 0, to calculate the Taylor expansion coefficients, and to perform the exponential resummation.For all of these, we used the reduced matrix formalism of Ref. [56], so obtaining the Taylor coefficients without the use of stochastic estimators and obtaining D n exactly for each configuration, and encountering no bias in the exponential resummation [35,36].
In Fig. 1 (left) we also show the results for ∆p from reweighting from µ B = 0.These results are again in perfect agreement with the two different ways of reweighting from the phase quenched ensemble.This is true for all temperatures in our study.
Thus, three different reweighting procedures -including one with no possible overlap problem -give identical results.This strongly supports the validity of our results for the equation of state.The full set of results for ∆p from phase reweighting are shown in Fig. 1 (right).We can safely use these results to test extrapolation schemes.
Without comparison with the phase reweighted results, we would have no way to guarantee the correctness of reweighting from µ B = 0, due to the uncontrolled systematics of the overlap problem.This problem is also present in the Taylor method, since for a finite ensemble the Taylor coefficients necessarily correspond to the Taylor coefficients of the pressure one would have obtained by reweighting from µ B = 0.By first showing that reweighting from µ B = 0 works for the region μB ≤ 3, we ensure that we are truly testing the convergence properties of the Taylor series, without being limited by an overlap problem in the higher order coefficients.
Comparison with extrapolation schemes To implement the resummation schemes based on shifting sigmoids, we perform simulations at imaginary chemical potentials, for Im μB 16 π = 0, 4, 6, 7, 8, 9, 10 and 12.We work to order κ 4 in Eq. ( 12) and to order λ 4 in Eq. ( 13), using a simplified version of the analysis of Refs.[34,38].The systematic error includes the fit range in imaginary µ B , the ansatz in µ 2 B , and the interpolation of the light quark susceptibility at µ B = 0. We also demonstrate a more straightforward use of the imaginary chemical potential data by performing a second determination of the Taylor expansion coefficients, fitting nL μB with a polynomial of order μ6 B .For the fits we also include We show the comparison of the different reweighting schemes with the direct data in Fig. 3, as a function of the chemical potential at fixed temperatures of T = 160 MeV and 170 MeV, and in Fig. 2, as a function of T at a fixed μB = √ 7.5 and 3 -the two largest values where we have direct data.The Taylor expansion at next-to-leading order -O(µ 4 B ) in the pressure -is not consistent with the direct data, systematically underestimating n L below 150 MeV, and systematically overestimating it above 150 MeV.This is due to a peak in p 4 (T ) slightly above the crossover temperature.Including the next term in the expansion, with the coefficient p 6 (T ), the Taylor method agrees with the direct data up to μB = 3 at T = 160 MeV and up to μB ≈ 1.2 at T = 170 MeV.Including the nextto-next-to-next-to-leading order term p 8 (T ), the expansion agrees with the direct data at all studied temperatures up to μB = 3.
In contrast, the exponential resummation scheme shows bad convergence properties from μ2 B ≈ 4 for all temperatures.While the N = 2 truncation of the scheme remains close to the direct results in the entire range, the higher orders make the agreement better only below this value, but not above.
At large T , the method based on shifting nL /n SBL L outperforms the method of shifting nL /μ B .This is not surprising, as the Stefan-Boltzmann correction was introduced in Ref. [38] as a way to improve the convergence properties of the scheme at high T .
In summary, we can say that both the Taylor expansion to order μ8 B and the resummation based on shifting nL /n SBL L to order λ 4 accurately describe the direct data for the equation of state in the range 0 ≤ μB ≤ 3 -which includes the entire range of the RHIC Beam Energy Scan.Note the faster convergence of the resummed expansion, as the calculation of the coefficient λ 4 only requires the determination of the Taylor coefficients up to order μ6 B .
On the other hand, the shifting nL /μ B method at order κ 4 has a slight systematic discrepancy with the direct data at large T , and exponential resummation shows bad convergence properties in N above μ2 B ≈ 4. Outlook We judged the reliability of different approximation schemes by comparing them with direct non-zero chemical potential results in the range 0 ≤ μB ≤ 3.While this gives a practical answer to the question of which approximation one can trust, a theoretical understanding of the reasons would also be welcome.For the schemes defined purely in terms of thermodynamic quantities, such as the Taylor expansion or the resummations based on shifting sigmoids, this requires knowledge of the position of partition function (Lee-Yang) zeros in the complex µ B plane [16,37,[60][61][62].The exponential resummation scheme, instead, is not defined in terms of thermodynamic quantities, but comes rather from manipulating the integrand of the path integral for the partition function.Understanding its convergence region might also require better understanding of the nuances of the path integral, in addition to the thermodynamic singularities.We speculate that the limited convergence region has to do with quark determinant zeros.In fact, the sum appearing as the argument of the exponential in Eq. ( 15) approximates the effective action of the quarks in a fixed gauge field background, with a radius of convergence determined by the determinant zeros, which correspond to logarithmic divergences of the effective action.These are not simply related to the Lee-Yang zeros of the partition function, and may provide stronger limitations on the convergence of the expansion.
An obvious challenge is to extend the range of validity of the methods studied in this paper to lower T and higher μB , so that the transition line [15,30,33,63,64], and the location of the conjectured critical endpoint [65][66][67][68][69] can be studied with first principle lattice calculations.Of course, the continuum and infinite volume limits will also have to be taken eventually.

FIG. 1 .
FIG. 1. Left panel: The pressure difference between the zero and finite chemical potential theory calculated with several different reweighting methods as a function of the quark chemical potential at T = 160 MeV.Results at finite isospin chemical potential are shown as purple points.Results reweighted from the phase quenched ensemble to finite baryon density with the two methods discussed in the text are shown by red and blue points.Results from reweighting from zero chemical potential are shown in grey.Right panel: The pressure difference between zero and non-zero baryochemical potentials, calculated with phase reweighting as a function of the temperature, for different fixed values of μ.

µ 2 B 4 B 6 B
FIG.2.The direct results for nL at non-zero µB compared with different approximations (or extrapolations): the Taylor method to different orders up to p8 and the exponential resummations[35] to order N = 2, 4 and 6, calculated from the ensemble at µB = 0, as well as the shifting nL/n SBL

BFIG. 3 .
FIG. 3. The density nL as a function of the temperature T for μB = √ 7.5 (left) and μB = 3 (right) for the ordinary Taylor expansion (bottom) and the resummation schemes based on shifting nL/μB and nL/n SBL L

d 2 p dμ 2 B and d 4 p dμ 4 B
at µ B = 0 as further datapoints.