# Corrections to the hadron resonance gas from lattice QCD and their effect on fluctuation-ratios at finite density

Rene Bellwied, Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Sandor D. Katz, Paolo Parotto, Attila Pasztor, David Pesznyak, Claudia Ratti, Kalman K. Szabo

CCorrections to the hadron resonance gas from lattice QCDand their eﬀect on ﬂuctuation-ratios at ﬁnite density

Rene Bellwied and Claudia Ratti

Department of Physics, University of Houston, Houston, TX 77204, USA

Szabolcs Bors´anyi and Paolo Parotto

Department of Physics, Wuppertal University, Gaussstr. 20, D-42119, Wuppertal, Germany

Zolt´an Fodor

Department of Physics, Wuppertal University, Gaussstr. 20, D-42119, Wuppertal, GermanyPennsylvania State University, Department of Physics, State College, PA 16801, USAInst. for Theoretical Physics, ELTE E¨otv¨os Lor´and University,P´azm´any P. s´et´any 1/A, H-1117 Budapest, Hungary andJ¨ulich Supercomputing Centre, Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany

Jana N. Guenther

Aix Marseille Univ, Universit´e de Toulon, CNRS, CPT, Marseille, France

S´andor D. Katz, Attila P´asztor, ∗ and D´avid Peszny´ak Inst. for Theoretical Physics, ELTE E¨otv¨os Lor´and University,P´azm´any P. s´et´any 1/A, H-1117 Budapest, Hungary

K´alm´an K. Szab´o

Department of Physics, Wuppertal University, Gaussstr. 20, D-42119, Wuppertal, Germany andJ¨ulich Supercomputing Centre, Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany (Dated: February 15, 2021)The hadron resonance gas (HRG) model is often believed to correctly describe the conﬁned phaseof QCD. This assumption is the basis of many phenomenological works on QCD thermodynamicsand of the analysis of hadron yields in relativistic heavy ion collisions. We use ﬁrst-principle latticesimulations to calculate corrections to the ideal HRG. Namely, we determine the sub-leading fugacityexpansion coeﬃcients of the grand canonical free energy, receiving contributions from processes likekaon-kaon or baryon-baryon scattering. We achieve this goal by performing a two dimensional scanon the imaginary baryon number chemical potential ( µ B ) - strangeness chemical potential ( µ S ) plane,where the fugacity expansion coeﬃcients become Fourier coeﬃcients. We carry out a continuumlimit estimation of these coeﬃcients by performing lattice simulations with temporal extents of N τ = 8 , ,

12 using the 4stout-improved staggered action. We then use the truncated fugacityexpansion to extrapolate ratios of baryon number and strangeness ﬂuctuations and correlations toﬁnite chemical potentials. Evaluating the fugacity expansion along the crossover line, we reproducethe trend seen in the experimental data on net-proton ﬂuctuations by the STAR collaboration.

I. INTRODUCTION

The study of the QCD phase diagram has been a veryactive area of research for the last few decades. Whilemuch is known about the thermodynamics of QCD atzero baryon number chemical potential, such as the tem-perature of the crossover transition [1–4] and the equa-tion of state [5–8], the properties of the theory at ﬁnitebaryon densities remain elusive. Eﬀective models pre-dict that the crossover transition turns into a real phasetransition at a critical endpoint [9–11]. However, conﬁr-mation of this feature is needed from a ﬁrst principlesapproach and/or experiment. The main goal of the cur-rently ongoing experimental eﬀort at the second Beam ∗ Corresponding author: [email protected]

Energy Scan program at RHIC in 2019-2021 is locatingthe supposed critical endpoint of QCD.Direct ﬁrst principle lattice simulations at ﬁnite chem-ical potential are hampered by the infamous sign prob-lem [12]. Methods to circumvent it include reweight-ing [13–18], Taylor expansion around zero chemical po-tential [19–30], and extrapolation from purely imaginarychemical potential [31–46]. The ﬁrst of these methodshas so far proved too expensive to apply on ﬁne lattices.Therefore, no continuum extrapolated results exist withthis approach so far. The latter two methods, on theother hand, involve analytic continuation, which is an ill-posed problem, regardless of whether the available data isa number of Taylor coeﬃcients at zero chemical potentialor the value of some observable at a number of points atimaginary chemical potential. In such a case, it is impor-tant to use physical insight to argue what the functional a r X i v : . [ h e p - l a t ] F e b form of a given observable could be as a function of thechemical potential.The conﬁned phase of QCD is often assumed to bewell described by the ideal hadron resonance gas (HRG)model [47–51]. The HRG model is based on the assump-tion that a gas of interacting hadrons can be describedas a gas of non-interacting hadrons and resonances. Theinclusion of the resonances as free particles is an approxi-mate way of taking into account resonant interactions be-tween the stable hadrons [49, 50]. The model describesbulk thermodynamic observables - like the pressure orthe energy density - obtained from ﬁrst principle latticecalculations rather well at zero chemical potential [6, 52–54]. However, when looking at observables probing ﬁnitechemical potentials, namely Taylor expansion coeﬃcientsin the chemical potentials µ B , µ S and µ Q near zero orin the respective fugacities e µ B /T , e µ S /T and e µ Q /T near1, some discrepancies start to emerge between the HRGmodel and lattice calculations. Some of these discrepan-cies can be traced back to the fact that some fugacityexpansion coeﬃcients with baryon number zero and oneare underestimated by the HRG model. Though this isnot the only possibility, this type of discrepancy can beinterpreted within the bounds of the HRG model itself,and has been used to try to infer the existence of as ofyet unobserved hadrons [41, 55–57]. The other possibil-ity is that instead of more resonances, a better treat-ment of resonances is needed, taking into account ﬁnitewidths and also non-resonant interactions [48–50, 58–60].Of course, both statements can be true at the same time.For a precise description of the thermodynamics, we mostlikely need better knowledge of the mass spectrum athigher energies, as well as a more accurate treatment ofresonances.Other discrepancies between the ideal HRG and thelattice are impossible to resolve by supposing the exis-tence of more resonances. These discrepancies can betraced back to the observation that even in the tempera-ture range below the crossover, the HRG fails to describesub-leading fugacity expansion coeﬃcients.In principle, the hadron resonance gas can be sys-tematically improved by the S-matrix formulation dueto Dashen, Bernstein and Ma [48–50], which allows forthe calculation of the fugacity expansion coeﬃcients, ifenough information is known about the scattering ma-trix of the hadrons. Applying the S-matrix formalismcan lead to a better description of QCD thermodynam-ics. Ref. [60] shows that the baryon-electric charge corre-lation χ BQ is particularly sensitive to pion-nucleon scat-tering phase shifts and that the inclusion of these phaseshifts into a hadron gas analysis leads to an improveddescription of the lattice data. This observable is some-what special though, in that if isospin symmetry is as-sumed | S | = 1 hyperons do not contribute at all to χ BQ ,and therefore it is only pion-nucleon scattering that dom-inates the non-resonant contributions.For other observables more scattering data, e.g. in-formation about baryon-baryon scattering would also be needed. This is especially the case at ﬁnite baryon den-sity. Unfortunately, information on these scattering pro-cesses is only partially available. While the nucleon-nucleon elastic scattering phase shifts are known exper-imentally [61–63], the inelastic part of the S-matrix isnot known. Even less is known about scattering be-tween hadrons other than nucleons. Hyperon-nucleonand hyperon-hyperon interactions have been studied inchiral eﬀective theory [64, 65]. In the last few years, theanalysis of momentum space correlations for hadron pairsmeasured in pp and p-Pb collisions has also been usedto infer properties of hadron-hadron interactions [66–68].There are also some lattice results available for baryon-baryon scattering [69–71], but not yet with a continuumextrapolation. While these mentioned research directionsshow a clear eﬀort from the community to learn aboutscattering between other hadrons, the preliminary natureof these results makes the use of the S-matrix formalismfor the fugacity expansion impractical at the moment.One simple way to nevertheless go beyond the idealhadron resonance gas is to use some kind of mean ﬁeldmodel for the short range repulsion and the long rangeattraction between the baryons. Such models were com-pared to lattice results in Refs [42, 72–74]. These worksin particular emphasized the importance of the hardcore repulsive interactions between hadrons when de-scribing thermodynamics at ﬁnite baryon chemical po-tential. This type of interaction is completely absent inthe ideal HRG and leads to a sizable negative contri-bution to the fugacity expansion coeﬃcients with baryonnumber two. Such approaches, while interesting, are veryfar from the ﬁrst principle approach that the S-matrixformulation could provide, were the necessary S-matrixelements known. In fact, the ﬂavour dependence of ex-cluded volume parameters used in the literature so farhave been quite arbitrary, often assuming the same ex-cluded volume for all hadrons. We believe the presentcalculation of the fugacity expansion coeﬃcients can leadto the construction of more realistic models.Going beyond equilibrium in the grand canonical en-semble, versions of the hadron resonance gas model havealso been used to interpret hadron yields in heavy ioncollision experiments. This approach is colloquially re-ferred to as thermal ﬁts as they involve the estimation ofthe temperature and chemical potential where the yieldsof hadrons are frozen, the so-called chemical freeze-outconditions. This approach was successful in describinghadron yields [75–80], which is quite remarkable, consid-ering that these yields at a single collision energy spanmany orders of magnitude. Though an important under-lying assumption here is the equilibration of the systemproduced in heavy ion collisions [81–83], the fact that theﬁts work also provides some evidence for this assump-tion. In this context, it has been realized that includingthe pion-proton phase shifts in the analysis changes thepredicted yields as compared to the ideal HRG at LHCenergies [84]. Of course, for consistency, one should ex-tend such and S-matrix treatment to strange hadrons aswell. In lack of the necessary scattering data, this ex-tension to strangeness is not straightforward [85]. Theinclusion of these corrections is important to preciselytest the assumption of a single freeze-out temperature.As a competitor to this assumption, in the context ofideal HRG, it was shown [86] that diﬀerent freeze-outtemperatures for light and strange hadrons, can signiﬁ-cantly improve the description of the experimental yieldsat LHC and the highest RHIC energies.Since comparisons with the available lattice data sug-gest that the agreement between full QCD and the idealhadron resonance gas model gets worse at ﬁnite chemicalpotential, we suspect that non-resonant scattering eﬀectswill be even more important at the RHIC Beam EnergyScan and future experiments at lower collision energies,like FAIR and NICA.In this work we calculate sub-leading fugacity expan-sion coeﬃcients with ﬁrst principle lattice simulations.To this end we perform simulations at imaginary chem-ical potentials, where the fugacity expansion coeﬃcientsturn into Fourier coeﬃcients in the imaginary values ofthe chemical potentials. This correspondence was al-ready exploited in our earlier works. In Ref. [41] wemade a detailed analysis of the fugacity expansion co-eﬃcients already appearing in the Boltzmann approxi-mation of the ideal HRG, to infer the existence of notyet discovered strange hadrons. In Ref. [42], some of usused the fugacity expansion to emphasize the importanceof repulsive baryonic interactions near the crossover re-gion. In Ref. [87] we compared the fugacity expansionwith the Taylor expansion in the chemical potentials forcross-correlators of conserved charges. Here we go be-yond our earlier works by performing lattice simulationson a 2 dimensional grid in the purely imaginary ( µ IB , µ IS )plane. This allows us for the ﬁrst time to separate thescattering contributions to QCD thermodynamics by thenet strangeness quantum number of the participants. Inaddition to giving insight on the origin of the discrep-ancies between full QCD and the ideal HRG model, webelieve our results on the fugacity expansion coeﬃcientswill also be useful to tune the parameters of the freeze-out models in heavy ion phenomenology.We also use the truncated fugacity expansion to ex-trapolate experimentally measured ratios of baryon num-ber and strangeness susceptibilities to ﬁnite baryonchemical potentials on the phenomenologically relevantstrangeness neutral line. This provides an alternativeextrapolation procedure to the standard Taylor method.When we extrapolate on the crossover line at strangenessneutrality, these sub-leading coeﬃcients approximatelyreproduces the trend seen in the experimental data ofthe STAR collaboration of net-proton ﬂuctuations [88–90].The structure of the paper is as follows. In the nextsection, we introduce the basic notation and observablesused in our study. In Sec. III we discuss our lattice setup.In Sec. IV we discuss our ﬁtting procedure for the sec-tors and we present the fugacity expansion coeﬃcients. In Sec.V we calculate the ﬂuctuation ratios using the fu-gacity expansion and extrapolate to small ﬁnite density.Finally in Sec. VI we give a brief summary and outlookfor future work. II. QCD IN THE GRAND CANONICALENSEMBLEA. Susceptibilities and the Taylor expansion

