Lattice QCD equation of state at finite chemical potential from an alternative expansion scheme
S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, K. K. Szabo
LLattice QCD equation of state at finite chemical potential from an alternativeexpansion scheme
S. Bors´anyi, Z. Fodor,
2, 1, 3, 4
J. N. Guenther, R. Kara, S. D.Katz, P. Parotto, ∗ A. P´asztor, C. Ratti, and K. K. Szab´o
1, 4 University of Wuppertal, Department of Physics, Wuppertal D-42119, Germany Pennsylvania State University, Department of Physics, State College, PA 16801, USA Inst. for Theoretical Physics, ELTE E¨otv¨os Lor´and University,P´azm´any P. s´et´any 1/A, H-1117 Budapest, Hungary J¨ulich Supercomputing Centre, Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany Aix Marseille Univ., Universit´e de Toulon, CNRS, CPT, Marseille, France E¨otv¨os University, Budapest 1117, Hungary Department of Physics, University of Houston, Houston, TX 77204, USA (Dated: February 15, 2021)Taylor expansion of the equation of state of QCD suffers from shortcomings at chemical potentials µ B ≥ (2 − . T . First, one faces difficulties inherent in performing such an expansion with alimited number of coefficients; second, higher order coefficients determined from lattice calculationssuffer from a poor signal-to-noise ratio. In this work, we present a novel scheme for extrapolatingthe equation of state of QCD to finite, real chemical potential that can extend its reach furtherthan previous methods. We present continuum extrapolated lattice results for the new expansioncoefficients and show the thermodynamic observables up to µ B /T ≤ . I. INTRODUCTION
The phase diagram of Quantum Chromodynamics(QCD) is an open field of investigation, which is cur-rently at the center of intense efforts from the theoret-ical and experimental communities alike. At vanishingbaryon density, the transition between confined and de-confined matter is known to be an analytic crossover [1].This knowledge was gained through lattice simulations,which represent a systematically improvable method tosolve equilibrium QCD. Although at finite baryon densitylattice QCD faces a sign problem, numerous results havebeen published for moderate chemical potentials in recentyears [2, 3]. New techniques that allow for direct simula-tions at finite chemical potential in the presence of a signproblem are the subject of intense investigation. Promis-ing results are available from Lefschetz thimbles [4–6], theComplex Langevin equation [7–9] or reweighting-basedmethods [10]. However, these novel approaches cannotbe applied to large scale QCD simulations with physicalquark masses yet.There are many indirect lattice methods to study QCDat finite density. These are based on the known analyticfeature of the QCD free energy at zero baryo-chemicalpotential ( µ B ). The conceptually simplest method is theTaylor expansion, where the leading µ B -derivatives of therelevant observables are calculated [11–14]. These deriva-tives were often calculated for chiral observables, and the µ B dependence of the transition temperature could beextracted [15–18].An additional twist to the Taylor method is the ob-servation that these coefficients can be efficiently calcu- ∗ Corresponding author: [email protected] lated by simulating at imaginary values of the chemicalpotential(s), besides µ B = 0 [19, 20]. The use of imag-inary chemical potentials is motivated by the analyticcrossover at µ B = 0, from which the smooth behavior ofthe thermodynamic observables as a function of µ B fol-lows. This approach has been popular for more than adecade [17, 21–26]. This possibility opened an avenuetowards finite density physics, which is often referredto as analytical continuation. This name suggests that,besides the computation of the Taylor expansion coeffi-cients, other extrapolation schemes can be established.The success of this method was most visible in the studyof the QCD transition line, where continuum extrapo-lated results are available for the leading µ B dependence[17, 27–29], and recently also for the next-to-leading co-efficient [30]. On the other hand, the Taylor series is justone of the possible expansion or continuation schemes.For example, the Pad´e summation was also considered inthe context of QCD thermodynamics [31–35] .The knowledge of the features of the QCD phase dia-gram from lattice simulations is currently limited to smallchemical potentials and data are mostly available onlyin the transition region. We have to mention, though,that at higher temperatures resummed perturbation the-ory has provided a quantitative description of the chem-ical potential dependence of several observables [36–38].Dedicated lattice studies have bridged the gap betweenthe transition region and perturbative temperatures andfound perfect agreement [39, 40].On the experimental side, several heavy-ion collisionprograms are in place, with the explicit intent of mappingout the phase structure of strongly interacting matter.An important tool for the theoretical interpretation of theresults from these experiments are hydrodynamic simula-tions, which describe the evolution of the system createdin the collisions. These simulations need the equation of a r X i v : . [ h e p - l a t ] F e b µ B /T = 1 χ B ( T ) T [MeV] N LONLOLO µ B /T = 2 χ B ( T ) T [MeV] N LONLOLO µ B /T = 3 χ B ( T ) T [MeV] N LONLOLO
FIG. 1. Baryon density from a Taylor expansion with the coefficients in Ref. [20], at µ B /T = 1 , ,
3, as a function of thetemperature. Different colors correspond to the order to which the expansion is carried out. state of QCD as an input, in the whole range of tempera-tures and densities covered in the experiments. Recently,a Bayesian analysis based on a systematic comparison be-tween heavy-ion data and theoretical predictions showedthat the posterior distribution over possible equations ofstates is compatible with the one calculated on the lat-tice [41]. For this reason, the equation of state at finitedensity is a crucial ingredient to understand the natureof strongly interacting matter, and to support the exper-imental program.The equation of state at vanishing chemical potentialhas been known now for several years over a broad rangeof temperatures [42–44]. The first continuum extrapo-lated extension to finite µ B using the Taylor methodin Ref. [13] was followed by several works with the in-tent of extending these results to higher µ B by addingmore terms in the Taylor series [14, 45]. Currently, eventhe eighth µ B -derivative of the QCD pressure is availablewith modest precision from lattice simulations [20, 46].Very recently, similar results were found by solving aQCD-assisted effective theory with functional methods[47].In this work, we propose a new extrapolation schemeto finite density QCD. We intend to remedy some short-comings of the Taylor-based equation of state, e.g. theextrapolation through a crossover boundary. We discussthis issue and suggest a solution in Section II. The formal-ism for our method is worked out in Section III. We thencompute the new observables from lattice simulations inSection IV and perform their continuum estimate. Usingthis lattice input, we construct the finite density thermo-dynamic functions in Section V. We conclude and discussfurther aspects in Section VI. II. MOTIVATION
The knowledge of the equation of state from latticesimulations commonly consists of the established µ B = 0result [43, 44] and the Taylor expansion coefficients of the pressure around µ B = 0 p ( T, µ B ) T = (cid:88) n =0 n )! χ B n ( T, (cid:16) µ B T (cid:17) n , (1)where χ Bj are the j -th derivatives of the normalized pres-sure: χ Bj ( T, µ B ) = (cid:18) ∂∂µ B /T (cid:19) j p ( T, µ B ) T . (2)Besides diagonal coefficients, one can also define off-diagonal correlators between different conserved chargesin QCD. Correlators between baryon number andstrangeness, which we will need in our procedure, aredefined as follows χ BSjk ( T, µ B ) = (cid:18) ∂∂µ B /T (cid:19) j (cid:18) ∂∂µ S /T (cid:19) k p ( T, µ B ) T . (3)Such correlators have phenomenological relevance [48]and they can also be used to extrapolate the equationof state of QCD in the full, four-dimensional phase dia-gram at finite T, µ B , µ S , µ Q [49, 50]. We will use theˆ µ i = µ i /T shorthand notation in this manuscript.Currently, results for the expansion coefficients areavailable up to order O ( µ B ) [20, 46]. The region of va-lidity of the resulting expansion is usually determined bythe range in chemical potential within which an appar-ent convergence is achieved with the available coefficients.This is stated to be µ B /T < ∼ − . µ B -dependence of the crossover temperature may explain thebasic structure.This may explain why including one more term ina truncated Taylor series will not always improve theconvergence. On the contrary, pathological behavior –namely non-monotonicity in the T - or µ B -dependence –appears in the extrapolated thermodynamic quantitiesat chemical potentials beyond µ B /T < ∼ − .
5. Thisis due to the fact that, for large enough µ B /T , the ob-servables at finite chemical potential are dictated by the µ B = 0 temperature dependence of the last coefficient in-cluded in the expansion. Hence, the structures appearingaround the QCD transition temperature in higher ordercoefficients are “translated” into the finite- µ B behaviorof e.g., the entropy, baryon density, etc.Another inherent problem with the Taylor expansionis the fact that it is carried out at constant temperature.This means that the values of the coefficients at µ B = 0and a certain temperature T , determine the equation ofstate at the same T at finite µ B , while the pseudo-criticaltemperature T pc might have varied considerably. While asufficiently large number of expansion coefficients wouldlead to smooth extrapolated functions, even though theTaylor coefficients themselves present a complex struc-ture around the transition temperature, the problem hereis rather practical. A scheme that could work with fewercoefficients would be much preferable from the numericalcost point of view.In Fig. 1 we show the baryon density n B ( T ) ob-tained from a Taylor expansion with the coefficients inRef. [20], at ˆ µ B = 1 , ,
3. The extrapolation is shownincluding an increasing number of coefficients, to showthe effect of higher-order ones. The leading-order (LO)and higher truncations refer to ∼ ˆ µ B ∂n B ( T ) /∂ ˆ µ B , or ∼ ˆ µ B ∂ n B ( T ) /∂ ˆ µ B , etc. being the last term in theexpansion. The derivatives are taken at µ B = 0.While at ˆ µ B = 1 apparent convergence is achieved atthe NLO level, for higher chemical potential this is notthe case. Especially at ˆ µ B = 3, the inclusion of all thecoefficients in Ref. [20] causes unphysical non-monotonicbehavior. Ultimately, a pathological behaviour could alsocome from a finite radius of convergence. Incidentally, re-cent estimates on coarse lattices [52, 53], but also univer-sality arguments [54] place the convergence in the sameball-park in µ B .In this work, we present an alternative summationscheme which can better cope with the fact that the QCDtransition temperature presents a µ B -dependence.We start from the observation that we made whileworking with imaginary values of the chemical poten-tials in an earlier work on analytical continuation. In theupper panel of Fig. 2 we show temperature scans of thequantity n B ( T ) / ˆ µ B = χ B ( T, ˆ µ B ) / ˆ µ B for several fixedimaginary µ B /T ratios. The 0 / µ B = 0 can beeasily resolved and equals χ B ( T ).The T -dependence of the normalized baryon den-sity at finite chemical potential appears to be simplyshifted/rescaled towards higher temperatures from the µ B = 0 results for χ B . This behavior is more apparentin the vicinity of the transition, where the slope of these T χ B / µ B ( T ) T [MeV] µ B = i 6 π /8 µ B = i 5 π /8 µ B = i 4 π /8 µ B = i 3 π /8 µ B = i 2 π /8 µ B = i π /8 µ B = 0 T rescaled using κ =0.0205 T χ B / µ B ( T ) T (1+ κ ( µ B /T) ) [MeV] µ B = i 6 π /8 µ B = i 5 π /8 µ B = i 4 π /8 µ B = i 3 π /8 µ B = i 2 π /8 µ B = i π /8 µ B = 0 FIG. 2. Upper panel: The (imaginary) baryon density atsimulated (imaginary) baryon chemical potentials, divided bythe chemical potential. The points at µ B = 0 (black) showthe second baryon susceptibility χ B ( T ). Lower panel: samecurves as in the upper panel, with a temperature rescaled inaccordance to Eq. (5) with κ = 0 . curves is larger. At very large, as well as at very low tem-peratures a simple shift cannot describe the physics, sincethe curves become extremely flat. A simple rescaling oftemperatures can be described as: χ B ( T, ˆ µ B )ˆ µ B = χ B ( T (cid:48) , , (4)where the actual temperature difference can be expressedthrough a µ B -dependent rescaling factor that we write forsimplicity as T (cid:48) = T (cid:0) κ ˆ µ B (cid:1) . (5)In the lower panel of Fig. 2 we show a version of the curvesin the upper panel, where all the finite-ˆ µ B curves havebeen shifted in accordance to Eq. (5) with κ = 0 . T -independent parameter governing the trans-formation. - T χ S / µ B ( T ) T [MeV] µ B = i 6 π /8 µ B = i 5 π /8 µ B = i 4 π /8 µ B = i 3 π /8 µ B = i 2 π /8 µ B = i π /8 µ B = 0 χ S ( T ) T [MeV] µ B = i 6 π /8 µ B = i 5 π /8 µ B = i 4 π /8 µ B = i 3 π /8 µ B = i 2 π /8 µ B = i π /8 µ B = 0 FIG. 3. The (imaginary) strangeness density divided bythe baryon chemical potential (upper panel) and the secondstrangeness susceptibility (lower panel) at simulated (imag-inary) baryon chemical potentials. The points at µ B = 0(black) show the baryon- strangeness correlator χ BS ( T ) (up-per panel) and the second strangeness susceptibility χ S ( T )(lower panel), respectively. A similar behavior is observed for other quantities too.We show in Fig. 3 the first and second order fluctuationsof strangeness at imaginary baryon chemical potentials.In analogy with Eq. (4) one has: χ S ( T, ˆ µ B )ˆ µ B = χ BS ( T (cid:48) , , (6) χ S ( T, ˆ µ B ) = χ S ( T (cid:48) , , where T (cid:48) is defined analogously to Eq. (5), albeit withdifferent κ parameters.In order to motivate our alternative summationscheme, let us first consider a crude approximation thatwe will later refine. We take Eq. (5) at face valueand use it together with Eq. (4) to obtain a well de-fined χ B ( T, ˆ µ B ) function. We need a χ B ( T,
0) functionas well, which we borrow from a deliberately simple fit f ( T ) = a + b arctan( c ( T − d )) to our data on a 48 ×
12 lat-tice. In principle, we could not only calculate χ B ( T, ˆ µ B ) at arbitrary ˆ µ B but, blindly believing Eq. (5), one couldcalculate the higher µ B -derivatives as well. While thiswill not describe Nature precisely, it can serve as a test forthe Taylor expansion method. To this end, we took sev-eral ˆ µ B -derivatives of our χ B ( T, ˆ µ B ) function and calcu-lated its truncated Taylor series for the lowest four orders.Here LO means just plotting χ B ( T, µ B = 1 ,
2, the summation up to LO and NLO issufficient to perfectly reproduce the function. However,as the chemical potential is increased, the Taylor expan-sion carried as far as the NNNLO does not reproducethe original function. On the one hand, convergence isachieved more slowly; on the other hand, spurious ef-fects appear, which generally manifest themselves in anon-monotonicity of the function. These spurious effectsin truncated Taylor series were pointed out also in Refs.[55–57].The picture emerging from this simple analysis israther suggestive, especially when compared to the re-sults shown in Fig. 1 (right panel) obtained from actuallattice data. We also note that this simple analysis doesnot suffer from the additional complications of signal ex-traction for higher order expansion coefficients, which inturn play a relevant role in the real-data analysis.
III. FORMALISM
At vanishing chemical potential, it is possible to ex-press the normalized baryon density as a Taylor expan-sion: χ B ˆ µ B ( T, ˆ µ B ) = χ B ( T,
0) + ˆ µ B χ B ( T,
0) + ˆ µ B χ B ( T,
0) + · · · (7)As we saw in Fig. 2, the behavior of χ B ˆ µ B ( T, ˆ µ B ) at finitechemical potential clearly resembles that of χ B ( T, ˆ µ B ),although shifted/rescaled in temperature. As long as χ B / ˆ µ B is a monotonic function of T , the finite densityphysics can be encoded into the T (cid:48) ( T, ˆ µ B ) function. Astraightforward, but systematic generalization of Eq. (5)reads: T (cid:48) ( T, ˆ µ B ) = T (cid:0) κ BB ( T )ˆ µ B + κ BB ( T )ˆ µ B + O (ˆ µ B ) (cid:1) . (8)In the above equation, we introduced the new parameters κ BB ( T ) and κ BB ( T ), which describe the shift/rescalingof the temperature of χ B / ˆ µ B at finite µ B . Analogousparameters will be introduced below to for the case of χ S / ˆ µ B ( κ BS and κ BS ) and of χ S ( κ SS and κ SS ) at finite µ B . In a way, this formalism replaces the fixed temper-ature µ B expansion by a fixed-observable temperatureexpansion.Having now two expressions, Eqs. (7) and (4), for thesame quantity we require their equality at each order in µ B /T = 1 T χ B / µ B ( T ) T [MeV]N LON LONLOLOfull µ B /T = 2 T χ B / µ B ( T ) T [MeV]N LON LONLOLOfull µ B /T = 3 T χ B / µ B ( T ) T [MeV]N LON LONLOLOfull
FIG. 4. Benchmarking various orders of the Taylor method assuming an equation of state, where the µ B -dependence of χ B / ˆ µ B consists of a simple shift in temperature. This equation of state is a somewhat simplified form of the observed behaviour. the ˆ µ B expansion at µ B = 0, having: χ B ( T ) = 6 T κ BB ( T ) dχ dT , (9) χ B ( T ) = 60 T ( κ BB ) ( T ) d χ dT + 120 T κ BB ( T ) dχ dT , which in turn yields: κ BB ( T ) = 16 T χ B ( T ) χ B (cid:48) ( T ) , (10) κ BB ( T ) = 1360 χ B (cid:48) ( T ) (cid:16) χ B (cid:48) ( T ) χ B ( T ) − χ B (cid:48)(cid:48) ( T ) χ B ( T ) (cid:17) . A similar treatment can be applied to the other observ-ables. For the second order fluctuations including baryonnumber and strangeness, one can consider: χ S ˆ µ B ( T, ˆ µ B ) = χ BS ( T, µ B χ BS ( T, µ B χ BS ( T, · · · (11)and: χ S ( T, ˆ µ B ) = χ S ( T, µ B χ BS ( T, µ B χ BS ( T, · · · (12)Similarly as before, one can show that: κ BS ( T ) = 16 T χ BS ( T ) χ BS (cid:48) ( T ) , (13) κ BS ( T ) = 1360 χ BS (cid:48) ( T ) (cid:16) χ BS (cid:48) ( T ) χ BS ( T ) − χ BS (cid:48)(cid:48) ( T ) χ BS ( T ) (cid:17) , and: κ SS ( T ) = 12 T χ BS ( T ) χ S (cid:48) ( T ) , (14) κ SS ( T ) = 124 χ S (cid:48) ( T ) (cid:16) χ S (cid:48) ( T ) χ BS ( T ) − χ S (cid:48)(cid:48) ( T ) χ BS ( T ) (cid:17) . κ x8 lattice χ B4 /60.02 T d χ B2 /dT FIG. 5. Evaluation of Eq. (10) on a coarse lattice with highstatistics. The resulting κ BB ( T ) shows a very mild tempera-ture dependence in the transition region. Before we embark into the discussion of the latticeanalysis, let us have an impression on the discussed quan-tities. Eq. (10) gives a way to directly determine κ BB ( T )and κ BB ( T ) using only µ B = 0 data.This approach might be subject to numerical problems,especially in the case of κ BB ( T ), which is obtained asthe difference of two competing terms. Notice, too, thatEq. (10) contains temperature-derivatives of the χ ( T ) co-efficients, which may pose a numerical challenge, unlessthe coefficients are known at sufficient statistics and res-olution in T .On lattices where high statistics data taking is feasible,we can investigate Eq. (10), at least for κ BB . In the toppanel of Fig. 5 we compare the numerator and (rescaled)denominator of Eq. (10), while their ratio κ BB ( T ) isshown in the bottom panel. In the entire transition re-gion the ratio is consistent with a constant, because thepeak in χ B ( T ) is replicated in the temperature depen-dence of χ B ( T ). As opposed to the Taylor coefficients, κ BB ( T ) shows a very mild temperature dependence.Finally we remark that very similar equations havebeen already used in Ref. [14] to calculate “lines of con-stant physics” to O ( µ B ) order. In this reference the pres-sure, energy density and entropy were calculated usingthe Taylor method, and in a further step lines were drawnon the µ B − T phase diagram, where these quantities areconstant in some normalization. The obtained κ coeffi-cients are closely related to ours.Contrary to Ref. [14], we use Eq. (8) as the definitionof a truncation scheme rather than to investigate a Tay-lor expanded result. In a way, we work in the oppositedirection: we will use lattice data at zero and imaginary µ B to obtain the coefficients in Eq. (8), which then canbe used to either calculate the Taylor coefficients or, evenbetter, to extrapolate the equation of state at finite µ B with no reference to the Taylor coefficients themselves.Also, we used imaginary chemical potentials notonly to calculate the coefficients, but also to studythe single observable first, on which the analysis isbased. We base our description of the entire chemicalpotential-dependence of the QCD free energy functionon χ B ( T, ˆ µ B ) / ˆ µ B . It is essential to rely on one observ-able only, in order to guarantee thermodynamic consis-tency: entropy, pressure and energy density will obey theknown thermodynamic relations (see Section V) only ifthey come from the same truncation scheme. IV. LATTICE DETERMINATION OF THEEXPANSION COEFFICIENTSA. Lattice details
In this work, we use the lattice action and the parame-ters described in Ref. [39]. The action benefits from tree-level Symanzik improvement in the gauge sector and fourlevels of stout smearing for the staggered flavors. The upand down quarks are degenerate. The resulting light pairof quarks, as well as the strange and charm quarks as-sume their respective physical mass.We performed simulations at µ B = 0 on 32 × ×
10, 48 ×
12 and 64 ×
16 lattices in a tempera-ture range of 130 −
300 MeV, and up to 500 MeV onlarger volumes. These simulations were complementedat imaginary values of the chemical potential in the tem-perature range 135 −
245 MeV for the lattice resolutions N τ = 8 , . . . ,
12. The µ B (cid:54) = 0 data were simulated at µ S = 0, some of these ensembles were already used inRef. [20].In addition, we performed a high-statistics run on thecheap and coarse 24 × B. The coefficients κ ij and κ ij For the determination of κ ij and κ ij , one can takeadvantage of simulations both at zero and finite imag-inary chemical potential. Using Eq. (10) we calculated -0.01 0 0.01 0.02 0.03 0.04 0.05 140 160 180 200 220 240 κ n BB ( T ) T [MeV] κ cont. est. κ cont. est. κ HRG κ HRG -0.01 0 0.01 0.02 0.03 0.04 140 160 180 200 220 240 κ n BS ( T ) T [MeV] κ cont. est. κ cont. est. κ HRG κ HRG -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 140 160 180 200 220 240 κ n SS ( T ) T [MeV] κ cont. est. κ cont. est. κ HRG κ HRG
FIG. 6. Continuum extrapolated result for the parameters κ BB ( T ) and κ BB ( T ) (top panel), κ BS ( T ) and κ BS ( T ) (cen-tral panel), and κ SS ( T ) and κ SS ( T ) (bottom panel). Theparameters κ ij are shown in blue, and the κ ij in red. HRGresults for all quantities are shown up to T = 160 MeV (ingreen for κ ij and orange for κ ij , respectively). κ BB ( T ), shown in Fig. 5. To extract κ BB ( T ) using thesame strategy, a precise result on χ B ( T ) would be neces-sary.Instead, we utilize imaginary chemical potential simu-lations as follows. We perform simulations at imaginaryvalues of the baryon chemical potential:ˆ µ B = i nπ , n ∈ { , , , } (15)with ˆ µ Q = ˆ µ S = 0.We simulate temperatures in the range T = 135 −
245 MeV. For each temperature T and chemical potentialˆ µ B we determine the temperature T (cid:48) for which Eq. (7)holds, hence defining a function T (cid:48) ( T, ˆ µ B ). Rearrangingthe terms in Eq (8), and in a similar way for the otherobservables, one can write:Π( T, ˆ µ B ) = κ ij ( T (cid:48) ) + κ ij ( T (cid:48) )ˆ µ B + κ ij ( T (cid:48) )ˆ µ B + . . . (16)with the proxy quantityΠ( T, ˆ µ B ) = T (cid:48) ( T, ˆ µ B ) − TT (cid:48) ( T, ˆ µ B )ˆ µ B . (17)The proxy Π( T, ˆ µ B ) can be determined using latticesimulations. The relatively precisely known function χ B ( T,
0) is first interpolated in temperature. Afterwards,we only need to measure χ B ( T, µ B ) / ˆ µ B on an ensemblewith imaginary µ B . Equating this to χ B ( T (cid:48) ,
0) gives us T (cid:48) ( T, µ B ), while we have to take care of the propagationof the statistical errors.Having determined Π( T, ˆ µ B ) for several imaginarychemical potentials and several lattice spacings for eachgiven temperature, one can perform a polynomial fitin ˆ µ B and obtain the expansion coefficients. This canbe done separately for each lattice. However, we pre-fer to combine the ˆ µ B and continuum fits in one two-dimensional fitting procedure. This combined fit is re-peated for every temperature, in steps of 5 MeV.In order to estimate the systematic uncertainties as-sociated to our results, we perform a number of analy-ses at each temperature. There are several ambiguouspoints that need to be considered. Most obviously, onecould choose to include the κ ij ( T ) term in the fit or not,or consider the fit of 1 / Π instead. Also the range inimaginary µ B is arbitrary: we consider Im µ B ≤ . µ B ≤ .
4. When different lattice spacings are fittedtogether in a continuum extrapolation, one selects thebare parameters such that the ensembles correspond tothe same physical temperature. This choice is, however,ambiguous too, as the scale setting may be based on vari-ous observables. In our case, we consider f π or w to thispurpose. As we mentioned before, χ B ( T,
0) is subject toan interpolation, performed through basis splines. Thesame is true for χ B ( T, ˆ µ B ) at finite chemical potentials.Since the location of the node points is also arbitrary, weinclude three versions at µ B = 0 and two at imaginary µ B , each with perfect fit quality. Finally, in the contin-uum extrapolation the coarsest lattice, 32 ×
8, has eitherbeen used or dropped. The listed options can be consid-ered in arbitrary combinations. In total we carry out all144 fits to perform a continuum extrapolation of κ ij ( T )and κ ij ( T ). After dropping the fits with a Q-value belowa percent, we use uniform weights to produce histogram out of these (somewhat less than) 144 results for eachtemperature. The width of the histogram defines thesystematic error. In the plots we show combined errors,where we assume that statistical and systematic errorsadd up in quadrature. This systematic error estimationprocedure has been used and described in more detail inseveral of our works, most recently in Ref. [20].In the top panel of Fig. 6 we show the results of thetemperature-by-temperature fit procedure for the param-eters κ BB ( T ) and κ BB ( T ), alongside the correspondingexpectations from the Hadron Resonance Gas (HRG)model. We find that, within errors, κ BB ( T ) has hardlyany dependence on the temperature, while κ BB ( T ) is ev-erywhere consistent with zero at our current level of pre-cision. Nonetheless, a clear separation of almost one or-der of magnitude appears between these two coefficients.We also note that good agreement with the HRG resultsis found up to at least T = 160 MeV.In the central and bottom panels of Fig. 6 we showour results for the parameters κ BS , κ BS and κ SS , κ SS respectively, together with their HRG determinations.While for χ BS barely any temperature dependence is ob-served, as in the case of χ B , for χ S a much stronger T -dependence is clearly visible. As in the case of κ BB , κ BS is consistent with zero throughout the whole tem-perature range we consider. On the other hand, κ SS risesabove zero for temperatures T > ∼
220 MeV.Note that the error bars in these plots are highly cor-related. The correlation is mostly systematic. The ap-parent ‘waves’ are often statistically not as significant asit seems at a first glance. For comparison we refer tothe direct result on our coarsest lattice in Fig. 5, wherethese ‘waves’ are absent. Before moving on to calculatethermodynamic observables at finite real chemical poten-tial, we construct a smoother version of our final resultsfor κ ij and κ ij , in order to limit the influence of numer-ical effects on the final observables, but also to facilitatethe temperature derivatives that enter the entropy for-mula. Thanks to the very mild dependence of the co-efficients on the temperature, we perform a polynomialfit of order 5 for the κ ij , and of order 2 for the evenless T -dependent κ ij . The very good fit qualities did notmotivate higher order polynomials, using fourth or sixthorder for κ ij hardly changes the result. In order to sta-bilize the low-temperature behavior, we include in the fittwo points from the HRG model for T = 120 ,
130 MeV,to which we associate an uncertainty of 5% for κ ij andof 300% for κ ij . The choice of these particular values forthe uncertainties is uniquely guided by the necessity ofplacing a constraint on the low-T behavior, while avoid-ing to drive the fit too strongly. For this reason, thesearbitrary errors are chosen to be smaller, but compara-ble to the lattice ones. We note that the fits we performtake fully into account the correlations between results atdifferent temperatures, systematic as well as statistical.Thus we encode all errors into the (correlated) errors ofthe coefficients of a polynomial.In Fig. 7 we show the results of the fits in darker color, -0.02-0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 120 140 160 180 200 220 240 κ n BB ( T ) T [MeV] κ cont. est. κ cont. est. κ fit κ fit κ HRG κ HRG -0.02-0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 120 140 160 180 200 220 240 κ n BS ( T ) T [MeV] κ cont. est. κ cont. est. κ fit κ fit κ HRG κ HRG -0.02-0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 120 140 160 180 200 220 240 κ n SS ( T ) T [MeV] κ cont. est. κ cont. est. κ fit κ fit κ HRG κ HRG
FIG. 7. Results of the polynomial fits to κ BB and κ BB (toppanel), κ BS ( T ) and κ BS ( T ) (central panel), and κ SS ( T ) and κ SS ( T ) (bottom panel). The parameters κ ij are shown inblue, and the κ ij in red. The fitted quantities from Fig. 6are shown by lighter blue and red points. Due to the lack ofpoints at low T , we continued the data set with HRG points,show as green and orange dots. Though the polynomials donot always follow all “excursions” of the T -by- T result, thereduced χ of the correlated temperature fit is always < T [MeV] χ T d χ /dT FIG. 8. Continuum extrapolated χ B ( T ) and T dχ B ( T ) /dT functions from N τ = 10, 12 and 16 lattices at µ B = 0. To-gether with κ BBn ( T ) these functions enter the thermodynamicanalysis. with the fitted data in lighter shades. The HRG pointsincluded in the fit are shown as well. C. Continuum result on χ B ( T ) and its temperaturederivative In order to determine the value of thermodynamicquantities at finite real chemical potential, a continuumresult for the basic quantity χ B ( T,
0) is required in addi-tion to κ BB ( T ) and κ BB ( T ). For some observables likethe entropy, the temperature derivative of χ B ( T ) has tobe calculated as well. We have already published thisquantity in Ref. [39], but not its temperature derivatives.To make this work self-contained, we update this deter-mination with the updated statistics using lattices up to64 × T ∈ [130 MeV ,
300 MeV] data points. For the quantity
T dχ B ( T ) /dT we extended the data range with the HRGcurve, so that the numerical derivative has a level armat the lowest temperatures. The basis splines are cubicsplines in the given temperature range. We fit the splinesin temperature and in 1 /N τ in one step. The fitted func-tion in this range is thus: χ B ( T, N τ ) = n (cid:88) i =1 α i b i ( T ) + 1 N τ n (cid:88) i =1 β i b i ( T ) , (18)where b i ( t j ) = δ ij and t j are the knots for the spline with j = 1 . . . n . n and the t j set are the only arbitrary in-puts. We combine the result from four sets of values forthem, so that we can estimate the systematics. Only thelattices with N τ = 10 ,
12 and 16 enter the continuum ex-trapolation. For χ B ( T, N τ ) in Eq. (18), the tree-levelimprovement was already applied [39]. The lattice reso-lutions and the spline knots are selected such that the fitqualities (Q value) are distributed as expected. We alsotook into account the scale setting ambiguity by combin-ing the results with w − and f π − based scale settings.The systematics of the T -derivative T dχ B ( T, /dT andof χ B ( T,
0) were determined separately.In the high temperature regime, the smooth monotonicbehavior is not well described with cubic splines. Insteadwe performed a high order polynomial fit in the inversetemperature 1 /T . The known analytical form of the con-vergence to the Stefan-Boltzmann limit is, of course, notthis polynomial. We do not wish to enforce the perturba-tive behavior at the intermediate temperatures that wedescribe. Thus, the constant in the 1 /T description is notexactly the Stefan-Boltzmann limit, and for this reason,we cannot regard this as a basis for an extrapolation.However, this simple approach allows the interpolationand the calculation of the derivative. We used latticedata in the range T = 180 −
450 MeV.High- and low-temperature regions overlap and thetwo fitting methods give consistent results between T =200 −
280 MeV for both quantities. Thus, we simply con-catenate the resulting functions at T = 260 MeV. Weshow the final version of the T dχ B ( T, /dT and χ B ( T, V. THERMODYNAMICS AT REAL CHEMICALPOTENTIAL
Once n B is determined, we have everything we need toextract the other thermodynamic quantities. The inte-gration constant for the pressure is obviously the pressureitself at µ B = 0.We note here that, on the lattice, we always deal withdimensionless thermodynamic quantities, which corre-spond to the physical ones divided by suitable powers ofthe temperature. E.g., the dimensionful baryon densityis n B = T ˆ n B (we will hereafter use the hat to indicatedimensionless quantities).From the baryon density ˆ n B (ˆ µ B , T ), the pressure isobtained through simple integration: p ( µ B , T ) T = ˆ p (ˆ µ B , T ) = ˆ p (0 , T ) + (cid:90) ˆ µ B dˆ µ (cid:48) B ˆ n B (ˆ µ (cid:48) B , T ) . (19)The entropy density is defined as: s ( µ B , T ) = ∂p ( µ B , T ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ . (20) For dimensionless quantities:ˆ s (ˆ µ B , T ) = 4 ˆ p (ˆ µ B , T ) + T ∂ ˆ p (ˆ µ B , T ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ = (21)= 4 ˆ p (ˆ µ B , T ) + T ∂ ˆ p (ˆ µ B , T ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ˆ µ − ˆ µ B ˆ n B (ˆ µ B , T )where in the last step we converted the derivative at con-stant µ B into a derivative at constant ˆ µ B .The T -derivative of the pressure involves the chain rule T ∂ ˆ p (ˆ µ B , T ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ˆ µ = T ∂ ˆ p (0 , T ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ˆ µ (22)+ 12 (cid:90) ˆ µ B T dχ B ( T (cid:48) ) dT (cid:48) (cid:12)(cid:12)(cid:12)(cid:12) T (cid:48) = T ( κ BB y + κ BB y ) ×× (cid:20) κ BB y + κ BB y + T (cid:18) dκ BB dT y + dκ BB dT y (cid:19)(cid:21) dy where dχ B ( T ) dT is calculated at µ B = 0 as already de-scribed.The dimensionless expression for the energy densityˆ (cid:15) = (cid:15)/T follows as:ˆ (cid:15) (ˆ µ B , T ) = ˆ s (ˆ µ B , T ) − ˆ p (ˆ µ B , T ) + ˆ µ B ˆ n B (ˆ µ B , T ) . (23)Finally, we present our results for the finite real chem-ical potential extrapolation of several thermodynamicquantities. The various panels of Fig. 9 show the baryondensity, pressure, entropy, energy density, strangenessdensity and χ S for ˆ µ B = 0 − .
5. Alongside our results, weshow predictions from the HRG model for
T <
150 MeV,which we find in very good agreement with our extrap-olation for all observables, at all values of the chemicalpotential.We also note that in all cases, the observables do notsuffer from pathological behavior. The uncertainties areunder control for our range of chemical potentials, whichhighly improves on the results currently achievable viaTaylor expansion.We devote the two panels of Fig. 10 to the compari-son of our results for the baryon density (left) and en-ergy density (right) to the simplified case where κ BB isneglected. We can appreciate how the inclusion of thenext-to-leading-order parameter came at the cost of anincreased uncertainty at larger chemical potential. Thisdoes not come unexpected, as we saw from our resultsthat κ BB was compatible with zero at all temperatures.In the case of the energy density, which is dominated bythe µ B = 0 contribution, hardly any effect is visible. VI. CONCLUSIONS
In this work, we proposed an alternative summationscheme for the equation of state of QCD at finite real0 n B / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 3 µ B /T = 2.5 µ B /T = 2 µ B /T = 1.5 µ B /T = 1 µ B /T = 0.5 µ B /T = 0 p / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 3 µ B /T = 2.5 µ B /T = 2 µ B /T = 1.5 µ B /T = 1 µ B /T = 0.5 µ B /T = 0 s / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 3 µ B /T = 2.5 µ B /T = 2 µ B /T = 1.5 µ B /T = 1 µ B /T = 0.5 µ B /T = 0 ε / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 3 µ B /T = 2.5 µ B /T = 2 µ B /T = 1.5 µ B /T = 1 µ B /T = 0.5 µ B /T = 0 - n s / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 3 µ B /T = 2.5 µ B /T = 2 µ B /T = 1.5 µ B /T = 1 µ B /T = 0.5 µ B /T = 0 χ S ( T ) T [MeV] µ B /T = 3.5 µ B /T = 3 µ B /T = 2.5 µ B /T = 2 µ B /T = 1.5 µ B /T = 1 µ B /T = 0.5 µ B /T = 0 FIG. 9. Baryon density, pressure, entropy, energy density, strangeness density and χ S at increasing values of ˆ µ B . With solidlines we show the results from the HRG model. chemical potential, designed to overcome the shortcom-ings which are characteristic of the Taylor expansion ap-proach. Combining simulations at both zero and imag-inary chemical potentials, we determined the LO andNLO parameters describing the chemical potential de-pendence of the baryon density, which we could then ex-trapolate to large real chemical potentials.By combining this new element, and previously pub- lished results for the EoS at vanishing density, we couldreconstruct all thermodynamic variables at chemical po-tentials as large as µ B /T = 3 . κ ij and κ ij , suggest that theavenue we pursue in this work is rather promising for1 n B / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 2.5 µ B /T = 1.5 µ B /T = 0.5 ε / T ( T ) T [MeV] µ B /T = 3.5 µ B /T = 2.5 µ B /T = 1.5 µ B /T = 0.5 µ B /T = 0 FIG. 10. Comparison of baryon density (left) and energy density (right) at different values of ˆ µ B in the case where a κ BB parameter is used (darker shades) or omitted (lighter shades). The HRG results are shown with solid lines. the description of QCD thermodynamics at finite chemi-cal potential. Moreover, our procedure is systematicallyimprovable with sufficient computing power, and mightprove to be a better strategy than existing “canonical”approaches.In this work we limited ourselves to the case where thestrange and electric chemical potentials are set to zero.We reserve for future work the exploration of the phe-nomenologically relevant case of strangeness neutralityand fixed electric charge to baryon ratio. Acknowledgments [1] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Sz-abo, The Order of the quantum chromodynamics transi-tion predicted by the standard model of particle physics,Nature , 675 (2006), arXiv:hep-lat/0611014 [hep-lat].[2] J. N. Guenther, Overview of the QCD phase dia-gram - Recent progress from the lattice, (2020),arXiv:2010.15503 [hep-lat].[3] C. Ratti, Lattice QCD and heavy ion collisions: a reviewof recent progress, Rept. Prog. Phys. , 084301 (2018),arXiv:1804.07810 [hep-lat].[4] M. Cristoforetti, F. Di Renzo, and L. Scorzato (Aurora-Science), New approach to the sign problem in quantumfield theories: High density QCD on a Lefschetz thimble,Phys. Rev. D , 074506 (2012), arXiv:1205.3996 [hep- lat].[5] A. Alexandru, G. Basar, and P. Bedaque, Monte Carloalgorithm for simulating fermions on Lefschetz thimbles,Phys. Rev. D , 014504 (2016), arXiv:1510.03258 [hep-lat].[6] M. Fukuma, N. Matsumoto, and N. Umeda, Implemen-tation of the HMC algorithm on the tempered Lefschetzthimble method, (2019), arXiv:1912.13303 [hep-lat].[7] G. Aarts, F. A. James, J. M. Pawlowski, E. Seiler,D. Sexty, and I.-O. Stamatescu, Stability of complexLangevin dynamics in effective models, JHEP , 073,arXiv:1212.5231 [hep-lat].[8] D. Sexty, Simulating full QCD at nonzero density usingthe complex Langevin equation, Phys. Lett. B , 108 (2014), arXiv:1307.7748 [hep-lat].[9] M. Scherzer, E. Seiler, D. Sexty, and I.-O. Stamatescu,Controlling Complex Langevin simulations of latticemodels by boundary term analysis, Phys. Rev. D ,014501 (2020), arXiv:1910.09427 [hep-lat].[10] M. Giordano, K. Kapas, S. D. Katz, D. Nogradi, andA. Pasztor, New approach to lattice QCD at finite den-sity; results for the critical end point on coarse lattices,JHEP , 088, arXiv:2004.10800 [hep-lat].[11] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek,F. Karsch, E. Laermann, C. Schmidt, and L. Scorzato,The QCD thermal phase transition in the presence of asmall chemical potential, Phys. Rev. D66 , 074507 (2002),arXiv:hep-lat/0204010 [hep-lat].[12] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kacz-marek, F. Karsch, E. Laermann, and K. Redlich, Ther-modynamics of two flavor QCD to sixth order in quarkchemical potential, Phys. Rev.
D71 , 054508 (2005),arXiv:hep-lat/0501030 [hep-lat].[13] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg,C. Ratti, and K. K. Szabo, QCD equation of stateat nonzero chemical potential: continuum results withphysical quark masses at order mu , JHEP , 053,arXiv:1204.6710 [hep-lat].[14] A. Bazavov et al. , The QCD Equation of State to O ( µ B )from Lattice QCD, Phys. Rev. D95 , 054504 (2017),arXiv:1701.04325 [hep-lat].[15] O. Kaczmarek, F. Karsch, E. Laermann, C. Miao,S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner,and W. Unger, Phase boundary for the chiral transitionin (2+1) -flavor QCD at small values of the chemical po-tential, Phys. Rev.
D83 , 014504 (2011), arXiv:1011.3130[hep-lat].[16] G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, TheQCD phase diagram at nonzero quark density, JHEP ,001, arXiv:1102.1356 [hep-lat].[17] C. Bonati, M. D’Elia, F. Negro, F. Sanfilippo, andK. Zambello, Curvature of the pseudocritical line inQCD: Taylor expansion matches analytic continuation,Phys. Rev. D98 , 054510 (2018), arXiv:1805.02960 [hep-lat].[18] A. Bazavov et al. (HotQCD), Chiral crossover in QCD atzero and non-zero chemical potentials, Phys. Lett.
B795 ,15 (2019), arXiv:1812.08235 [hep-lat].[19] M. D’Elia, G. Gagliardi, and F. Sanfilippo, Higher orderquark number fluctuations via imaginary chemical po-tentials in N f = 2 + 1 QCD, Phys. Rev. D , 094503(2017), arXiv:1611.08285 [hep-lat].[20] S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K.Szabo, A. Pasztor, I. Portillo, and C. Ratti, Higher orderfluctuations and correlations of conserved charges fromlattice QCD, JHEP , 205, arXiv:1805.04445 [hep-lat].[21] P. de Forcrand and O. Philipsen, The QCD phasediagram for small densities from imaginary chemicalpotential, Nucl. Phys. B642 , 290 (2002), arXiv:hep-lat/0205016 [hep-lat].[22] M. D’Elia and M.-P. Lombardo, Finite density QCD viaimaginary chemical potential, Phys. Rev.
D67 , 014505(2003), arXiv:hep-lat/0209146 [hep-lat].[23] M. D’Elia, F. Di Renzo, and M. P. Lombardo, TheStrongly interacting quark gluon plasma, and the criticalbehaviour of QCD at imaginary mu, Phys. Rev.
D76 ,114509 (2007), arXiv:0705.3814 [hep-lat]. [24] P. Cea, L. Cosmai, M. D’Elia, C. Manneschi, andA. Papa, Analytic continuation of the critical line: Sug-gestions for QCD, Phys. Rev. D , 034501 (2009),arXiv:0905.1292 [hep-lat].[25] P. Cea, L. Cosmai, M. D’Elia, C. Manneschi, andA. Papa, On the analytic continuation of the critical line,PoS LAT2009 , 161 (2009), arXiv:1001.4439 [hep-lat].[26] C. Bonati, M. D’Elia, P. de Forcrand, O. Philipsen,and F. Sanfilippo, The chiral phase transition for two-flavour QCD at imaginary and zero chemical potential,PoS
LATTICE2013 , 219 (2014), arXiv:1311.0473 [hep-lat].[27] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro,and F. Sanfilippo, Curvature of the chiral pseudocriticalline in QCD: Continuum extrapolated results, Phys. Rev.
D92 , 054503 (2015), arXiv:1507.03571 [hep-lat].[28] P. Cea, L. Cosmai, and A. Papa, Critical line of 2+1 fla-vor QCD: Toward the continuum limit, Phys. Rev.
D93 ,014507 (2016), arXiv:1508.07599 [hep-lat].[29] R. Bellwied, S. Borsanyi, Z. Fodor, J. Guenther, S. D.Katz, C. Ratti, and K. K. Szabo, The QCD phase dia-gram from analytic continuation, Phys. Lett.
B751 , 559(2015), arXiv:1507.07510 [hep-lat].[30] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D.Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo,QCD Crossover at Finite Chemical Potential from Lat-tice Simulations, Phys. Rev. Lett. , 052001 (2020),arXiv:2002.02821 [hep-lat].[31] F. Karsch, B.-J. Schaefer, M. Wagner, and J. Wambach,Towards finite density QCD with Taylor expansions,Phys. Lett. B , 256 (2011), arXiv:1009.5211 [hep-ph].[32] P. Cea, L. Cosmai, M. D’Elia, A. Papa, and F. Sanfil-ippo, The critical line of two-flavor QCD at finite isospinor baryon densities from imaginary chemical potentials,Phys. Rev. D , 094512 (2012), arXiv:1202.5700 [hep-lat].[33] S. Datta, R. V. Gavai, and S. Gupta, Quark number sus-ceptibilities and equation of state at finite chemical po-tential in staggered QCD with Nt=8, Phys. Rev. D95 ,054512 (2017), arXiv:1612.06673 [hep-lat].[34] A. P´asztor, Z. Sz´ep, and G. Mark´o, Apparent conver-gence of Pad \ ’e approximants for the crossover line infinite density QCD, (2020), arXiv:2010.00394 [hep-lat].[35] C. Schmidt, J. Goswami, G. Nicotra, F. Ziesch´e, P. Di-mopoulos, F. Di Renzo, S. Singh, and K. Zambello,Net-baryon number fluctuations, in Criticality in QCDand the Hadron Resonance Gas (2021) arXiv:2101.02254[hep-lat].[36] S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, andA. Vuorinen, Equation of State of hot and dense QCD:Resummed perturbation theory confronts lattice data,JHEP , 055, arXiv:1307.8098 [hep-ph].[37] N. Haque, J. O. Andersen, M. G. Mustafa, M. Strickland,and N. Su, Three-loop HTLpt Pressure and Susceptibili-ties at Finite Temperature and Density, Phys.Rev. D89 ,061701 (2014), arXiv:1309.3968 [hep-ph].[38] N. Haque and M. Strickland, NNLO HTLpt predictionsfor the curvature of the QCD phase transition line,(2020), arXiv:2011.06938 [hep-ph].[39] R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasz-tor, C. Ratti, and K. K. Szabo, Fluctuations and correla-tions in high temperature QCD, Phys. Rev.
D92 , 114505(2015), arXiv:1507.04627 [hep-lat]. [40] H. T. Ding, S. Mukherjee, H. Ohno, P. Petreczky, andH. P. Schadler, Diagonal and off-diagonal quark numbersusceptibilities at high temperatures, Phys. Rev. D92 ,074043 (2015), arXiv:1507.06637 [hep-lat].[41] S. Pratt, E. Sangaline, P. Sorensen, and H. Wang, Con-straining the Eq. of State of Super-Hadronic Matterfrom Heavy-Ion Collisions, Phys. Rev. Lett. , 202301(2015), arXiv:1501.04042 [nucl-th].[42] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D.Katz, S. Krieg, C. Ratti, and K. K. Szabo, The QCDequation of state with dynamical quarks, JHEP , 077,arXiv:1007.2580 [hep-lat].[43] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg,and K. K. Szabo, Full result for the QCD equation ofstate with 2+1 flavors, Phys. Lett. B730 , 99 (2014),arXiv:1309.5258 [hep-lat].[44] A. Bazavov et al. (HotQCD), Equation of state in (2+1 )-flavor QCD, Phys. Rev.
D90 , 094503 (2014),arXiv:1407.6387 [hep-lat].[45] J. Gunther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D.Katz, A. Pasztor, and C. Ratti, The QCD equationof state at finite density from analytical continuation,
Proceedings, 12th Conference on Quark Confinementand the Hadron Spectrum (Confinement XII): Thes-saloniki, Greece , EPJ Web Conf. , 07008 (2017),arXiv:1607.02493 [hep-lat].[46] A. Bazavov et al. , Skewness, kurtosis, and the fifth andsixth order cumulants of net baryon-number distributionsfrom lattice QCD confront high-statistics STAR data,Phys. Rev. D , 074502 (2020), arXiv:2001.08530 [hep-lat].[47] W.-j. Fu, X. Luo, J. M. Pawlowski, F. Rennecke, R. Wen,and S. Yin, Hyper-order baryon number fluctuations atfinite temperature and density, (2021), arXiv:2101.06035[hep-ph].[48] R. Bellwied, S. Borsanyi, Z. Fodor, J. N. Guenther,J. Noronha-Hostler, P. Parotto, A. Pasztor, C. Ratti,and J. M. Stafford, Off-diagonal correlators of con-served charges from lattice QCD and how to relatethem to experiment, Phys. Rev. D , 034506 (2020), arXiv:1910.14592 [hep-lat].[49] J. Noronha-Hostler, P. Parotto, C. Ratti, and J. M.Stafford, Lattice-based equation of state at finite baryonnumber, electric charge and strangeness chemical poten-tials, Phys. Rev. C , 064910 (2019), arXiv:1902.06723[hep-ph].[50] A. Monnai, B. Schenke, and C. Shen, Equation of stateat finite densities for QCD matter in nuclear collisions,Phys. Rev. C , 024907 (2019), arXiv:1902.05095[nucl-th].[51] B. Friman, F. Karsch, K. Redlich, and V. Skokov, Fluctu-ations as probe of the QCD phase transition and freeze-out in heavy ion collisions at LHC and RHIC, Eur.Phys.J.
C71 , 1694 (2011), arXiv:1103.3511 [hep-ph].[52] M. Giordano and A. P´asztor, Reliable estimation of theradius of convergence in finite density QCD, Phys. Rev.D , 114510 (2019), arXiv:1904.01974 [hep-lat].[53] M. Giordano, K. Kapas, S. D. Katz, D. Nogradi, andA. Pasztor, Radius of convergence in lattice QCD at finite µ B with rooted staggered fermions, Phys. Rev. D ,074511 (2020), arXiv:1911.00043 [hep-lat].[54] A. Connelly, G. Johnson, S. Mukherjee, and V. Skokov,Universality driven analytic structure of QCD crossover:radius of convergence and QCD critical point, Nucl. Phys.A , 121834 (2021), arXiv:2004.05095 [hep-ph].[55] C. Ratti, S. Roessner, and W. Weise, Quark number sus-ceptibilities: Lattice QCD versus PNJL model, Phys.Lett. B , 57 (2007), arXiv:hep-ph/0701091.[56] R. Critelli, J. Noronha, J. Noronha-Hostler, I. Por-tillo, C. Ratti, and R. Rougemont, Critical point in thephase diagram of primordial quark-gluon matter fromblack hole physics, Phys. Rev. D , 096026 (2017),arXiv:1706.00455 [nucl-th].[57] P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang,J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Sch¨afer,and M. Stephanov, QCD equation of state matched tolattice data and exhibiting a critical point singularity,Phys. Rev. C101