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

The hadron resonance gas (HRG) model is often believed to correctly describe the confined phase of QCD. This assumption is the basis of many phenomenological works on QCD thermodynamics and of the analysis of hadron yields in relativistic heavy ion collisions. We use first-principle lattice simulations to calculate corrections to the ideal HRG. Namely, we determine the sub-leading fugacity expansion coefficients of the grand canonical free energy, receiving contributions from processes like kaon-kaon or baryon-baryon scattering. We achieve this goal by performing a two dimensional scan on the imaginary baryon number chemical potential ($\mu_B$) - strangeness chemical potential ($\mu_S$) plane, where the fugacity expansion coefficients become Fourier coefficients. We carry out a continuum limit estimation of these coefficients by performing lattice simulations with temporal extents of $N_\tau=8,10,12$ using the 4stout-improved staggered action. We then use the truncated fugacity expansion to extrapolate ratios of baryon number and strangeness fluctuations and correlations to finite chemical potentials. Evaluating the fugacity expansion along the crossover line, we reproduce the trend seen in the experimental data on net-proton fluctuations by the STAR collaboration.


INTRODUCTION
The study of the QCD phase diagram has been a very active area of research for the last few decades. While much is known about the thermodynamics of QCD at zero baryon number chemical potential, such as the temperature of the crossover transition [1][2][3][4] and the equation of state [5][6][7][8], the properties of the theory at finite baryon densities remain elusive. Effective models predict that the crossover transition turns into a real phase transition at a critical endpoint [9][10][11]. However, confirmation of this feature is needed from a first principles approach and/or experiment. The main goal of the currently ongoing experimental effort at the second Beam * Corresponding author:apasztor@bodri.elte.hu Energy Scan program at RHIC in 2019-2021 is locating the supposed critical endpoint of QCD.
The confined phase of QCD is often assumed to be well described by the ideal hadron resonance gas (HRG) model [47][48][49][50][51]. The HRG model is based on the assumption that a gas of interacting hadrons can be described as a gas of non-interacting hadrons and resonances. The inclusion of the resonances as free particles is an approximate way of taking into account resonant interactions between the stable hadrons [49,50]. The model describes bulk thermodynamic observables -like the pressure or the energy density -obtained from first principle lattice calculations rather well at zero chemical potential [6,[52][53][54]. However, when looking at observables probing finite chemical potentials, namely Taylor expansion coefficients in the chemical potentials µ B , µ S and µ Q near zero or in the respective fugacities e µ B /T , e µ S /T and e µ Q /T near 1, some discrepancies start to emerge between the HRG model and lattice calculations. Some of these discrepancies can be traced back to the fact that some fugacity expansion coefficients with baryon number zero and one are underestimated by the HRG model. Though this is not the only possibility, this type of discrepancy can be interpreted within the bounds of the HRG model itself, and has been used to try to infer the existence of as of yet unobserved hadrons [41,[55][56][57]. The other possibility is that instead of more resonances, a better treatment of resonances is needed, taking into account finite widths and also non-resonant interactions [48][49][50][58][59][60]. Of course, both statements can be true at the same time. For a precise description of the thermodynamics, we most likely need better knowledge of the mass spectrum at higher energies, as well as a more accurate treatment of resonances.
Other discrepancies between the ideal HRG and the lattice are impossible to resolve by supposing the existence of more resonances. These discrepancies can be traced back to the observation that even in the temperature range below the crossover, the HRG fails to describe sub-leading fugacity expansion coefficients.
In principle, the hadron resonance gas can be systematically improved by the S-matrix formulation due to Dashen, Bernstein and Ma [48][49][50], which allows for the calculation of the fugacity expansion coefficients, if enough information is known about the scattering matrix of the hadrons. Applying the S-matrix formalism can lead to a better description of QCD thermodynamics. Ref. [60] shows that the baryon-electric charge correlation χ BQ 11 is particularly sensitive to pion-nucleon scattering phase shifts and that the inclusion of these phase shifts into a hadron gas analysis leads to an improved description of the lattice data. This observable is somewhat special though, in that if isospin symmetry is assumed |S| = 1 hyperons do not contribute at all to χ BQ 11 , and therefore it is only pion-nucleon scattering that dominates the non-resonant contributions.
For other observables more scattering data, e.g. information about baryon-baryon scattering would also be needed. This is especially the case at finite baryon density. Unfortunately, information on these scattering processes is only partially available. While the nucleonnucleon elastic scattering phase shifts are known experimentally [61][62][63], the inelastic part of the S-matrix is not known. Even less is known about scattering between hadrons other than nucleons. Hyperon-nucleon and hyperon-hyperon interactions have been studied in chiral effective theory [64,65]. In the last few years, the analysis of momentum space correlations for hadron pairs measured in pp and p-Pb collisions has also been used to infer properties of hadron-hadron interactions [66][67][68]. There are also some lattice results available for baryonbaryon scattering [69][70][71], but not yet with a continuum extrapolation. While these mentioned research directions show a clear effort from the community to learn about scattering between other hadrons, the preliminary nature of these results makes the use of the S-matrix formalism for the fugacity expansion impractical at the moment.
One simple way to nevertheless go beyond the ideal hadron resonance gas is to use some kind of mean field model for the short range repulsion and the long range attraction between the baryons. Such models were compared to lattice results in Refs [42,[72][73][74]. These works in particular emphasized the importance of the hard core repulsive interactions between hadrons when describing thermodynamics at finite baryon chemical potential. This type of interaction is completely absent in the ideal HRG and leads to a sizable negative contribution to the fugacity expansion coefficients with baryon number two. Such approaches, while interesting, are very far from the first principle approach that the S-matrix formulation could provide, were the necessary S-matrix elements known. In fact, the flavour dependence of excluded volume parameters used in the literature so far have been quite arbitrary, often assuming the same excluded volume for all hadrons. We believe the present calculation of the fugacity expansion coefficients can lead to the construction of more realistic models.
Going beyond equilibrium in the grand canonical ensemble, versions of the hadron resonance gas model have also been used to interpret hadron yields in heavy ion collision experiments. This approach is colloquially referred to as thermal fits as they involve the estimation of the temperature and chemical potential where the yields of hadrons are frozen, the so-called chemical freeze-out conditions. This approach was successful in describing hadron yields [75][76][77][78][79][80], which is quite remarkable, considering that these yields at a single collision energy span many orders of magnitude. Though an important underlying assumption here is the equilibration of the system produced in heavy ion collisions [81][82][83], the fact that the fits work also provides some evidence for this assumption. In this context, it has been realized that including the pion-proton phase shifts in the analysis changes the predicted yields as compared to the ideal HRG at LHC energies [84]. Of course, for consistency, one should extend such and S-matrix treatment to strange hadrons as well. In lack of the necessary scattering data, this extension to strangeness is not straightforward [85]. The inclusion of these corrections is important to precisely test the assumption of a single freeze-out temperature. As a competitor to this assumption, in the context of ideal HRG, it was shown [86] that different freeze-out temperatures for light and strange hadrons, can significantly improve the description of the experimental yields at LHC and the highest RHIC energies.
Since comparisons with the available lattice data suggest that the agreement between full QCD and the ideal hadron resonance gas model gets worse at finite chemical potential, we suspect that non-resonant scattering effects will be even more important at the RHIC Beam Energy Scan and future experiments at lower collision energies, like FAIR and NICA.
In this work we calculate sub-leading fugacity expansion coefficients with first principle lattice simulations. To this end we perform simulations at imaginary chemical potentials, where the fugacity expansion coefficients turn into Fourier coefficients in the imaginary values of the chemical potentials. This correspondence was already exploited in our earlier works. In Ref. [41] we made a detailed analysis of the fugacity expansion coefficients already appearing in the Boltzmann approximation of the ideal HRG, to infer the existence of not yet discovered strange hadrons. In Ref. [42], some of us used the fugacity expansion to emphasize the importance of repulsive baryonic interactions near the crossover region. In Ref. [87] we compared the fugacity expansion with the Taylor expansion in the chemical potentials for cross-correlators of conserved charges. Here we go beyond our earlier works by performing lattice simulations on a 2 dimensional grid in the purely imaginary (µ I B , µ I S ) plane. This allows us for the first time to separate the scattering contributions to QCD thermodynamics by the net strangeness quantum number of the participants. In addition to giving insight on the origin of the discrepancies between full QCD and the ideal HRG model, we believe our results on the fugacity expansion coefficients will also be useful to tune the parameters of the freezeout models in heavy ion phenomenology.
We also use the truncated fugacity expansion to extrapolate experimentally measured ratios of baryon number and strangeness susceptibilities to finite baryon chemical potentials on the phenomenologically relevant strangeness neutral line. This provides an alternative extrapolation procedure to the standard Taylor method. When we extrapolate on the crossover line at strangeness neutrality, these sub-leading coefficients approximately reproduces the trend seen in the experimental data of the STAR collaboration of net-proton fluctuations [88][89][90].
The structure of the paper is as follows. In the next section, we introduce the basic notation and observables used in our study. In Sec. III we discuss our lattice setup. In Sec. IV we discuss our fitting procedure for the sectors and we present the fugacity expansion coefficients.
In Sec.V we calculate the fluctuation ratios using the fugacity expansion and extrapolate to small finite density. Finally in Sec. VI we give a brief summary and outlook for future work.

A. Susceptibilities and the Taylor expansion
There is a conserved charge corresponding to each quark flavor of QCD. Working with three flavors, the grand canonical partition function can be then written in terms of three quark number chemical potentials µ u , µ d and µ s . The generalized susceptibilities are defined to be derivatives of the grand potential (or pressure) with respect to these chemical potentials: with the dimensionless chemical potentialsμ X = µ X /T . For the purpose of hadronic phenomenology it is more convenient to work with the conserved charges B (baryon number), 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 basis of µ B , µ Q , µ S with a simple linear transformation, whose coefficients are given by the B, Q and S charges of the individual quarks: Analogously to the case of the quark number chemical potentials, the susceptibilities are then defined as It is straightforward to express the susceptibilities defined in Eq. (5) in terms of the coefficients in Eq. (1) [24,91,92]. The susceptibilities at µ B = µ S = µ Q = 0 are (up to a trivial factorial factor) the Taylor expansion coefficients of the pressure near that point. Due to charge conjugation symmetry, only the even derivatives contribute. In the present study, we always take µ Q = 0 and only consider derivatives with respect to µ B and µ S . The Taylor expansion therefore reads: where χ BS 00 is just the dimensionless pressure at zero chemical potential.
We note that the Taylor expansion is probably the most natural expansion to work within the plasma phase of QCD. As an exhibit of this, the pressure in the Stefan-Boltzmann (or infinite temperature) limit reads: In this approximation, all derivatives above 4th order are zero, and therefore the Taylor expansion is rapidly convergent. Calculating corrections to this free gas behavior in resummed perturbation theory leads to a non-zero, but small, value for the sixth-order derivatives [93], leaving the qualitative conclusion of the fast convergence of the Taylor series in the plasma phase unchanged.

B. Fugacity expansion of the free energy
An alternative to the Taylor expansion discussed in the previous subsection is a Laurent expansion in the fugacity parameters eμ B and eμ S near 1. Due to charge conjugation symmetry, a combination e mμ B +nμ S and its reciprocal have the same expansion coefficients, making the Laurent expansion an expansion in hyperbolic cosines: The coefficients P BS jk are also called fugacity expansion or sector coefficients, alluding to the fact that they get contributions from the Hilbert subspace corresponding to the fixed values of the conserved charges B = j and S = k. In the ideal HRG model, the expansion coefficients P BS 00 , P BS 01 , P BS 10 , P BS 11 , P BS 12 , P BS 13 all get sizable contributions from known hadrons and hadron resonances. In the Boltzmann approximation to the ideal HRG, coefficients like P BS 20 are zero, while in the full HRG they are non-zero, but very small in magnitude and essentially negligible.
At purely imaginary chemical potentials µ q = iµ I q , where the sign problem is absent and lattice simulations can be performed, we have a Fourier expansion of the form: Differentiation with respect to the original chemical potentials µ B = iµ I B and µ S = iµ I S gives: Im . (14) These formulas and the higher order derivatives of these will be used in our fitting procedure, to be described in Section 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 known hadronic resonances H: with: In the HRG model, the χ BQS ijk susceptibilities of Eq. (5) can be expressed as: where the phase space integral at order i + j + k reads: The fugacity expansion coefficients P BS 00 , P BS 01 , P BS 10 , P BS 11 , P BS 12 and P BS 13 can be obtained via the expansion of equation (16) in terms of the modified Bessel functions K 2 : The Boltzmann approximation consists of taking only the n = 1 term in the above expansion, which accounts for the lowest order in the fugacity parameters. In the Boltzmann approximation, the sectors read: In the full ideal HRG, a hadron with B H = 1 and S H = 0 will also give contributions to the higher order sectors, such as P BS 20 and P BS 30 , due to the terms n = 2 and n = 3 in Eq. (19), respectively. These are, however, exponentially suppressed due to the behavior of the Bessel function K 2 (x) ∼ π 2x e −x as x → ∞. These contributions are orders of magnitude smaller than the full weight of the respective sectors as obtained from the lattice.
The HRG model is an approximation of the more general formula by Dashen, Bernstein and Ma, which gives the fugacity expansion coefficients in terms of the Smatrix: where M BS jk is the mass threshold for the B = j, S = k channel, the trace is taken over this Hilbert subspace, and the subscript c signifies that only connected S-matrix elements are to be taken. For the specific case of elastic 2 → 2 body scattering, where the δ J,I are the scattering phase shifts for angular momentum J and isospin I and the isospin singlet and triplet contributions have been written separately. After integration by parts with respect to E, we get to the conclusion that the contribution of elastic scattering is given by the integral of the phase shift with an exponential weight. This leads to the expectation that dominantly repulsive interactions will lead to a negative sub-leading fugacity expansion coefficient. This fact was exploited when constructing repulsive core HRG models and comparing them with lattice data in Refs. [42,74]. It is also reasonable to expect, due to the exponential suppression of the K 2 Bessel functions, that in the hadronic phase there will be a strong hierarchy of the fugacity expansion coefficients with increasing quantum numbers, so e.g. P BS It is thus a reasonable expectation that, in the hadronic phase, the fugacity expansion will converge faster than the Taylor expansion. This is the opposite situation as in the plasma phase, where the Taylor expansion converges quickly, while the fugacity expansion converges slowly. This makes the fugacity expansion particularly useful for modelling the hadronic phase, and therefore also for the study of chemical freezeout in heavy ion collisions.
Here we do not utilize any S-matrix formula or any mean field approximation thereof, but rather calculate the sub-leading sector coefficients P BS 20 , P BS 21 , P BS 22 , P BS 02 etc. directly from lattice simulations.

III. LATTICE SETUP
We use a staggered fermion action with 4 steps of stout smearing [94] with the smearing parameter ρ = 0.125 and a tree-level Symanzik-improved gauge action. This combination was first used in Ref. [24], where information about the line of constant physics can be found. For the scale setting we use the pion decay constant f π = 130.41MeV [95]. We use lattices of temporal extent N τ = 8, 10 and 12 to perform an estimation of the continuum value of our observables. The spatial extent of the lattice is given by the aspect ratio LT ≈ 3. Due to technical reasons, some lattices had slightly different values for this ratio, as given in Table. I. Given the error bars on the final results, we did not optimize this further.
For the continuum extrapolations, we assume a linear scaling in 1/N 2 τ . Since taste breaking effects are still rather large on these lattices, we only call our results continuum estimates, as opposed to fully controlled continuum extrapolations, e.g. when N τ = 16 is part of the extrapolation.
For all values of N τ , we use simulations at four different temperatures T = 145 MeV, 150 MeV, 155 MeV and 160 MeV. At each temperature and each lattice spacing, we perform a two-dimensional scan in the imaginary chemical potentials µ I B and µ I S , with the chemical potentials taking the values (µ I B , µ I S ) = π 8 (i, j), with i = 0, 1, . . . , 15 and j = 0, 1, . . . , 9, for a total of 9 × 16 = 144 simulation points. In each µ i = 0 point, we simulated one Rational Hybrid Monte Carlo stream with several thousand trajectories, evaluating every fifth configuration for the fluctuation observables as detailed in Ref. [24]. Our statistics is summarized in Table I.
The statistical errors are calculated using the jackknife method. The estimation of the systematic errors is a more elaborate process. Ambiguities appear at various points of the analysis, e.g. in the way the continuum extrapolation is calculated, or how many fit parameters we use for the extraction of the fugacity expansion coefficients. We consider all combinations of the possibilities and take the spread of the results as systematic error.

IV. FUGACITY EXPANSION COEFFICIENTS
The estimation of the coefficients P BS ij proceeds through a correlated fit. On the µ B = µ S = 0 ensembles, the fluctuations χ BS   Only the leftmost five are accounted for by ideal HRG. The next seven appear in the next order, which we use later for phenomenology. For the next order (last four coefficients with B = 3) we see no stable signal. There are two symbols per coefficient, triangles for the complete fit and circles for the case without the B = 3 part. The data refer to our finest lattice spacing. Right: Example for the combined continuum extrapolation of the extracted coefficients.
use Im χ BS 10 and Im χ BS 01 . This leads to a block-diagonal covariance matrix with one 8 × 8 block corresponding to the µ = 0 ensemble, and 143 blocks of size 2 × 2, corresponding to the ensembles with a non-zero value of at least one of the chemical potentials. The covariance matrix blocks are estimated by the jackknife method with 24 jackknife samples. The truncation of the fugacity expansion is somewhat ambiguous, as there is no single small parameter in which we actually perform this expansion. To estimate systematic errors coming from the choice of the ansatz, we therefore perform two fits for each ensemble, for which we introduce the shorthand notations B max = 2 and B max = 3. The sectors included in the B max = 2 analysis are: The first five of these correspond to sectors that are already present in the ideal HRG in the Boltzmann approximation. They also set contributions from interactions though, e.g. non-resonant pion-nucleon interactions contribute to P BS 10 , while K-Λ interactions contribute to P BS 12 . The P BS 1,−1 is an exotic (pentaquark) sector. It gets contributions for example for valence quark content uudds which corresponds e.g. to p + K 0 scattering. The sectors P BS 2i get contributions from baryon-baryon scattering: P BS 20 from N − N , P BS 21 from N − Λ, P BS 22 from N − Ξ or Λ − Λ and finally P BS 23 from N − Ω or Λ − Ξ. In each case, Σ can replace Λ. The coefficients P BS 02 and P BS 03 get contributions from two-and three-kaon scattering, respectively. The inclusion of the P BS 03 sector with the omission of the P BS 30 sector is motivated by the lower mass threshold of 3 kaon scattering as compared to 3 baryon scattering. In addition to these, we also performed an analysis were four sectors with B = 3 were added: These get contributions from three-baryon scattering, with various strangeness contents. Our data is not yet sufficiently accurate to obtain a reliable estimate of the sectors with B = 3, and the inclusion of these sectors does not improve the χ 2 of the fits. Whether we include these, or not, the results for the B = 2 sectors remain consistent, as we show in the left panel of Fig. 1, where the sector coefficients from the two different fits on the N τ = 12 lattices are shown. This way we demonstrate the stability of the sectors included in the B max = 2 set. Only at the highest temperature T = 160MeV, and only for one sector, P BS 20 , is the systematic error coming from including the B = 3 sectors comparable to the statistical error of the fits.
The continuum limit estimation of the sectors proceeds through a combined fit in temperature and lattice spacing (or equivalently N τ ) via the ansatz (25) For the systematic error we compare this ansatz with and without the coefficient b 2 . The continuum extrapolation of the beyond-ideal-HRG sectors for the case of the B max = 2 fits at fixed T and N τ and b 2 kept as a free parameter, is shown in Fig 1 (right). The other 3 fits look quantitatively similar. All of the continuum fits have acceptable fit quality, with Q values over 1%. As a conservative estimate of the systematics, we combine them with uniform weights. As can be seen in the right panel of Fig 1, the slopes of the continuum extrapolations of all beyond-ideal-HRG sectors appear to be mild, except for the sector P BS 02 , which corresponds to kaon-kaon scattering, and changes its sign during the continuum extrapolation. As expected, this sector -being related to kaons -suffers from relatively large taste-breaking effects.
The final results for the beyond-ideal-HRG sectors can be seen in Fig. 2. Within the statistical precision of our results, P BS 20 is roughly the same as P BS 21 , while P BS 22 is smaller than the previous two. As a comparison, the ideal HRG model prediction for the sum k P BS 2k at T = 155MeV is of the order 10 −5 , orders of magnitude lower than what we see here. The two-kaon scattering sector P BS 02 goes slightly below zero at around 155MeV within 1σ uncertainty. The three-kaon sector P BS 03 is consistent FIG. 2. Our continuum estimates of the beyond-ideal-HRG sector coefficients. We show the results both from our combined temperature and continuum fit (green bands) and of a T -by-T continuum limit extrapolation (blue points). Systematic errors are included, in the first case by varying Bmax = 2 vs 3 and b2 = 0 or b2 = 0, and in the second case by varying Bmax = 2 vs 3.
with zero in the entire temperature range and is therefore not included in the plot. An upper limit on its magnitude with 1σ uncertainty is 2 · 10 −3 . The exotic sector P BS 1,−1 is rather large, consistently with our earlier, statistically independent finding in Ref. [87] on N τ = 12 lattices. We have already published the leading sector coefficients for which the ideal HRG has a prediction in the Boltzmannapproximation, namely P BS 01 , P BS 10 , P BS 11 , P BS 12 , P BS 13 in Ref. [41]. We will not repeat the results for those sectors here.

V. FLUCTUATION-RATIOS AT FINITE BARYON DENSITY
Having the coefficients of the fugacity expansion, the thermodynamics can be readily obtained. In this section we calculate the baryon number and strangeness fluctuations and their ratios in the studied temperature range. There is no difficulty in evaluating Eq. 8 and its µ B -and µ S -derivatives at any chemical potential.
Heavy ion collisions involving lead or gold atoms correspond to the conditions χ S 1 = 0 and χ B 1 = 0.4χ Q 1 . For the purposes of the present study, we impose strangeness neutrality and leave the second conditions for future work. In fact, we use the simplified form χ B 1 = 0.5χ Q 1 , which is realized at vanishing electric charge chemical potential.
To show the magnitude of the cut-off effects, we start with a pair of quantities at µ B = 0. The fugacity expansion, and more generally imaginary chemical potential simulations offer an efficient way to calculate susceptibilities at µ B = 0. In Ref. [44] we calculated the ratios χ B 3 /χ B 1 and χ B 4 /χ B 2 at a finite lattice spacing of N τ = 12 with the same lattice action used here. Here we show continuum estimates of the fluctuation ratios χ B 3 /χ B 1 and χ B 4 /χ B 2 at µ B = 0, together with the data at finite N τ in Fig. 3. In the ideal HRG model χ B 4 /χ B 2 = 1 for all temperatures, meaning that above T = 150MeV our results show a clear deviation from the HRG prediction, due to presence of the non-zero beyond-ideal-HRG sectors. The difference between the two ratios is also shown. In the µ S = 0 case the two ratios at µ B = 0 are identical. The difference between the two ratios comes from imposing the strangeness neutrality condition χ S 1 = 0. This difference also shows mild cut-off effects.
After the sectors are obtained, we perform extrapolations to real chemical potentials using the ansatz of Eq. (8) truncated at the B max = 2 level. We extrapolate first at fixed T and N τ . We consider the fluctuation ratios χ B 1 /χ B 2 , χ B 3 /χ B 1 , χ B 4 /χ B 2 and χ BS 11 /χ S 2 on the strangeness-neutral line χ S 1 = 0, which determines µ S as a function of µ B . While the extrapolation always uses the 12 sectors of the B max = 2 level, the values of these sectors are taken both from the B max = 2 and B max = 3 fits, to estimate the systematic errors. We then perform a continuum estimation at fixed values of µ B /T with the same combined T and N τ fit as in the case of the baryon and strangeness sectors.
We had one sector, P BS 02 , with a steep continuum extrapolation. Should we expect additional systematic errors coming from the non-trivial continuum scaling in the phenomenology? The answer is no, as we demonstrate in Fig. 4. The multi-kaon sectors do not contribute to the baryon fluctuations. The only ratio with phenomenological relevance where the P BS 02 may be important is χ BS 11 /χ S 2 . We calculated this ratio with and without the multi-kaon sectors and compared the results in Fig. 4 at T = 160 MeV. Although at this temperature and this FIG. 3. Fluctuation-ratios χ B 3 /χ B 1 and χ B 4 /χ B 2 obtained from the fugacity expansion truncated at the Bmax = 2 level at µB = 0 on our Nτ = 8, 10 and 12 lattices and our continuum estimates from these data. The points at finite lattice spacing include a systematic error coming from whether we used the 12-or the 16-parameter fit to determine the Bmax = 2 sectors. The continuum results include systematic error from 4 fits, in addition for the 12 vs 16 parameter fits at fixed Nτ and T we also include a 5-vs 6-parameter combined T and Nτ continuum fit. lattice spacing we have the largest value for this difficult sector, we see hardly any significant effect from it, especially not when compared to the statistical errors after the continuum extrapolation step.
The final results for the fluctuation ratios on the strangeness neutral line can be seen in Fig. 5. The first of these ratio is strongly dependent on the chemical potential, but not on the temperature, making it a proxy of the chemical potential, at least for small values of µ B . On has to remember though that if the critical endpoint exists, the fluctuation χ B 2 diverges there, leading to χ B 1 /χ B 2 → 0 at the critical point, and therefore making this quantity a nonmonotonic function of µ B .
The other three are more strongly dependent on the temperature and less strongly on the chemical potential, therefore making them possible proxies for the temperature. The ratios χ B 3 /χ B 1 and χ B 4 /χ B 2 can be regarded as a baryon thermometer, while the ratio χ BS 11 /χ S 2 as a strangeness-related one. This latter ratio is of large phenomenological interest, as experimental net-lambda and net-kaon fluctuations can be used to construct the ratio σ 2 Λ /(σ 2 Λ + σ 2 K ). It was shown in Ref. [87] that this is a good experimental proxy of χ BS 11 /χ S 2 , not strongly affected by experimental effects, which makes it a prime target for comparison with experiments.
We show the fluctuation ratios in Fig. 5 as functions of the dimensionless chemical potential µ B /T at a few values of the temperature, as well as on the crossover line calculated to order µ 2 B in Ref. [45]: with T 0 c = 158.0 ± 0.6 and κ 2 = 0.0153 ± 0.0018. The errors on these numbers are included in the error estimation, but are negligible. Note that, since the crossover temperature changes very little in the chemical potential range of our study, the 1σ bands on the T c (µ B ) line always overlap with the 1σ bands for a fixed T = T 0 c = 158MeV.
Our results on the ratios χ B 3 /χ B 1 and χ B 4 /χ B 2 are also shown as a function of χ B 1 /χ B 2 in Fig. 6. Within 1σ, our results are consistent with recent lattice results on the susceptibility ratios using the Taylor method [30]. In Fig. 6 we also compare to the data on net-proton fluctuations from the STAR collaboration [88][89][90], the net-proton skewness-to-mean ratio C 3 /C 1 and the netproton kurtosis-to-variance ratio C 4 /C 2 , as functions of the net-proton mean-to-variance ratio C 1 /C 2 at chemical freeze-out. The advantage of using these variables in the comparison is that is does not involve any modelling

FIG. 5. Continuum estimates of the fluctuation ratios
from a fugacity expansion truncated at the Bmax = 2 level, shown as a function of the dimensionless chemical potential µB/T for fixed temperatures, as well as on the crossover line Tc(µB).
of the freeze-out conditions, other than assuming that chemical freeze-out happens close to the QCD crossover on the phase diagram. Our results are consistent with the experimental results. While such a direct comparison suffers from many caveats [87,[96][97][98][99][100], the similarity in the trends supports the idea that experimentally observed net-proton fluctuation ratios reflect with some accuracy the thermal fluctuations in an equilibrated QCD medium.

VI. SUMMARY AND OUTLOOK
We have calculated fugacity expansion coefficients of the QCD pressure beyond the ideal HRG model, separating contributions to the QCD free energy coming Our continuum estimates of the fluctuation ratios χ B 3 /χ B 1 and χ B 4 /χ B 2 compared with STAR data on net-proton fluctuations from Ref. [88]. Both the values of the fluctuations and the trend as a function of baryon density are consistent.
from Hilbert subspaces with different values of the baryon number and strangeness quantum numbers. This allows one to quantify the importance of processes like kaonkaon and baryon-baryon scattering, when modelling the QCD medium in the hadronic phase, but close to the crossover. We estimated the continuum value of these coefficients with lattice simulations of temporal extent N τ = 8, 10 and 12 using the staggered discretization. We observed large cut-off effects in the kaon and multi-kaon sectors, only. To study these, and to make the continuum extrapolation more robust, the inclusion of finer lattices is desirable.
Note that our study was limited to an aspect ratio of LT ≈ 3, future studies should also investigate finite volume effects in the baryon-strangeness sectors. However, a strong volume dependence is more likely to be observed in correspondence with the electric charge sectors. The full picture of baryon interactions in the hadronic phase will emerge from the three-dimensional mapping of the µ B − µ S − µ Q space. In this work we restricted the space to µ Q = 0. The sectors we obtained are sums of various charge sectors, P BS ij = k P BSQ ijk , and we cannot differentiate between the terms, though it would be possible in a more elaborate setup. Still, the level of separation achieved in this work already provides plenty of new information for hadronic modelling of the QCD medium.
We used the truncated fugacity expansion to calculate phenomenologically relevant fluctuation ratios on the strangeness neutral line both as a function of the chemical potential and as a function of the baryon number meanto-variance ratio, which can be regarded as a proxy of the baryo-chemical potential. While a direct comparison is by no means trivial, the fugacity expansion coefficients appear to describe the trend in the STAR data on netproton fluctuations [88][89][90].
It has been pointed out in the literature, that the modifications of the ideal HRG model that include the effects of global baryon number conservation lead to a reduction in higher order fluctuations and therefore describe the experimental data better. For a recent study of these effects see [99]. In a recent paper, it was also pointed out that the decrease in the kurtosis with increasing chemical potential observed in the data is likely not due to critical point effects [101]. The corrections to ideal HRG studied in this work have a similar magnitude as the experimental effects like canonical effects and volume fluctuations. Since both baryon interactions and the global conservation laws appear to push the χ B 4 /χ B 2 and χ B 3 /χ B 1 ratios down, it is important to have an estimate of the relative size of these types of effects under realistic conditions. Performing a study of this kind is an important task for the near future, as it will guide the correct interpretation of STAR data.

ACKNOWLEDGMENTS
This project was partly funded by the DFG grant SFB/TR55 and also supported by the Hungarian National Research, Development and Innovation Office, NKFIH grants KKP126769. The project leading to this publication has received funding fromExcellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme, AMX-18-ACE-005. The project also received support from the BMBF grant 05P18PXFCA. Parts of this work were supported by the National Science Foundation under grant no.