There is a conserved charge corresponding to eachquark ﬂavor of QCD. Working with three ﬂavors, thegrand canonical partition function can be then written interms of three quark number chemical potentials µ u , µ d and µ s . The generalized susceptibilities are deﬁned tobe derivatives of the grand potential (or pressure) withrespect to these chemical potentials: χ udsijk = ∂ i + j + k ( p/T )( ∂ ˆ µ u ) i ( ∂ ˆ µ d ) j ( ∂ ˆ µ s ) k , (1)with the dimensionless chemical potentials ˆ µ X = µ X /T .For the purpose of hadronic phenomenology it is moreconvenient to work with the conserved charges B (baryonnumber), Q (electric charge) and S (strangeness) instead,with chemical potentials µ B , µ Q and µ S , respectively.The basis of µ u , µ d , µ s can be transformed into a basisof µ B , µ Q , µ S with a simple linear transformation, whosecoeﬃcients are given by the B , Q and S charges of theindividual quarks: µ u = 13 µ B + 23 µ Q , (2) µ d = 13 µ B − µ Q , (3) µ s = 13 µ B − µ Q − µ S . (4)Analogously to the case of the quark number chemicalpotentials, the susceptibilities are then deﬁned as χ BQSijk = ∂ i + j + k (cid:0) p/T (cid:1) ∂ ˆ µ iB ∂ ˆ µ jQ ∂ ˆ µ kS . (5)It is straightforward to express the susceptibilities de-ﬁned in Eq. (5) in terms of the coeﬃcients in Eq. (1)[24, 91, 92]. The susceptibilities at µ B = µ S = µ Q = 0are (up to a trivial factorial factor) the Taylor expan-sion coeﬃcients of the pressure near that point. Due tocharge conjugation symmetry, only the even derivativescontribute. In the present study, we always take µ Q = 0and only consider derivatives with respect to µ B and µ S .The Taylor expansion therefore reads: pT = ∞ (cid:88) i =0 ∞ (cid:88) j =0 i ! j ! χ BSij ˆ µ iB ˆ µ jS , (6)where χ BS is just the dimensionless pressure at zerochemical potential.We note that the Taylor expansion is probably themost natural expansion to work within the plasma phaseof QCD. As an exhibit of this, the pressure in the Stefan-Boltzmann (or inﬁnite temperature) limit reads: pT = 8 π

45 + 7 π N f + 12 (cid:88) f (cid:32) µ f T + µ f π T (cid:33) (7)In this approximation, all derivatives above 4th order arezero, and therefore the Taylor expansion is rapidly con-vergent. Calculating corrections to this free gas behaviorin resummed perturbation theory leads to a non-zero, butsmall, value for the sixth-order derivatives [93], leavingthe qualitative conclusion of the fast convergence of theTaylor series in the plasma phase unchanged. B. Fugacity expansion of the free energy

