# 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 ﬁnite 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 suﬀers from shortcomings at chemical potentials µ B ≥ (2 − . T . First, one faces diﬃculties inherent in performing such an expansion with alimited number of coeﬃcients; second, higher order coeﬃcients determined from lattice calculationssuﬀer from a poor signal-to-noise ratio. In this work, we present a novel scheme for extrapolatingthe equation of state of QCD to ﬁnite, real chemical potential that can extend its reach furtherthan previous methods. We present continuum extrapolated lattice results for the new expansioncoeﬃcients and show the thermodynamic observables up to µ B /T ≤ . I. INTRODUCTION

The phase diagram of Quantum Chromodynamics(QCD) is an open ﬁeld of investigation, which is cur-rently at the center of intense eﬀorts from the theoret-ical and experimental communities alike. At vanishingbaryon density, the transition between conﬁned and de-conﬁned 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 ﬁnite 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 ﬁnite 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 ﬁnite 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 coeﬃcients can be eﬃciently 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 ﬁnite density physics, which is often referredto as analytical continuation. This name suggests that,besides the computation of the Taylor expansion coeﬃ-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-eﬃcient [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 coeﬃcients in Ref. [20], at µ B /T = 1 , ,

3, as a function of thetemperature. Diﬀerent 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 ﬁnitedensity 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 ﬁrst continuum extrapo-lated extension to ﬁnite µ 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 eﬀective theory with functional methods[47].In this work, we propose a new extrapolation schemeto ﬁnite 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 ﬁnite 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 coeﬃcients 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 coeﬃcients, one can also deﬁne oﬀ-diagonal correlators between diﬀerent conserved chargesin QCD. Correlators between baryon number andstrangeness, which we will need in our procedure, aredeﬁned 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 ﬁnite T, µ B , µ S , µ Q [49, 50]. We will use theˆ µ i = µ i /T shorthand notation in this manuscript.Currently, results for the expansion coeﬃcients 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 coeﬃcients.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 ﬁnite chemical potential are dictated by the µ B = 0 temperature dependence of the last coeﬃcient in-cluded in the expansion. Hence, the structures appearingaround the QCD transition temperature in higher ordercoeﬃcients are “translated” into the ﬁnite- µ 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 coeﬃcients at µ B = 0and a certain temperature T , determine the equation ofstate at the same T at ﬁnite µ B , while the pseudo-criticaltemperature T pc might have varied considerably. While asuﬃciently large number of expansion coeﬃcients wouldlead to smooth extrapolated functions, even though theTaylor coeﬃcients themselves present a complex struc-ture around the transition temperature, the problem hereis rather practical. A scheme that could work with fewercoeﬃcients 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 coeﬃcients inRef. [20], at ˆ µ B = 1 , ,

3. The extrapolation is shownincluding an increasing number of coeﬃcients, to showthe eﬀect 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 thecoeﬃcients in Ref. [20] causes unphysical non-monotonicbehavior. Ultimately, a pathological behaviour could alsocome from a ﬁnite 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 ﬁxedimaginary µ B /T ratios. The 0 / µ B = 0 can beeasily resolved and equals χ B ( T ).The T -dependence of the normalized baryon den-sity at ﬁnite 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 ﬂat. A simple rescaling oftemperatures can be described as: χ B ( T, ˆ µ B )ˆ µ B = χ B ( T (cid:48) , , (4)where the actual temperature diﬀerence 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 ﬁnite-ˆ µ 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 ﬁrst and second order ﬂuctuationsof 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 deﬁned analogously to Eq. (5), albeit withdiﬀerent κ parameters.In order to motivate our alternative summationscheme, let us ﬁrst consider a crude approximation thatwe will later reﬁne. We take Eq. (5) at face valueand use it together with Eq. (4) to obtain a well de-ﬁned χ B ( T, ˆ µ B ) function. We need a χ B ( T,

0) functionas well, which we borrow from a deliberately simple ﬁt 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 issuﬃcient 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 eﬀectsin 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 suﬀer from the additional complications of signal ex-traction for higher order expansion coeﬃcients, 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 ﬁnitechemical 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 ﬁnite 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 ﬁnite µ B . Analogousparameters will be introduced below to for the case of χ S / ˆ µ B ( κ BS and κ BS ) and of χ S ( κ SS and κ SS ) at ﬁnite µ B . In a way, this formalism replaces the ﬁxed temper-ature µ B expansion by a ﬁxed-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 simpliﬁed 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 ﬂuctuations 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 diﬀerence of two competing terms. Notice, too, thatEq. (10) contains temperature-derivatives of the χ ( T ) co-eﬃcients, which may pose a numerical challenge, unlessthe coeﬃcients are known at suﬃcient 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 coeﬃcients, κ 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 κ coeﬃ-cients are closely related to ours.Contrary to Ref. [14], we use Eq. (8) as the deﬁnitionof 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 coeﬃcients in Eq. (8), which then canbe used to either calculate the Taylor coeﬃcients or, evenbetter, to extrapolate the equation of state at ﬁnite µ B with no reference to the Taylor coeﬃcients themselves.Also, we used imaginary chemical potentials notonly to calculate the coeﬃcients, but also to studythe single observable ﬁrst, 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 beneﬁts from tree-level Symanzik improvement in the gauge sector and fourlevels of stout smearing for the staggered ﬂavors. 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 coeﬃcients κ ij and κ ij For the determination of κ ij and κ ij , one can takeadvantage of simulations both at zero and ﬁnite 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 deﬁning 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 ﬁrst 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 ﬁtin ˆ µ B and obtain the expansion coeﬃcients. This canbe done separately for each lattice. However, we pre-fer to combine the ˆ µ B and continuum ﬁts in one two-dimensional ﬁtting procedure. This combined ﬁt 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 ﬁt or not,or consider the ﬁt of 1 / Π instead. Also the range inimaginary µ B is arbitrary: we consider Im µ B ≤ . µ B ≤ .

4. When diﬀerent lattice spacings are ﬁttedtogether 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 ﬁnite 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 ﬁt 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 ﬁts to perform a continuum extrapolation of κ ij ( T )and κ ij ( T ). After dropping the ﬁts 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 deﬁnes 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 ﬁt procedure for the param-eters κ BB ( T ) and κ BB ( T ), alongside the correspondingexpectations from the Hadron Resonance Gas (HRG)model. We ﬁnd 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 coeﬃcients.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 signiﬁcant asit seems at a ﬁrst 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 ﬁnite real chemical poten-tial, we construct a smoother version of our ﬁnal resultsfor κ ij and κ ij , in order to limit the inﬂuence of numer-ical eﬀects on the ﬁnal observables, but also to facilitatethe temperature derivatives that enter the entropy for-mula. Thanks to the very mild dependence of the co-eﬃcients on the temperature, we perform a polynomialﬁt of order 5 for the κ ij , and of order 2 for the evenless T -dependent κ ij . The very good ﬁt 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 ﬁttwo 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 ﬁt too strongly. For this reason, thesearbitrary errors are chosen to be smaller, but compara-ble to the lattice ones. We note that the ﬁts we performtake fully into account the correlations between results atdiﬀerent temperatures, systematic as well as statistical.Thus we encode all errors into the (correlated) errors ofthe coeﬃcients of a polynomial.In Fig. 7 we show the results of the ﬁts 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 ﬁts 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 ﬁtted 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 ﬁt 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 ﬁtted data in lighter shades. The HRG pointsincluded in the ﬁt are shown as well. C. Continuum result on χ B ( T ) and its temperaturederivative In order to determine the value of thermodynamicquantities at ﬁnite 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 ﬁt the splinesin temperature and in 1 /N τ in one step. The ﬁtted 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 ﬁtqualities (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 ﬁt 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 ﬁtting 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 ﬁnal 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 deﬁned 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 ﬁnite 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 ﬁnd 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 notsuﬀer 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 simpliﬁed 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 eﬀect is visible. VI. CONCLUSIONS

In this work, we proposed an alternative summationscheme for the equation of state of QCD at ﬁnite 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 diﬀerent 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 ﬁnite chemi-cal potential. Moreover, our procedure is systematicallyimprovable with suﬃcient 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 ﬁxed 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 quantumﬁeld 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 eﬀective 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 ﬁnite 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 ﬂavor 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) -ﬂavor 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. Sanﬁlippo, 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. Sanﬁlippo, Higher orderquark number ﬂuctuations 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 orderﬂuctuations 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. Sanﬁlippo, The chiral phase transition for two-ﬂavour 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. Sanﬁlippo, 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 ﬂa-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 ﬁnite 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. Sanﬁl-ippo, The critical line of two-ﬂavor QCD at ﬁnite 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 ﬁnite 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 inﬁnite 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 ﬂuctuations, 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 oﬀ-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 ﬂavors, Phys. Lett. B730 , 99 (2014),arXiv:1309.5258 [hep-lat].[44] A. Bazavov et al. (HotQCD), Equation of state in (2+1 )-ﬂavor 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 ﬁnite density from analytical continuation,

Proceedings, 12th Conference on Quark Conﬁnementand the Hadron Spectrum (Conﬁnement XII): Thes-saloniki, Greece , EPJ Web Conf. , 07008 (2017),arXiv:1607.02493 [hep-lat].[46] A. Bazavov et al. , Skewness, kurtosis, and the ﬁfth 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 ﬂuctuations atﬁnite 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. Staﬀord, Oﬀ-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.Staﬀord, Lattice-based equation of state at ﬁnite 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 ﬁnite 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 ﬁnite 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 ﬁnite µ 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