An alternative to the Taylor expansion discussed in theprevious subsection is a Laurent expansion in the fugacityparameters e ˆ µ B and e ˆ µ S near 1. Due to charge conjuga-tion symmetry, a combination e m ˆ µ B + n ˆ µ S and its recip-rocal have the same expansion coeﬃcients, making theLaurent expansion an expansion in hyperbolic cosines: P ( T, ˆ µ B , ˆ µ S ) = (cid:88) j,k P BSjk ( T ) cosh( j ˆ µ B − k ˆ µ S ) . (8)The coeﬃcients P BSjk are also called fugacity expansionor sector coeﬃcients, alluding to the fact that they getcontributions from the Hilbert subspace correspondingto the ﬁxed values of the conserved charges B = j and S = k . In the ideal HRG model, the expansioncoeﬃcients P BS , P BS , P BS , P BS , P BS , P BS all get siz-able contributions from known hadrons and hadron res-onances. In the Boltzmann approximation to the idealHRG, coeﬃcients like P BS are zero, while in the full HRGthey are non-zero, but very small in magnitude and es-sentially negligible.At purely imaginary chemical potentials µ q = iµ Iq ,where the sign problem is absent and lattice simulationscan be performed, we have a Fourier expansion of theform: P ( T, ˆ µ IB , ˆ µ IS ) = (cid:88) j,k P BSjk ( T ) cos( j ˆ µ IB − k ˆ µ IS ) . (9)Diﬀerentiation with respect to the original chemical po-tentials µ B = iµ IB and µ S = iµ IS gives:Im χ BS = (cid:80) j,k jP BSjk ( T ) sin( j ˆ µ IB − k ˆ µ IS ) (10)Im χ BS = (cid:80) j,k ( − k ) P BSjk ( T ) sin( j ˆ µ IB − k ˆ µ IS ) (11) χ BS = (cid:80) j,k j P BSjk ( T ) cos( j ˆ µ IB − k ˆ µ IS ) (12) χ BS = (cid:80) j,k ( − jk ) P BSjk ( T ) cos( j ˆ µ IB − k ˆ µ IS ) (13) χ BS = (cid:80) j,k k P BSjk ( T ) cos( j ˆ µ IB − k ˆ µ IS ) . (14)These formulas and the higher order derivatives of thesewill be used in our ﬁtting procedure, to be described inSection IV. C. The hadron resonance gas and its extensions

In the ideal HRG model the free energy (or pressure)is written as a sum of ideal gas contributions of all knownhadronic resonances H : pT = 1 T (cid:88) H p H = 1 V T (cid:88) H ln Z H ( T, (cid:126)µ ) , (15)with:ln Z H = η H V d H π T (cid:90) ∞ dp p log [1 − η H z H exp ( − (cid:15) H /T )] , (16)where the subscript H indicates dependence on the spe-ciﬁc hadron or hadron resonance in the sum. The rela-tivistic energy is (cid:15) H = (cid:112) p + m H , where m H is the massof the given hadron. The fugacity is z H = exp ( µ H /T ),where the chemical potential associated to H is µ H = µ B B H + µ Q Q H + µ S S H , and the conserved charges B H , Q H and S H are the baryon number, electric charge andstrangeness, respectively. d H is the spin degeneracy, andthe factor η H is 1 for (anti)baryons (fermions) and − χ BQSijk susceptibilities of Eq. (5)can be expressed as: χ BQSijk ( T, ˆ µ B , ˆ µ Q , ˆ µ S ) = (cid:88) H B iH Q jH S kH I Hi + j + k , (17)where the phase space integral at order i + j + k reads: I Hl ( T, ˆ µ B , ˆ µ Q , ˆ µ S ) = ∂ l p H /T ∂ ˆ µ lH . (18)The fugacity expansion coeﬃcients P BS , P BS , P BS , P BS , P BS and P BS can be obtained via the expansionof equation (16) in terms of the modiﬁed Bessel functions K :ln Z H = V T m H d H π ∞ (cid:88) n =1 ( − η H ) n +1 z nH n K (cid:16) nm H T (cid:17) . (19)The Boltzmann approximation consists of taking only the n = 1 term in the above expansion, which accounts forthe lowest order in the fugacity parameters. In the Boltz-mann approximation, the sectors read: P BSjk = (cid:88) H δ B H ,j δ S H ,k d H m H π T K (cid:16) m H T (cid:17) . (20)In the full ideal HRG, a hadron with B H = 1 and S H = 0will also give contributions to the higher order sectors,such as P BS and P BS , due to the terms n = 2 and n = 3in Eq. (19), respectively. These are, however, exponen-tially suppressed due to the behavior of the Bessel func-tion K ( x ) ∼ (cid:112) π x e − x as x → ∞ . These contributionsare orders of magnitude smaller than the full weight ofthe respective sectors as obtained from the lattice.The HRG model is an approximation of the more gen-eral formula by Dashen, Bernstein and Ma, which givesthe fugacity expansion coeﬃcients in terms of the S-matrix: P BSjk = 1 π T (cid:90) ∞ M BSjk dEE K (cid:18) ET (cid:19) i Tr B = j,S = k (cid:18) S † dSdE − dS † dE S (cid:19) c , (21)where M BSjk is the mass threshold for the B = j , S = k channel, the trace is taken over this Hilbert subspace,and the subscript c signiﬁes that only connected S-matrixelements are to be taken. For the speciﬁc case of elastic2 → i Tr B = j,S = k (cid:18) S † dSdE − dS † dE S (cid:19) c → (cid:88) J (2 J + 1) (cid:18) dδ J,I =0 dE + 3 dδ J,I =1 dE (cid:19) , (22)where the δ J,I are the scattering phase shifts for angu-lar momentum J and isospin I and the isospin singletand triplet contributions have been written separately.After integration by parts with respect to E , we get tothe conclusion that the contribution of elastic scatter-ing is given by the integral of the phase shift with anexponential weight. This leads to the expectation thatdominantly repulsive interactions will lead to a negativesub-leading fugacity expansion coeﬃcient. This fact wasexploited when constructing repulsive core HRG modelsand comparing them with lattice data in Refs. [42, 74]. Itis also reasonable to expect, due to the exponential sup-pression of the K Bessel functions, that in the hadronicphase there will be a strong hierarchy of the fugacity ex-pansion coeﬃcients with increasing quantum numbers, soe.g. P BS (cid:29) P BS (cid:29) P BS as well as P BS (cid:29) P BS (cid:29) P BS and P BS (cid:29) P BS (cid:29) P BS etc. It is thus a reasonable ex-pectation that, in the hadronic phase, the fugacity expan-sion will converge faster than the Taylor expansion. Thisis the opposite situation as in the plasma phase, wherethe Taylor expansion converges quickly, while the fugac-ity expansion converges slowly. This makes the fugacityexpansion particularly useful for modelling the hadronicphase, and therefore also for the study of chemical freeze-out in heavy ion collisions.Here we do not utilize any S-matrix formula or anymean ﬁeld approximation thereof, but rather calculatethe sub-leading sector coeﬃcients P BS , P BS , P BS , P BS etc. directly from lattice simulations. III. LATTICE SETUP

We use a staggered fermion action with 4 steps of stoutsmearing [94] with the smearing parameter ρ = 0 .

145 MeV 150 MeV 155 MeV 160 MeV40 × µ I = 0 10348 10520 10345 1161132 × µ I = 0 8518 8461 1695 917424 × µ I = 0 40247 39996 19953 2001536 × µ I (cid:54) = 0 146968 154479 153513 14416932 × µ I (cid:54) = 0 124915 81814 300779 26464724 × µ I (cid:54) = 0 184896 171224 166034 161454TABLE I. Number of evaluated conﬁgurations on the variouslattices and temperatures. The µ I (cid:54) = 0 statistics is distributedover 143 pairs of imaginary strange and baryon chemical po-tentials. combination was ﬁrst used in Ref. [24], where informationabout the line of constant physics can be found. Forthe scale setting we use the pion decay constant f π =130 . N τ =8 ,

10 and 12 to perform an estimation of the continuumvalue of our observables. The spatial extent of the latticeis given by the aspect ratio LT ≈

3. Due to technicalreasons, some lattices had slightly diﬀerent values for thisratio, as given in Table. I. Given the error bars on theﬁnal results, we did not optimize this further.For the continuum extrapolations, we assume a lin-ear scaling in 1 /N τ . Since taste breaking eﬀects are stillrather large on these lattices, we only call our resultscontinuum estimates, as opposed to fully controlled con-tinuum extrapolations, e.g. when N τ = 16 is part of theextrapolation.For all values of N τ , we use simulations at four dif-ferent temperatures T = 145 MeV ,

150 MeV ,

155 MeVand 160 MeV. At each temperature and each latticespacing, we perform a two-dimensional scan in the imag-inary chemical potentials µ IB and µ IS , with the chem-ical potentials taking the values ( µ IB , µ IS ) = π ( i, j ),with i = 0 , , . . . ,

15 and j = 0 , , . . . ,

9, for a total of9 ×

16 = 144 simulation points. In each µ i (cid:54) = 0 point,we simulated one Rational Hybrid Monte Carlo streamwith several thousand trajectories, evaluating every ﬁfthconﬁguration for the ﬂuctuation observables as detailedin Ref. [24]. Our statistics is summarized in Table I.The statistical errors are calculated using the jackknifemethod. The estimation of the systematic errors is amore elaborate process. Ambiguities appear at variouspoints of the analysis, e.g. in the way the continuumextrapolation is calculated, or how many ﬁt parameterswe use for the extraction of the fugacity expansion coef-ﬁcients. We consider all combinations of the possibilitiesand take the spread of the results as systematic error. IV. FUGACITY EXPANSION COEFFICIENTS

The estimation of the coeﬃcients P BSij proceedsthrough a correlated ﬁt. On the µ B = µ S = 0 ensem-bles, the ﬂuctuations χ BS , χ BS , χ BS , χ BS , χ BS , χ BS , χ BS and χ BS are included, while for the other ensembles we × , 4stout, f π scale T = 145 MeV χ red ( B max = 2) = 1 . χ red ( B max = 3) = 1 . T = 150 MeV χ red ( B max = 2) = 1 . χ red ( B max = 3) = 1 . T = 155 MeV χ red ( B max = 2) = 0 . χ red ( B max = 3) = 0 . P BS P BS P BS P BS P BS P BS , − P BS P BS P BS P BS P BS P BS P BS P BS P BS P BS T = 160 MeV χ red ( B max = 2) = 1 . χ red ( B max = 3) = 1 . | P B S i j | B max = 2 ,+ B max = 2 ,- B max = 3 ,+ B max = 3 ,- | P B S i j || P B S i j || P B S i j | − . − . − . χ /N dof ≈ / − . − . . χ /N dof ≈ . / − . − . − . . χ /N dof ≈ / − . . . χ /N dof ≈ . / − . − . − . − .

145 150 155 160 χ /N dof ≈ . / P B S N τ = 8 N τ = 10 N τ = 12 cont. estim. P B S P B S P B S P B S , − T [ MeV ] FIG. 1. Left: The four panels refer to four temperatures, each showing the obtained coeﬃcients of the fugacity expansion on alogarithmic scale (negative values are shown in blue). Only the leftmost ﬁve are accounted for by ideal HRG. The next sevenappear in the next order, which we use later for phenomenology. For the next order (last four coeﬃcients with B = 3) we seeno stable signal. There are two symbols per coeﬃcient, triangles for the complete ﬁt and circles for the case without the B = 3part. The data refer to our ﬁnest lattice spacing. Right: Example for the combined continuum extrapolation of the extractedcoeﬃcients. use Im χ BS and Im χ BS . This leads to a block-diagonalcovariance matrix with one 8 × µ = 0 ensemble, and 143 blocks of size 2 ×

2, cor-responding to the ensembles with a non-zero value of atleast one of the chemical potentials. The covariance ma-trix blocks are estimated by the jackknife method with 24jackknife samples. The truncation of the fugacity expan-sion is somewhat ambiguous, as there is no single smallparameter in which we actually perform this expansion.To estimate systematic errors coming from the choice ofthe ansatz, we therefore perform two ﬁts for each en-semble, for which we introduce the shorthand notations B max = 2 and B max = 3. The sectors included in the B max = 2 analysis are: P BS , P BS , P BS , P BS ,P BS , P BS , − , P BS , P BS ,P BS , P BS , P BS , P BS . (23)The ﬁrst ﬁve of these correspond to sectors that are al-ready present in the ideal HRG in the Boltzmann ap-proximation. They also set contributions from interac-tions though, e.g. non-resonant pion-nucleon interactionscontribute to P BS , while K -Λ interactions contribute to P BS . The P BS , − is an exotic (pentaquark) sector. Itgets contributions for example for valence quark content uudd ¯ s which corresponds e.g. to p + K scattering. Thesectors P BS i get contributions from baryon-baryon scat-tering: P BS from N − N , P BS from N − Λ, P BS from N − Ξ or Λ − Λ and ﬁnally P BS from N − Ω or Λ − Ξ.In each case, Σ can replace Λ. The coeﬃcients P BS and P BS get contributions from two- and three-kaon scatter-ing, respectively. The inclusion of the P BS sector withthe omission of the P BS sector is motivated by the lowermass threshold of 3 kaon scattering as compared to 3baryon scattering. In addition to these, we also per-formed an analysis were four sectors with B = 3 wereadded: P BS , P BS31 , P BS32 , P BS33 . (24)These get contributions from three-baryon scattering,with various strangeness contents. Our data is not yetsuﬃciently accurate to obtain a reliable estimate of thesectors with B = 3, and the inclusion of these sectorsdoes not improve the χ of the ﬁts. Whether we includethese, or not, the results for the B = 2 sectors remainconsistent, as we show in the left panel of Fig. 1, wherethe sector coeﬃcients from the two diﬀerent ﬁts on the N τ = 12 lattices are shown. This way we demonstratethe stability of the sectors included in the B max = 2 set.Only at the highest temperature T = 160MeV, and onlyfor one sector, P BS , is the systematic error coming fromincluding the B = 3 sectors comparable to the statisticalerror of the ﬁts.The continuum limit estimation of the sectors proceedsthrough a combined ﬁt in temperature and lattice spacing(or equivalently N τ ) via the ansatz f ( T, N τ ) = (cid:0) a + a T + a T (cid:1) + (cid:0) b + b T + b T (cid:1) N τ . (25)For the systematic error we compare this ansatz withand without the coeﬃcient b . The continuum extrap-olation of the beyond-ideal-HRG sectors for the case ofthe B max = 2 ﬁts at ﬁxed T and N τ and b kept as afree parameter, is shown in Fig 1 (right). The other 3ﬁts look quantitatively similar. All of the continuum ﬁtshave acceptable ﬁt quality, with Q values over 1%. Asa conservative estimate of the systematics, we combinethem with uniform weights. As can be seen in the rightpanel of Fig 1, the slopes of the continuum extrapolationsof all beyond-ideal-HRG sectors appear to be mild, ex-cept for the sector P BS , which corresponds to kaon-kaonscattering, and changes its sign during the continuum ex-trapolation. As expected, this sector - being related tokaons - suﬀers from relatively large taste-breaking eﬀects.The ﬁnal results for the beyond-ideal-HRG sectors canbe seen in Fig. 2. Within the statistical precision of ourresults, P BS is roughly the same as P BS , while P BS issmaller than the previous two. As a comparison, theideal HRG model prediction for the sum (cid:80) k P BS k at T =155MeV is of the order 10 − , orders of magnitude lowerthan what we see here. The two-kaon scattering sector P BS goes slightly below zero at around 155MeV within1 σ uncertainty. The three-kaon sector P BS is consistent − . − . . − . − . . − . − . . − . − . . − . − . .

145 150 155 160 P B S combinedT-by-T P B S P B S P B S P B S , − T [ MeV ] FIG. 2. Our continuum estimates of the beyond-ideal-HRGsector coeﬃcients. We show the results both from our com-bined temperature and continuum ﬁt (green bands) and of a T -by- T continuum limit extrapolation (blue points). System-atic errors are included, in the ﬁrst case by varying B max = 2vs 3 and b = 0 or b (cid:54) = 0, and in the second case by varying B max = 2 vs 3. with zero in the entire temperature range and is thereforenot included in the plot. An upper limit on its magnitudewith 1 σ uncertainty is 2 · − . The exotic sector P BS , − is rather large, consistently with our earlier, statisticallyindependent ﬁnding in Ref. [87] on N τ = 12 lattices. Wehave already published the leading sector coeﬃcients forwhich the ideal HRG has a prediction in the Boltzmann-approximation, namely P BS , P BS , P BS , P BS , P BS inRef. [41]. We will not repeat the results for those sectorshere. V. FLUCTUATION-RATIOS AT FINITEBARYON DENSITY

Having the coeﬃcients of the fugacity expansion, thethermodynamics can be readily obtained. In this sectionwe calculate the baryon number and strangeness ﬂuctu-ations and their ratios in the studied temperature range.There is no diﬃculty in evaluating Eq. 8 and its µ B - and µ S -derivatives at any chemical potential.Heavy ion collisions involving lead or gold atoms cor-respond to the conditions χ S = 0 and χ B = 0 . χ Q . Forthe purposes of the present study, we impose strangenessneutrality and leave the second conditions for futurework. In fact, we use the simpliﬁed form χ B = 0 . χ Q ,which is realized at vanishing electric charge chemicalpotential.To show the magnitude of the cut-oﬀ eﬀects, we startwith a pair of quantities at µ B = 0. The fugacity expan-sion, and more generally imaginary chemical potentialsimulations oﬀer an eﬃcient way to calculate suscepti-bilities at µ B = 0. In Ref. [44] we calculated the ratios χ B /χ B and χ B /χ B at a ﬁnite lattice spacing of N τ = 12with the same lattice action used here. Here we showcontinuum estimates of the ﬂuctuation ratios χ B /χ B and χ B /χ B at µ B = 0, together with the data at ﬁnite N τ inFig. 3. In the ideal HRG model χ B /χ B = 1 for all tem-peratures, meaning that above T = 150MeV our resultsshow a clear deviation from the HRG prediction, dueto presence of the non-zero beyond-ideal-HRG sectors.The diﬀerence between the two ratios χ B /χ B − χ B /χ B is also shown. In the µ S = 0 case the two ratios at µ B = 0 are identical. The diﬀerence between the tworatios comes from imposing the strangeness neutralitycondition χ S = 0. This diﬀerence also shows mild cut-oﬀeﬀects.After the sectors are obtained, we perform extrapo-lations to real chemical potentials using the ansatz ofEq. (8) truncated at the B max = 2 level. We extrapo-late ﬁrst at ﬁxed T and N τ . We consider the ﬂuctua-tion ratios χ B /χ B , χ B /χ B , χ B /χ B and χ BS /χ S on thestrangeness-neutral line χ S = 0, which determines µ S asa function of µ B . While the extrapolation always usesthe 12 sectors of the B max = 2 level, the values of thesesectors are taken both from the B max = 2 and B max = 3ﬁts, to estimate the systematic errors. We then performa continuum estimation at ﬁxed values of µ B /T with thesame combined T and N τ ﬁt as in the case of the baryonand strangeness sectors.We had one sector, P BS , with a steep continuum ex-trapolation. Should we expect additional systematic er-rors coming from the non-trivial continuum scaling in thephenomenology? The answer is no, as we demonstrate inFig. 4. The multi-kaon sectors do not contribute to thebaryon ﬂuctuations. The only ratio with phenomeno-logical relevance where the P BS may be important is χ BS /χ S . We calculated this ratio with and without themulti-kaon sectors and compared the results in Fig. 4 at T = 160 MeV. Although at this temperature and this . . . . . . . . − . . . . .

145 150 155 160 χ B / χ B ( T , µ B = ) N τ = 8 N τ = 10 N τ = 12 comb. cont. est. χ B / χ B ( T , µ B = ) χ B / χ B − χ B / χ B T [ MeV ] FIG. 3. Fluctuation-ratios χ B /χ B and χ B /χ B obtained fromthe fugacity expansion truncated at the B max = 2 level at µ B = 0 on our N τ = 8 ,

10 and 12 lattices and our continuumestimates from these data. The points at ﬁnite lattice spacinginclude a systematic error coming from whether we used the12- or the 16-parameter ﬁt to determine the B max = 2 sectors.The continuum results include systematic error from 4 ﬁts,in addition for the 12 vs 16 parameter ﬁts at ﬁxed N τ and T we also include a 5- vs 6-parameter combined T and N τ continuum ﬁt. lattice spacing we have the largest value for this diﬃcultsector, we see hardly any signiﬁcant eﬀect from it, espe-cially not when compared to the statistical errors afterthe continuum extrapolation step.The ﬁnal results for the ﬂuctuation ratios χ B /χ B , χ B /χ B , χ B /χ B , χ BS /χ S on the strangeness neutral linecan be seen in Fig. 5. The ﬁrst of these ratio is stronglydependent on the chemical potential, but not on the tem-perature, making it a proxy of the chemical potential,at least for small values of µ B . On has to rememberthough that if the critical endpoint exists, the ﬂuctua-tion χ B diverges there, leading to χ B /χ B → µ B .The other three are more strongly dependent on the − . − . − . − . − . − . − . − .

26 0 0 . . . . . × , 4stout, T=160MeV, f π scale χ B S / χ S µ B /T full B max = 2 expansiondropping P BS and P BS FIG. 4. Multi-kaon interactions have a negligible impact onthe ﬂuctuations ratios studied in this work. Here we show theeﬀect of dropping the multi-kaon sectors from the fugacity ex-pansion, demonstrated on the N τ = 8 data, at T = 160MeV,where they are the largest in our temperature range. Notealso that, on the N τ = 8 lattices, the magnitude of the multi-kaon sectors is larger than our continuum estimate for them. temperature and less strongly on the chemical potential,therefore making them possible proxies for the temper-ature. The ratios χ B /χ B and χ B /χ B can be regardedas a baryon thermometer, while the ratio χ BS /χ S as astrangeness-related one. This latter ratio is of large phe-nomenological interest, as experimental net-lambda andnet-kaon ﬂuctuations can be used to construct the ra-tio σ / ( σ + σ K ). It was shown in Ref. [87] that thisis a good experimental proxy of χ BS /χ S , not stronglyaﬀected by experimental eﬀects, which makes it a primetarget for comparison with experiments.We show the ﬂuctuation ratios in Fig. 5 as functionsof the dimensionless chemical potential µ B /T at a fewvalues of the temperature, as well as on the crossoverline calculated to order µ B in Ref. [45]: T c ( µ B ) ≈ T c (cid:0) − κ ˆ µ B (cid:1) , (26)with T c = 158 . ± . κ = 0 . ± . σ bands on the T c ( µ B ) line al-ways overlap with the 1 σ bands for a ﬁxed T = T c =158MeV.Our results on the ratios χ B /χ B and χ B /χ B are alsoshown as a function of χ B /χ B in Fig. 6. Within 1 σ ,our results are consistent with recent lattice results onthe susceptibility ratios using the Taylor method [30].In Fig. 6 we also compare to the data on net-protonﬂuctuations from the STAR collaboration [88–90], thenet-proton skewness-to-mean ratio C /C and the net-proton kurtosis-to-variance ratio C /C , as functions ofthe net-proton mean-to-variance ratio C /C at chemi-cal freeze-out. The advantage of using these variables inthe comparison is that is does not involve any modelling . . . . . . . . . . . . . . . . . . . − . − . − . − . − . − . χ B / χ B MeV

MeV T c ( µ ) χ B / χ B MeV

MeV T c ( µ ) χ B / χ B MeV

MeV T c ( µ ) χ B S / χ S µ B /T MeV

MeV

MeV T c ( µ ) FIG. 5. Continuum estimates of the ﬂuctuation ratios χ B /χ B , χ B /χ B , χ B /χ B , χ BS /χ S from a fugacity expansion truncatedat the B max = 2 level, shown as a function of the dimension-less chemical potential µ B /T for ﬁxed temperatures, as wellas on the crossover line T c ( µ B ). of the freeze-out conditions, other than assuming thatchemical freeze-out happens close to the QCD crossoveron the phase diagram. Our results are consistent withthe experimental results. While such a direct compari-son suﬀers from many caveats [87, 96–100], the similarityin the trends supports the idea that experimentally ob-served net-proton ﬂuctuation ratios reﬂect with some ac-curacy the thermal ﬂuctuations in an equilibrated QCDmedium. VI. SUMMARY AND OUTLOOK

We have calculated fugacity expansion coeﬃcients ofthe QCD pressure beyond the ideal HRG model, sep-arating contributions to the QCD free energy coming0 . . . . . √ s NN /GeV =

200 62.4 54.4 39 27 χ B n + / χ B n o r C n + / C n χ B /χ B or C /C µ B /T Fugacity expansion, χ B /χ B Fugacity expansion, χ B /χ B STAR, C /C STAR, C /C FIG. 6. Our continuum estimates of the ﬂuctuation ratios χ B /χ B and χ B /χ B compared with STAR data on net-protonﬂuctuations from Ref. [88]. Both the values of the ﬂuctuationsand the trend as a function of baryon density are consistent. from Hilbert subspaces with diﬀerent values of the baryonnumber and strangeness quantum numbers. This allowsone to quantify the importance of processes like kaon-kaon and baryon-baryon scattering, when modelling theQCD medium in the hadronic phase, but close to thecrossover. We estimated the continuum value of thesecoeﬃcients with lattice simulations of temporal extent N τ = 8 ,

10 and 12 using the staggered discretization. Weobserved large cut-oﬀ eﬀects in the kaon and multi-kaonsectors, only. To study these, and to make the continuumextrapolation more robust, the inclusion of ﬁner latticesis desirable.Note that our study was limited to an aspect ratio of LT ≈

3, future studies should also investigate ﬁnite vol-ume eﬀects in the baryon-strangeness sectors. However,a strong volume dependence is more likely to be observedin correspondence with the electric charge sectors. Thefull picture of baryon interactions in the hadronic phasewill emerge from the three-dimensional mapping of the µ B − µ S − µ Q space. In this work we restricted the spaceto µ Q = 0. The sectors we obtained are sums of variouscharge sectors, P BSij = (cid:80) k P BSQijk , and we cannot diﬀer-entiate between the terms, though it would be possiblein a more elaborate setup. Still, the level of separationachieved in this work already provides plenty of new in-formation for hadronic modelling of the QCD medium.We used the truncated fugacity expansion to calcu-late phenomenologically relevant ﬂuctuation ratios on thestrangeness neutral line both as a function of the chemicalpotential and as a function of the baryon number mean-to-variance ratio, which can be regarded as a proxy ofthe baryo-chemical potential. While a direct comparisonis by no means trivial, the fugacity expansion coeﬃcientsappear to describe the trend in the STAR data on net-proton ﬂuctuations [88–90]. It has been pointed out in the literature, that the mod-iﬁcations of the ideal HRG model that include the eﬀectsof global baryon number conservation lead to a reduc-tion in higher order ﬂuctuations and therefore describethe experimental data better. For a recent study of theseeﬀects see [99]. In a recent paper, it was also pointed outthat the decrease in the kurtosis with increasing chemicalpotential observed in the data is likely not due to criticalpoint eﬀects [101]. The corrections to ideal HRG studiedin this work have a similar magnitude as the experimen-tal eﬀects like canonical eﬀects and volume ﬂuctuations.Since both baryon interactions and the global conserva-tion laws appear to push the χ B /χ B and χ B /χ B ratiosdown, it is important to have an estimate of the relativesize of these types of eﬀects under realistic conditions.Performing a study of this kind is an important task forthe near future, as it will guide the correct interpretationof STAR data. 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 , 443:675–678, 2006.[2] Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo. TheQCD transition temperature: Results with physicalmasses in the continuum limit.

Phys. Lett. B , 643:46–54,2006.[3] Szabolcs Borsanyi, Zoltan Fodor, Christian Hoelbling,Sandor D Katz, Stefan Krieg, Claudia Ratti, andKalman K. Szabo. Is there still any T c mystery in latticeQCD? Results with physical masses in the continuumlimit III. JHEP , 09:073, 2010.[4] A. Bazavov et al. The chiral and deconﬁnement aspectsof the QCD transition.

Phys. Rev. D , 85:054503, 2012.[5] Szabolcs Borsanyi, Zoltan Fodor, Christian Hoelbling,Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo.Full result for the QCD equation of state with 2+1 ﬂa-vors.

Phys. Lett. B , 730:99–104, 2014.[6] A. Bazavov et al. Equation of state in ( 2+1 )-ﬂavorQCD.

Phys. Rev. D , 90:094503, 2014.[7] Sz. Borsanyi et al. Calculation of the axion mass basedon high-temperature lattice quantum chromodynamics.

Nature , 539(7627):69–71, 2016.[8] A. Bazavov et al. Update on the 2+1+1 Flavor QCDEquation of State with HISQ.

PoS , LATTICE2013:154,2014.[9] Kenji Fukushima and Chihiro Sasaki. The phase dia-gram of nuclear and quark matter at high baryon den-sity.

Prog. Part. Nucl. Phys. , 72:99–154, 2013.[10] Peter Kov´acs, Zsolt Sz´ep, and Gy¨orgy Wolf. Existenceof the critical endpoint in the vector meson extendedlinear sigma model.

Phys. Rev. D , 93(11):114014, 2016.[11] Renato Critelli, Jorge Noronha, Jacquelyn Noronha-Hostler, Israel Portillo, Claudia Ratti, and RomuloRougemont. Critical point in the phase diagram ofprimordial quark-gluon matter from black hole physics.

Phys. Rev. D , 96(9):096026, 2017.[12] Philippe de Forcrand. Simulating QCD at ﬁnite density.

PoS , LAT2009:010, 2009.[13] A. Hasenfratz and D. Toussaint. Canonical ensemblesand nonzero density quantum chromodynamics.

Nucl.Phys. B , 371:539–549, 1992.[14] Z. Fodor and S. D. Katz. A New method to study lat-tice QCD at ﬁnite temperature and chemical potential.

Phys. Lett. B , 534:87–92, 2002.[15] Z. Fodor and S. D. Katz. Critical point of QCD at ﬁ-nite T and mu, lattice results for physical quark masses.

JHEP , 04:050, 2004.[16] Z. Fodor and S. D. Katz. Lattice determination of thecritical point of QCD at ﬁnite T and mu.

JHEP , 03:014,2002.[17] Matteo Giordano, Kornel Kapas, Sandor D. Katz,Daniel Nogradi, and Attila Pasztor. Eﬀect of stoutsmearing on the phase diagram from multiparam-eter reweighting in lattice QCD.

Phys. Rev. D ,102(3):034503, 2020.[18] Matteo Giordano, Kornel Kapas, Sandor D. Katz,Daniel Nogradi, and Attila Pasztor. New approach tolattice QCD at ﬁnite density; results for the critical endpoint on coarse lattices.

JHEP , 05:088, 2020. [19] 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.[20] 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.[21] R. V. Gavai and Sourendu Gupta. QCD at ﬁnitechemical potential with six time slices.

Phys. Rev. ,D78:114503, 2008.[22] Szabolcs Bors´anyi, Zolt´an Fodor, S´andor D. Katz, Ste-fan Krieg, Claudia Ratti, and K´alm´an Szab´o. Fluctu-ations of conserved charges at ﬁnite temperature fromlattice QCD.

JHEP , 01:138, 2012.[23] Sz. Bors´anyi, G. Endr˝odi, Z. Fodor, S.D. Katz, S. Krieg,et al. QCD equation of state at nonzero chemical po-tential: continuum results with physical quark massesat order mu . JHEP , 1208:053, 2012.[24] R. Bellwied, S. Bors´anyi, Z. Fodor, S. D. Katz,A. P´asztor, C. Ratti, and K. K. Szab´o. Fluctuationsand correlations in high temperature QCD.

Phys. Rev. ,D92(11):114505, 2015.[25] H. T. Ding, Swagato Mukherjee, H. Ohno, P. Petreczky,and H. P. Schadler. Diagonal and oﬀ-diagonal quarknumber susceptibilities at high temperatures.

Phys.Rev. , D92(7):074043, 2015.[26] A. Bazavov et al. The QCD Equation of State to O ( µ B )from Lattice QCD. Phys. Rev. D , 95(5):054504, 2017.[27] Zolt´an Fodor, Matteo Giordano, Jana N. G¨unther,Korn´el Kap´as, S´andor D. Katz, Attila P´asztor, IsraelPortillo, Claudia Ratti, D´enes Sexty, and K´alman K.Szab´o. Trying to constrain the location of the QCDcritical endpoint with lattice simulations.

Nucl. Phys.A , 982:843–846, 2019.[28] Matteo Giordano and Attila P´asztor. Reliable estima-tion of the radius of convergence in ﬁnite density QCD.

Phys. Rev. D , 99(11):114510, 2019.[29] Matteo Giordano, Kornel Kapas, Sandor D. Katz,Daniel Nogradi, and Attila Pasztor. Radius of conver-gence in lattice QCD at ﬁnite µ B with rooted staggeredfermions. Phys. Rev. D , 101(7):074511, 2020.[30] A. Bazavov et al. Skewness, kurtosis, and the ﬁfth andsixth order cumulants of net baryon-number distribu-tions from lattice QCD confront high-statistics STARdata.

Phys. Rev. D , 101(7):074502, 2020.[31] Philippe de Forcrand and Owe Philipsen. The QCDphase diagram for small densities from imaginary chem-ical potential.

Nucl. Phys. , B642:290–306, 2002.[32] Massimo D’Elia and Maria Paola Lombardo. Finite den-sity QCD via imaginary chemical potential.

Phys. Rev. ,D67:014505, 2003.[33] Massimo D’Elia and Francesco Sanﬁlippo. Thermody-namics of two ﬂavor QCD from imaginary chemical po-tentials.

Phys. Rev. , D80:014502, 2009.[34] Paolo Cea, Leonardo Cosmai, and Alessandro Papa.Critical line of 2+1 ﬂavor QCD.

Phys. Rev. ,D89(7):074512, 2014.[35] Claudio Bonati, Philippe de Forcrand, Massimo D’Elia,Owe Philipsen, and Francesco Sanﬁlippo. Chiral phasetransition in two-ﬂavor QCD from an imaginary chem- ical potential. Phys. Rev. , D90(7):074030, 2014.[36] Paolo Cea, Leonardo Cosmai, and Alessandro Papa.Critical line of 2+1 ﬂavor QCD: Toward the continuumlimit.

Phys. Rev. , D93(1):014507, 2016.[37] Claudio Bonati, Massimo D’Elia, Marco Mariti, MicheleMesiti, Francesco Negro, and Francesco Sanﬁlippo. Cur-vature of the chiral pseudocritical line in QCD: Contin-uum extrapolated results.

Phys. Rev. , D92(5):054503,2015.[38] R. Bellwied, S. Bors´anyi, Z. Fodor, J. G¨unther, S. D.Katz, C. Ratti, and K. K. Szab´o. The QCD phase dia-gram from analytic continuation.

Phys. Lett. , B751:559–564, 2015.[39] Massimo D’Elia, Giuseppe Gagliardi, and FrancescoSanﬁlippo. Higher order quark number ﬂuctuations viaimaginary chemical potentials in N f = 2+1 QCD. Phys.Rev. D , 95(9):094503, 2017.[40] J. N. G¨unther, R. Bellwied, S. Bors´anyi, Z. Fodor, S. D.Katz, A. P´asztor, C. Ratti, and K. K. Szab´o. The QCDequation of state at ﬁnite density from analytical con-tinuation.

Nucl. Phys. , A967:720–723, 2017.[41] Paolo Alba et al. Constraining the hadronic spectrumthrough QCD thermodynamics on the lattice.

Phys.Rev. D , 96(3):034517, 2017.[42] Volodymyr Vovchenko, Attila Pasztor, Zoltan Fodor,Sandor D. Katz, and Horst Stoecker. Repulsive baryonicinteractions and lattice QCD observables at imaginarychemical potential.

Phys. Lett. B , 775:71–78, 2017.[43] Claudio Bonati, Massimo D’Elia, Francesco Negro,Francesco Sanﬁlippo, and Kevin Zambello. Curva-ture of the pseudocritical line in QCD: Taylor ex-pansion matches analytic continuation.

Phys. Rev. ,D98(5):054510, 2018.[44] Szabolcs Bors´anyi, Zolt´an Fodor, Jana N. G¨unther,S´andor K. Katz, K´alm´an K. Szab´o, Attila P´asztor, Is-rael Portillo, and Claudia Ratti. Higher order ﬂuctua-tions and correlations of conserved charges from latticeQCD.

JHEP , 10:205, 2018.[45] Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther,Ruben Kara, Sandor D. Katz, Paolo Parotto, AttilaPasztor, Claudia Ratti, and Kalman K. Szabo. QCDCrossover at Finite Chemical Potential from LatticeSimulations.

Phys. Rev. Lett. , 125(5):052001, 2020.[46] Attila P´asztor, Zsolt Sz´ep, and Gergely Mark´o. Ap-parent convergence of Pad \ ’e approximants for thecrossover line in ﬁnite density QCD. 10 2020.[47] R. Hagedorn. Statistical thermodynamics of strong in-teractions at high-energies. Nuovo Cim. Suppl. , 3:147–186, 1965.[48] Roger Dashen, Shang-Keng Ma, and Herbert J. Bern-stein. S Matrix formulation of statistical mechanics.

Phys. Rev. , 187:345–370, 1969.[49] R. F. Dashen and R. Rajaraman. Narrow Resonancesin Statistical Mechanics.

Phys. Rev. D , 10:694, 1974.[50] R. F. Dashen and R. Rajaraman. Eﬀective Elementarityof Resonances and Bound States in Statistical Mechan-ics.

Phys. Rev. D , 10:708, 1974.[51] R. Venugopalan and M. Prakash. Thermal properties ofinteracting hadrons.

Nucl. Phys. A , 546:718–760, 1992.[52] F. Karsch, K. Redlich, and A. Tawﬁk. Hadron resonancemass spectrum and lattice QCD thermodynamics.

Eur.Phys. J. C , 29:549–556, 2003.[53] Pasi Huovinen and Pter Petreczky. QCD Equationof State and Hadron Resonance Gas.

Nucl. Phys. A , 837:26–53, 2010.[54] Szabolcs Borsanyi, Gergely Endrodi, Zoltan Fodor, An-tal Jakovac, Sandor D. Katz, Stefan Krieg, ClaudiaRatti, and Kalman K. Szabo. The QCD equation ofstate with dynamical quarks.

JHEP , 11:077, 2010.[55] Abhijit Majumder and Berndt Muller. Hadron MassSpectrum from Lattice QCD.

Phys. Rev. Lett. ,105:252002, 2010.[56] A. Bazavov et al. Additional Strange Hadrons fromQCD Thermodynamics and Strangeness Freezeout inHeavy Ion Collisions.

Phys. Rev. Lett. , 113(7):072001,2014.[57] Pok Man Lo, Micha(cid:32)l Marczenko, Krzysztof Redlich, andChihiro Sasaki. Matching the Hagedorn mass spectrumwith Lattice QCD results.

Phys. Rev. C , 92(5):055206,2015.[58] W. Weinhold, B. Friman, and W. Norenberg. Thermo-dynamics of Delta resonances.

Phys. Lett. B , 433:236–242, 1998.[59] Bengt Friman, Pok Man Lo, Micha(cid:32)l Marczenko,Krzysztof Redlich, and Chihiro Sasaki. Strangenessﬂuctuations from K − π interactions. Phys. Rev. D ,92(7):074003, 2015.[60] Pok Man Lo, Bengt Friman, Krzysztof Redlich, andChihiro Sasaki. S-matrix analysis of the baryon electriccharge correlation.

Phys. Lett. B , 778:454–458, 2018.[61] Richard A. Arndt, John S. Hyslop, III, and L. DavidRoper. Nucleon-Nucleon Partial Wave Analysis to 1100-MeV.

Phys. Rev. D , 35:128, 1987.[62] R. A. Arndt, W. J. Briscoe, I. I. Strakovsky, and R. L.Workman. Updated analysis of NN elastic scattering to3-GeV.

Phys. Rev. C , 76:025209, 2007.[63] Ron L. Workman, William J. Briscoe, and Igor I.Strakovsky. Partial-Wave Analysis of Nucleon-NucleonElastic Scattering Data.

Phys. Rev. C , 94(6):065203,2016.[64] Henk Polinder, Johann Haidenbauer, and Ulf-G. Meiss-ner. Hyperon-nucleon interactions: A Chiral eﬀectiveﬁeld theory approach.

Nucl. Phys. A , 779:244–266,2006.[65] J. Haidenbauer, S. Petschauer, N. Kaiser, U. G. Meiss-ner, A. Nogga, and W. Weise. Hyperon-nucleon inter-action at next-to-leading order in chiral eﬀective ﬁeldtheory.

Nucl. Phys. A , 915:24–58, 2013.[66] Shreyasi Acharya et al. Unveiling the strong interactionamong hadrons at the LHC.

Nature , 588:232–238, 2020.[67] Shreyasi Acharya et al. Measurement of strangebaryon–antibaryon interactions with femtoscopic corre-lations.

Phys. Lett. B , 802:135223, 2020.[68] Laura Fabbietti, Valentina Mantovani Sarti, andOt´on V´azquez Doce. Hadron-hadron interactions mea-sured by ALICE at the LHC. 12 2020.[69] Takumi Doi et al. Baryon interactions from lattice QCDwith physical masses – Overview and S = 0 , − PoS , LATTICE2016:110, 2017.[70] Kenji Sasaki et al. ΛΛ and NΞ interactions from Lat-tice QCD near the physical point.

Nucl. Phys. A ,998:121737, 2020.[71] Takumi Iritani et al. N Ω dibaryon from lattice QCDnear the physical point.

Phys. Lett. B , 792:284–289,2019.[72] Volodymyr Vovchenko, Mark I. Gorenstein, and HorstStoecker. van der Waals Interactions in Hadron Reso-nance Gas: From Nuclear Matter to Lattice QCD.

Phys. Rev. Lett. , 118(18):182301, 2017.[73] Paolo Alba, Wanda Maria Alberico, Alessandro Nada,Marco Panero, and Horst St¨ocker. Excluded-volume ef-fects for a hadron gas in Yang-Mills theory.

Phys. Rev.D , 95(9):094511, 2017.[74] Pasi Huovinen and Peter Petreczky. Hadron resonancegas with repulsive interactions and ﬂuctuations of con-served charges.

Phys. Lett. B , 777:125–130, 2018.[75] Jean Cleymans, Helmut Oeschler, and KrzysztofRedlich. Inﬂuence of impact parameter on thermal de-scription of relativistic heavy ion collisions at (1-2) A-GeV.

Phys. Rev. C , 59:1663, 1999.[76] A. Andronic, P. Braun-Munzinger, and J. Stachel. Ther-mal hadron production in relativistic nuclear collisions:The Hadron mass spectrum, the horn, and the QCDphase transition.

Phys. Lett. B , 673:142–145, 2009. [Er-ratum: Phys.Lett.B 678, 516 (2009)].[77] Anton Andronic, Peter Braun-Munzinger, KrzysztofRedlich, and Johanna Stachel. Decoding the phasestructure of QCD via particle production at high en-ergy.

Nature , 561(7723):321–330, 2018.[78] J. Manninen and F. Becattini. Chemical freeze-out inultra-relativistic heavy ion collisions at s(NN)**(1/2) =130 and 200-GeV.

Phys. Rev. C , 78:054901, 2008.[79] B. I. Abelev et al. Identiﬁed particle production, az-imuthal anisotropy, and interferometry measurementsin Au+Au collisions at s(NN)**(1/2) = 9.2- GeV.

Phys.Rev. C , 81:024911, 2010.[80] M. M. Aggarwal et al. Scaling properties at freeze-out in relativistic heavy ion collisions.

Phys. Rev. C ,83:034910, 2011.[81] Ralf Rapp and Edward V. Shuryak. Resolving the anti-baryon production puzzle in high-energy heavy ion col-lisions.

Phys. Rev. Lett. , 86:2980–2983, 2001.[82] P. Braun-Munzinger, J. Stachel, and Christof Wetterich.Chemical freezeout and the QCD phase transition tem-perature.

Phys. Lett. B , 596:61–69, 2004.[83] J. Noronha-Hostler, C. Greiner, and I. A. Shovkovy.Fast equilibration of hadrfons in an expanding ﬁreball.

Phys. Rev. Lett. , 100:252301, 2008.[84] Anton Andronic, Peter Braun-Munzinger, BengtFriman, Pok Man Lo, Krzysztof Redlich, and JohannaStachel. The thermal proton yield anomaly in Pb-Pbcollisions at the LHC and its resolution.

Phys. Lett. B ,792:304–309, 2019.[85] Jean Cleymans, Pok Man Lo, Krzysztof Redlich,and Natasha Sharma. Multiplicity dependence of(multi)strange baryons in the canonical ensemble withphase shift corrections.

Phys. Rev. C , 103(1):014904,2021.[86] Fernando Antonio Flor, Gabrielle Olinger, and ReneBellwied. Flavour and Energy Dependence of Chemi-cal Freeze-out Temperatures in Relativistic Heavy IonCollisions from RHIC-BES to LHC Energies.

Phys. Lett.B , 814:136098, 2021.[87] Rene Bellwied, Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Jacquelyn Noronha-Hostler, PaoloParotto, Attila Pasztor, Claudia Ratti, and Jamie M.Staﬀord. Oﬀ-diagonal correlators of conserved chargesfrom lattice QCD and how to relate them to experiment.

Phys. Rev. D , 101(3):034506, 2020.[88] J. Adam et al. Net-proton number ﬂuctuations and theQuantum Chromodynamics critical point. 1 2020.[89] Toshihiro Nonaka. Measurement of the Sixth-OrderCumulant of Net-Proton Distributions in Au+Au Col-lisions from the STAR Experiment.

Nucl. Phys. A ,1005:121882, 2021.[90] Ashish Pandav. Measurement of cumulants of con-served charge multiplicity distributions in Au+Au col-lisions from the STAR experiment.

Nucl. Phys. A ,1005:121936, 2021.[91] C. Bernard, T. Burch, E. B. Gregory, D. Toussaint,Carleton E. DeTar, J. Osborn, Steven Gottlieb, U. M.Heller, and R. Sugar. QCD thermodynamics with threeﬂavors of improved staggered quarks.

Phys. Rev. D ,71:034504, 2005.[92] A. Bazavov et al. Fluctuations and Correlations ofnet baryon number, electric charge, and strangeness:A comparison of lattice QCD results with the hadronresonance gas model.

Phys. Rev. D , 86:034509, 2012.[93] Najmul Haque, Aritra Bandyopadhyay, Jens O. Ander-sen, Munshi G. Mustafa, Michael Strickland, and NanSu. Three-loop HTLpt thermodynamics at ﬁnite tem-perature and chemical potential.

JHEP , 05:027, 2014.[94] Colin Morningstar and Mike J. Peardon. Analyticsmearing of SU(3) link variables in lattice QCD.

Phys.Rev. D , 69:054501, 2004.[95] M. Tanabashi et al. Review of Particle Physics.

Phys.Rev. D , 98(3):030001, 2018.[96] Masakiyo Kitazawa and Masayuki Asakawa. Revealingbaryon number ﬂuctuations from proton number ﬂuc-tuations in relativistic heavy ion collisions.

Phys. Rev.C , 85:021901, 2012.[97] V. Skokov, B. Friman, and K. Redlich. Volume Fluctu-ations and Higher Order Cumulants of the Net BaryonNumber.

Phys. Rev. C , 88:034911, 2013.[98] Adam Bzdak, Shinichi Esumi, Volker Koch, JinfengLiao, Mikhail Stephanov, and Nu Xu. Mapping thePhases of Quantum Chromodynamics with Beam En-ergy Scan.

Phys. Rept. , 853:1–87, 2020.[99] Peter Braun-Munzinger, Bengt Friman, KrzysztofRedlich, Anar Rustamov, and Johanna Stachel. Rel-ativistic nuclear collisions: Establishing the non-criticalbaseline for ﬂuctuation measurements. 7 2020.[100] Volodymyr Vovchenko, Oleh Savchuk, Roman V.Poberezhnyuk, Mark I. Gorenstein, and Volker Koch.Connecting ﬂuctuation measurements in heavy-ion col-lisions with the grand-canonical susceptibilities.