Second order cumulants of conserved charge fluctuations revisited I. Vanishing chemical potentials

We update lattice QCD results for second order cumulants of conserved charge fluctuations and correlations at non-zero temperature and vanishing values of the conserved charge chemical potentials. We compare these results to hadron resonance gas calculations with and without excluded volume terms as well as S-matrix results in the hadronic phase of QCD, and comment on their current limitations. We, furthermore, use these results to characterize thermal conditions in the vicinity of the pseudo-critical line of the chiral transition in QCD. We argue that the ratio of strange to baryon chemical potentials is a robust observable that, on the one hand, deviates only little from hadron resonance gas results, but, on the other hand, is very sensitive to the spectrum of strange baryon resonances.

We update lattice QCD results for second order cumulants of conserved charge fluctuations and correlations at non-zero temperature and vanishing values of the conserved charge chemical potentials. We compare these results to hadron resonance gas calculations with and without excluded volume terms as well as S-matrix results in the hadronic phase of QCD, and comment on their current limitations. We, furthermore, use these results to characterize thermal conditions in the vicinity of the pseudo-critical line of the chiral transition in QCD. We argue that the ratio of strange to baryon chemical potentials is a robust observable that, on the one hand, deviates only little from hadron resonance gas results, but, on the other hand, is very sensitive to the spectrum of strange baryon resonances.

I. INTRODUCTION
The thermal state of strong-interaction matter described by Quantum Chromodynamics (QCD) with two light up and down quarks and a heavier strange quark is fixed by four external control parameters, the temperature (T ) and three chemical potentials (µ B , µ Q , µ S ) that couple to the conserved currents for net baryonnumber (B), electric-charge (Q) and strangeness (S), respectively. In a grand canonical description of equilibrium thermodynamics these external control parameters are Lagrange multiplier, appearing in the statistical operator, that characterize the thermodynamic conditions of, e.g., matter described by the strong force. The set of external parameters, (T, µ B , µ Q , µ S ) may be used to calculate various bulk thermodynamic observables, e.g. energy and number densities as well as generalized susceptibilities that are obtained from higher order derivatives of the grand canonical partition function. Such observables are, in principle, measurable in experiments, while the set of Lagrange multiplier, (T, µ B , µ Q , µ S ), is not. They are specific to a given model to the extent that the value of a certain bulk thermodynamic observable, e.g. the net baryon-number density, which is related to a certain temperature value T in QCD (or nature) will correspond to another temperature value in a model calculation. The latter will be close to the former only when the model approximations provide a realistic description of QCD.
Lattice QCD calculations provide continuum extrapolated results for bulk thermodynamic observables at small values of the chemical potential with percent-level accuracy. This provides, in particular, the pseudo-critical temperature for the chiral transition at vanishing values of the baryon chemical potential with better than 1% accuracy [1] for a specific set of observables. This result differs by less than 2% from calculations using different observables and discretization schemes to char-acterize the crossover transition in QCD [1,2]. Similarly the temperature dependence of the pseudo-critical line, T pc (µ B ), is known with better than 2% accuracy at least for baryon chemical potentials µ B ≤ T pc,0 ≡ T pc (0) [1][2][3]. This also includes a determination of constraints on the strangeness chemical potential, µ S (µ B ), required to insure strangeness neutrality in strongly interacting matter. At fixed value of temperature this constraint is known to better than 5% for a range of baryon chemical potentials µ B ≤ T pc (µ B ).
In heavy-ion experiments at the Large Hadron Collider (LHC) and the Relativistic Heavy Ion Collider (RHIC) the thermal properties of strongly interacting matter as it existed at the time of hadronization (freeze-out of various hadrons) is probed. The underlying thermal parameters, (T, µ B , µ Q , µ S ), are extracted from particle yields or higher order cumulants characterizing the distributions of these hadrons. These parameters, however, generally are not obtained through a direct comparison with QCD but with statistical models, e.g. models that use the statistical operator of hadron resonance gas (HRG) models [4] for point-like non-interacting hadrons, which are variants of the Hagedorn resonance gas model [5][6][7]. This too leads to a determination of thermal parameters with errors that are on the percent level. The approach, however, may suffer from systematic uncertainties, as the statistical operator used in model calculations obviously will differ from that of QCD in certain ranges of the temperature and baryon chemical potential. Differences between the statistical operator of QCD and that of HRG models based on a spectrum of point-like, non-interacting hadron resonances become apparent in properties of higher order cumulants, e.g. deviations from (generalized) Skellam distributions that are inherent to HRG models using point-like, non-interacting hadron resonances. A proper inclusion of interactions in a hadronic medium, using a relativistic virial expansion [8][9][10][11][12] or more phenomenol-ogy oriented excluded volume [13][14][15][16][17][18] as well as repulsive mean approaches [19][20][21] thus may be necessary. These approaches generally need to adjust parameters which is done by comparing to lattice QCD calculations of thermodynamic observables.
Second order cumulants can be used as important benchmark observables that allow to establish the range of validity of hadronic models that are needed to provide an interface between experimental observables, e.g. hadron yields and fluctuations, and thermal observables obtained in field theoretic calculations, e.g. QCD. In this paper we will provide high precision, continuum extrapolated lattice QCD results for second order cumulants of conserved charge fluctuations and their correlations at vanishing values of the baryon chemical potential. We will extend these calculations in a follow-up paper to nonvanishing values of the chemical potential.
The calculations presented here at vanishing values of the chemical potential allow to constrain the basic parameters that enter in model calculations, e.g. the excluded volume parameters introduced in hadron resonance gas (HRG) models to mimic the effect of repulsive interactions or the modeling of higher order corrections in S-matrix calculations that are needed to go beyond the calculation of second order coefficients in a virial expansion. Preliminary results from this study have been presented by us previously [22,23]. A similar comparison of lattice QCD results on cumulants of conserved charge fluctuations with excluded volume HRG models also appeared recently [24]. This paper is organized as follows. In the next section we describe the computational set-up for our calculations and describe the scale setting used by us to define temperature scales. Section III gives a short description of the basic observables that we analyze. Section IV describes our determination of continuum extrapolated results for all second order cumulants. Section V is devoted to a comparison of the QCD results for these cumulants with hadron resonance gas model calculations. Finally we give our conclusions in Section VI. Three appendices are devoted to the presentation of further details on the scale setting, the fits to second order cumulants and a description of differences between two lists of hadron resonances used in HRG model calculations.

II. COMPUTATIONAL SET-UP
The framework for our calculations with the Highly Improved Staggered Quark (HISQ) [25] discretization scheme for (2 + 1)-flavor QCD with a physical strange quark mass and two degenerate, physical light quark masses is well established and has been used by us for several studies of higher order cumulants of conserved charge fluctuations and correlations. The specific set-up used in our current study has been described in [26].
In addition to the data-sets used for that study we increased the statistics on lattices with temporal extent N τ = 12 and 16. This allows us to present continuum extrapolated results for all second order cumulants, for which the statistical error on the data themselves amounts to less than 50%, the reminder coming from systematic errors related to the fit ansätze used for continuum extrapolations and uncertainties on the zero temperature observables used to set the scale for the temperature.

A. Data sets and statistics
We follow here the notation and conventions used in [27] for the calculation of the equation of state in (2 + 1)flavor QCD with non-vanishing chemical potentials. In that work first continuum extrapolated results for second order cumulants, obtained with the HISQ action, had been presented. The statistics collected for our current analysis is more than a factor 10 larger than that used in [27] for the calculation of the equation of state on lattices with temporal extent N τ = 8 and 12. Moreover, compared to the statistics used previously for a determination of the pseudo-critical temperature T pc of the chiral transition in (2 + 1)-flavor QCD [1], we increased the statistics by a factor (3)(4) for existing data-sets on lattices with temporal extent N τ = 16 in the temperature range T ∈ [135 MeV : 178 MeV]. We also make use of data-sets obtained in earlier calculations on lattices with temporal extent N τ = 6 [1,27]. Details on the simulation parameters and our current statistics are given in Table I. The gauge field configurations stored in our data-sets are separated by 10 time units in a simulation with the rational Hybrid Monte Carlo (RHMC) algorithm.
All our calculations have been performed on lattices with spatial extent N σ = 4N τ . Additional data on second order cumulants on lattices with temporal extent N τ = 6 are taken from [27]. In order to stabilize the asymptotic behavior of our fits to second order cumulants we also use data at higher temperatures, which also are taken from [27].
In this new analysis we focus on a temperature range of 20 MeV above and below the pseudo-critical temperature for the chiral transition at vanishing value of the chemical potentials, T pc,0 = (156.5 ± 1.5) MeV [1].

B. Scale setting
In order to set the scale for the temperature used in lattice QCD calculations we closely follow the strategy laid down for the definition of a line of constant physics and the calculation of the equation of state in [28]. There we introduced temperature scales based on a calculation of r 1 , which characterizes the short distance part of the heavy quark potential, and the kaon decay constant f K . We use parametrizations of both observables to set the scale for the temperature, The two temperature scales are related through the value of r 1 f K , which is solely determined through a lattice calculation, and the physical value of r 1 , which requires input from experiment, e.g. the pion decay constant f π . At non-zero lattice spacing the temperature scales derived from different observables differ. At finite values of the gauge coupling β they are related through The parametrizations of a/r 1 and af K used by us at nonvanishing lattice spacing are given in Appendix A. In all figures, that show lattice QCD results obtained at nonvanishing lattice spacing, we use, for definiteness, the temperature scale based on calculations of af K , as has been done by us also in the past [29]. For these figures and the continuum extrapolations at fixed temperature values we use as basic input the central value of the MILC results for r 1 , i.e. r 1 = 0.3106 fm. The fits presented in Appendix A then fix the central value of the kaon decay constant to f K = 155.7/ √ 2 MeV. As discussed in more detail in Section IV A, we treat the error on r 1 and f K as an overall systematic error on the temperature scale that will enter the final error budget in our analysis.

III. FIRST AND SECOND ORDER MOMENTS OF CONSERVED NET CHARGE FLUCTUATIONS
We focus here on a discussion of second order cumulants of conserved charge fluctuations at vanishing chemical potentials for the conserved charges of (2 + 1)-flavor QCD, i.e. net baryon-number (B), electric charge (Q) and strangeness (S). They are obtained from the QCD partition function, Z(T, V, µ), as derivatives with respect to the associated chemical potentials µ = (µ B , µ Q , µ S ), This set of six second order cumulants are leading order terms in Taylor expansion of various thermodynamic quantities derived from Taylor expansions of the pressure of (2 + 1)-flavor QCD, withμ X ≡ µ X /T and arbitrary natural numbers i, j, k. In particular, they provide the leading order expansion coefficients of mean (n X /T 3 ≡ χ X 1 ) and variance (σ 2 X /T 2 ≡ χ X 2 ) of conserved charge distributions. At leading order in the chemical potentials the former are given by Of particular interest, for a discussion of properties of strongly interacting matter, created in heavy ion collisions, is the case of strangeness neutral matter (n S ≡ 0). The second order cumulants provide important information on the relation between strangeness and baryon chemical potentials in strongly interacting matter [30].
To leading order this is given by, with

11
(10) Here r = n Q /n S . Note that in the isospin symmetric case, r = 1/2, one has q 1 = 0. We will discuss Section V A 4 the sensitivity of this ratio of chemical potentials to the hadron resonance spectrum contributing to the thermodynamics of strongly interacting matter.

IV. SECOND ORDER CUMULANTS AT VANISHING CHEMICAL POTENTIAL
A. Continuum extrapolation of second order cumulants Second order cumulants of conserved charge fluctuations and cross-correlations among them have been calculated in lattice QCD using various discretization schemes. This led to continuum extrapolated results in (2 + 1)flavor QCD with physical (degenerate) light quark masses (m u = m d ) and a physical strange quark mass (m s ) that have been presented, previously [31,32]. The HotQCD Collaboration presented such results obtained in calculations with the HISQ action. Making use of recently obtained high statistics data we update here the continuum extrapolation of the second order cumulants and present fits to a set of four independent second order cumulants that are suitable for analyzing all six 2 nd order cumulants in the conserved charge (B, Q, S) basis as well as the (u, d, s) flavor basis [27] for the chemical potentials.
Among the six second order cumulants only four are independent at vanishing values of the chemical potential. The other two are constrained by isospin symmetry as our calculations, as well as most other lattice QCD calculations, are performed for two degenerate light up and down quarks. This leads to two constraints [27], In the quark flavor basis the four independent observables and the two constraints relate to the three diagonal and three off-diagonal cumulants for the light and strange quark number fluctuations and correlations. The two constraints merely reflect that u and d-quark fluctuations are identical as are their correlations with the strange quarks, i.e. χ u 2 = χ d 2 and χ us 11 = χ ds 11 . We also note that the above constraints are fulfilled to better than 1% also in HRG model calculations that utilize the experimentally determined, physical hadron resonance spectrum [33]. Here small mass differences among isospin partners arise from differences in the up and down quark masses as well as electromagnetic interactions. The above constraints, of course, also hold for HRG model calculations using hadron spectra calculated e.g. within relativistic quark models [34][35][36] as well as for spectra obtained in lattice QCD calculations [37,38].
All possible second order cumulants and ratios thereof, calculated either in the (B, Q, S) or (u, d, s) basis for the chemical potentials [39], thus can be obtained from four independent observables; we focus on the set of observables χ Q 2 , χ QS 11 , χ BQ 11 , χ BS 11 . While the first two cumulants are dominated by contributions from the nonstrange (χ Q 2 ) and strange (χ QS 11 ) meson spectrum the latter two are dominated by the non-strange (χ BQ 11 ) and strange (χ BS 11 ) baryon sectors of the hadron spectrum. The second order cumulants have been obtained in recent high statistics calculations on lattices with temporal extent N τ = 8, 12, 16 that focused on a calculation of higher order cumulants [26,40]. These calculations have been extended, focusing on the temperature range in the vicinity of the pseudo-critical temperature, T pc,0 , for the chiral transition. The current statistics, accumulated in these calculations, is given in Table I. Moreover, we used data sets obtained earlier on lattices with temporal extent N τ = 6 [1,27]. For each N τ we used cubic spline interpolations of our data in the interval T ∈ [134 MeV : 178 MeV] (see Appendix B). These interpolations allow us to obtain for each of the four different lattice sizes cumulants at the same temperature. This has been done for each of the two schemes used to define a temperature scale at non-vanishing values of the lattice spacing.
Continuum extrapolations of these observables have been performed in both schemes, using linear and quadratic extrapolations in 1/N 2 τ , For some temperature values the linear and quadratic extrapolations are shown in Appendix B. Results obtained with linear extrapolations on the N τ = 8, 12, 16 data sets based on the r 1 and f K temperature scales are compared in Fig. 1. Differences in the continuum extrapolation shown in Fig. 1, that arise from the usage of two different temperature scales, reflect systematic errors arising from the parametrization of a/r 1 and af K at finite values of the gauge coupling. Moreover, the restriction to linear fit ansätze reacts differently to truncated higher order corrections in both schemes. Our final continuum extrapolation result for second order cumulants is obtained by averaging over the two different linear fit results. The difference in the two extrapolations at each temperature is taken as a systematic error and has been added linearly to the statistical errors of the two extrapolations based on the r 1 and f K temperature scales, respectively. Further details on our continuum extrapolations are given in Appendix B. (bottom, right) at several values of the temperature, which has been obtained from afK (colored lines) and a/r1 (black lines), respectively. Extrapolations are linear in 1/N 2 τ and for data obtained on lattices with temporal extent Nτ = 8, 12, 16. Crosses indicate the corresponding QMHRG2020 value at that temperature. Note that for χ BQ 11 and χ Q 2 these HRG results are not shown for all temperature values as the deviations from the corresponding QCD results are too large. For χ Q 2 QMHRG2020 with the finite-volume corrected contributions for pions and kaons in a volume LT ≡ Nσ/Nτ = 4 has been used (see Section V B).
The presentation of the continuum extrapolations in Fig. 1 makes use of a specific choice for the value of r 1 to define the temperature scale, i.e. we used 1 r 1 = 0.3106 fm, which is the central value for r 1 quoted by the MILC Collaboration [41], r 1 = 0.3106(8)(14)(4) fm.
Here the errors are statistical, systematic and experimental, respectively. The statistical error and part of the systematic error on the value for r 1 , quoted by the MILC collaboration, corresponds to errors also arising in our calculation when analyzing results obtained with two different T -scales, different fit ansätze and fit ranges used for the continuum extrapolations.
In addition to this error we obviously need to add the systematic uncertainty in the T -scale arising from the error on the experimental value of the pion decay constant, f π , quoted by MILC as ∆r 1 = 0.0004 fm. This gives rise to an uncertainty of the T -scale of about 0.2 MeV in the T -range considered by us. There also is a substantial spread in values for r 1 obtained by several groups and quoted by FLAG [42], for instance in Ta-ble 56. We estimate this overall systematic uncertainty on r 1 as ∆r 1 = 0.001 fm, which for the temperature range considered by us amounts to a scale uncertainty of ∆T 0.6 MeV. The resulting systematic error is shown separately for all our continuum extrapolated results on second order cumulants.
Continuum extrapolated results for all four independent second order cumulants, calculated on lattices with temporal extent N τ = 6, 8, 12 and 16 are shown in Fig. 2 in the temperature range 130 MeV < T < 180 MeV. In that figure we also show separately the systematic error arising from the scale uncertainties discussed in the previous paragraph (red band), and the combined statistical and systematic error arising from the continuum extrapolation using two different T -scales (grey band). The insets in these figures show a comparison between the lattice QCD results and a specific HRG model calculation that is based on the QMHRG2020 list of resonances (see Section V A 1). This will be discussed further in the next section.
The uncertainty of the T -scale determination has been propagated into a systematic error for our continuum extrapolated observables at temperature T . This is quoted as a systematic error on the observables at temperature  T . Results for some values of the temperature are given in Table II and compared with corresponding results obtained in Ref. [32]. As can be seen agreement between these two analyses is quite good. Further details on our continuum extrapolated results and comparisons with various HRG model calculations as well as results obtained from virial expansions will be presented in the next sections.

B. Parametrization of the continuum extrapolated second order cumulants
When performing the continuum extrapolations of second order cumulants we did not assume any specific ansatz for the temperature dependence of the cumulants. It, however, will be useful also for comparisons with other observables and model calculations to have at hand a parametrization of the four independent second order cumulants as well as the two dependent observables, χ B 2 and χ S 2 . For this purpose we provide a parametrization in terms of rational polynomials, where X, Y ∈ B, Q, S, and it is understood that χ XY 11 should be replaced by χ X 2 for X = Y . This parametrization corresponds to the center of the error bands shown in Fig. 2. The coefficients for these parametrizations are given in Table III. As a consistency check we use the parametrizations of the four independent cumulants and compare with the diagonal susceptibilities χ S 2 and χ B 2 , respectively. Note that for fixed N τ these data fulfill exactly the constraints given in Eqs. 11 and 12. The data for χ BQ 11 , χ BS 11 and χ QS 11 , however, are correlated, which may lead to slight differences in the continuum extrapolated results. We show the continuum extrapolations for χ B 2 and χ S 2 in Fig. 3. Here the error band has been obtained in the same way as for the off-diagonal cumulants. The continuum extrapolated results for these diagonal second order cumulants are given in Table IV for the same temperature values as those given in Table II for the set of four independent second order cumulants. Although the continuum extrapolated results for χ B 2 and χ S 2 have been obtained without imposing explicitly the relations given in Eqs. 11 and 12 we find that the continuum extrapolated cumulants are consistent with these constraints within errors and also the parametrization of second order cumulants given in Table III do so to better than 1%. [32] this work [32] this work [32] this work 135 0.0114 (5)(2)

V. SECOND ORDER CUMULANTS IN QCD AND THE HADRON RESONANCE GAS MODEL
In this section we compare QCD results for second order cumulants with HRG model calculations. We will discuss the importance of additional resonances contributing to correlations of net strangeness fluctuations with  net baryon-number and net electric-charge fluctuations, respectively. We also discuss in detail the quite different behavior of correlations between net baryon-number and net electric-charge fluctuations on the one hand and net baryon-number and net strangeness fluctuation on the other hand, which is apparent when comparing QCD and HRG model calculations. Furthermore, we comment on the significance of finite volume effects visible in the second order cumulant of net electric-charge fluctuations.
A. Sensitivity of second order cumulants to details of the hadron resonance gas spectrum Obviously, conclusions drawn from a comparison of HRG model calculations with QCD results crucially depend on the hadron spectrum which is input to the HRG model calculations. Although such models can rely on a lot of information from experimentally determined resonances [33], it has been noted that this information is not sufficient to constrain interactions in a hadron resonance gas to such an extent that these models do provide satisfactory comparisons with QCD results. As noted earlier [43], in particular in the baryon sector of the spectrum additional strange hadron resonances seem to be needed to obtain reasonable agreement between HRG model calculations and QCD results for strangeness fluctuations and their correlations with electric charge and baryon number fluctuations, respectively.
In particular, when using in HRG model calculations with point-like, non-interacting resonances only well established mesons and 3-and 4-star baryon resonances listed in the summary tables of the particle data group (PDG) [33], the comparison with QCD results yields only poor agreement in the strangeness sector (see Fig. 4 (top)). Including additional resonances that have been predicted in QCD based quark model (QM) calculations, e.g. in [34][35][36], as well as lattice QCD calculations [37,38], significantly improves the comparison between HRG model and QCD calculations. However, such an approach is not unique. It depends on details of the relativistic quark model calculations (for a recent compilation of results see [44]) as well as the treatment of the not well established resonances listed by the PDG. Moreover, it is well understood since the early work of Hagedorn [5] and Dashen, Ma and Bernstein [8] that a modeling of strong interactions in terms of point-like, interacting resonances is not sufficient to account for all interactions in strongly interacting matter. The need for a proper treatment of additional repulsive interactions, for instance by assigning an intrinsic volume to each hadron, has been pointed out early on [13,14]. However, it also has been noted that at the same time the interplay between repulsive and attractive interaction needs to be taken into account.
A more rigorous treatment of interactions among hadrons in medium can be achieved in the S-matrix formalism [8] which is the starting point for a systematic virial expansion of relativistic quantum gases [9,45]. The subtle interplay of attractive and repulsive interactions is, at least in principle, taken care of in the virial expansion. In the case of the strange meson K * 0 (700), for instance, an analysis of the effect of attractive and repulsive con-tributions in the S-matrix, partial wave analysis led to the conclusion that the contribution of this resonance to the thermodynamics of strong interaction matter is strongly suppressed. Effectively it does not contribute at all and the resonance thus should not be included in point-like, non-interacting HRG models [46,47], despite the fact that it is listed as a well established resonance in the PDG tables.
Setting up a HRG model for the description of the thermodynamics of strong interaction matter in the low temperature, hadronic regime thus is subject to a certain amount of ambiguity. Nonetheless, such models are a good starting point for the comparison to QCD calculations and can serve as a mediator between QCD calculations and experimental observations.

QMHRG2020
Previously we used a list of hadron resonances (QMHRG) that included only established mesons and 3and 4-star baryon resonances listed in the PDG summary tables and had been augmented with a list of QM states in the strange and non-strange baryon sectors [43]. We now updated this list of hadron resonances by taking into account the 1-and 2-star baryon resonances as well as mesons not listed as being well established in the PDG 2020 summary tables [33]. We, however, left out the K * 0 (700) (see below). This defines the list of hadron resonances, QMHRG2020 2 .
To a large extent the resonances, listed in the PDG tables, have counterparts in the hadron spectrum calculated in relativistic quark models [34][35][36]. In order to avoid double counting we used from the QM calculations only states that have no identified counterparts in the PDG tables. QMHRG2020 differs from the QMHRG2016+ list [49] by only a few resonances. In the strange baryon sector this leads to a few percent differences in the HRG model calculation of χ BS 11 , that are of significance when comparing HRG model calculations with QCD results, while in all other cases differences in the two lists are negligible. For instance, for temperatures below the QCD pseudo-critical temperature, i.e. for T ∼ (140 − 155) MeV, HRG model results for χ BS 11 , obtained from both lists, differ by about 6%. Eliminating the QM counterparts [34,36] of two 4-star as well as two 3-star resonances from the QMHRG2016+ list reduces this difference to 3%. In Table VI of Appendix C we give a list of 8 strange baryon resonances that are listed in the PDG tables and also have been calculated in relativistic quark models. Eliminating the latter from QMHRG2016+ reproduces results obtained with QMHRG2020 within 1% accuracy 3 .
In the following we use the QMHRG2020 list of hadron resonances as baseline model for comparisons with QCD calculations. In the insets of Fig. 2 we have compared the continuum extrapolated results for four independent second order cumulants, obtained in (2 + 1)-flavor QCD, with HRG model calculations that make use of the QMHRG2020 list. The insets show the ratio of results obtained in HRG model calculations and QCD, respectively. As can be seen, at low temperatures the agreement is fairly good. In detail, however, the differences seen in QCD and HRG model calculations are quite different for the four second order cumulants and reflect different physics. We will discuss this in more detail in the following subsections. We note already here that the difference between QCD results and HRG calculations is particularly striking for correlations between net baryonnumber and electric charge, χ BQ 11 . . Shown is a comparison to HRG model calculations based on different lists for hadron resonances as discussed in the text. Also shown are results obtained with excluded volume HRG models, using an excluded volume parameter b = 0.4 fm 3 , and virial expansions [11,50], respectively. 3 Both lists still differ to the extent that masses of strange baryons, obtained in quark model calculations, are taken from [36] in QMHRG2020 whereas QMHRG2016+ uses results from [34].

Net baryon-number fluctuations and correlations
In Fig. 4 we compare HRG model calculations using different hadron resonance spectra with continuum extrapolated lattice QCD results for (2+1)-flavor QCD. As can be seen correlations between net baryon-number and strangeness fluctuations (χ BS 11 ) are particularly sensitive to contributions from the strange baryon resonances, while correlations between net baryon-number and electric charge fluctuations (χ BQ 11 ) show only little sensitivity to contributions from additional baryon resonances; this cumulant only depends mildly on the strange baryon sector as only |S| = 2, 3 baryons contribute to χ BQ 11 . The small differences seen in the magnitude of χ BQ 11 when calculated with PDGHRG and QMHRG spectra, respectively, (Fig. 4 (bottom)) mainly arise from additional non-strange baryons obtained in QM calculations.
We note that |χ BS 11 |, calculated in a HRG model using QMHRG2020, is larger by about 30% relative to calculations based only on 3-and 4-star resonances listed by the PDG. HRG model calculations of χ BS 11 , using QMHRG2020, are consistent with QCD results at and below the pseudo-critical temperature. HRG model results for χ BQ 11 , on the other hand, are clearly larger than the corresponding QCD results for temperatures T > ∼ 145 MeV. They are about 20% larger at T pc,0 . This is a quite robust result, as HRG model calculations for χ BQ 11 are not very sensitive to the details of the HRG resonance spectrum used and even the inclusion of 1and 2-star resonances has little effect, as can be seen in Fig. 4 (bottom).
In the case of net baryon-number fluctuations, χ B 2 , which are related to χ BQ 11 and χ BS 11 through Eq. 12, this still leads to 10% larger results in HRG model calculations than QCD results obtained at T pc,0 .
In Fig. 5 we show for temperatures close to and below T pc,0 a comparison between HRG model calculations, performed with the QMHRG2020 and QMHRG2016+ resonance lists, and QCD. In this figure we also compare our results with results obtained by Bellwied et al. [32]. Within the current statistical and systematic uncertainties the results of both works are consistent with each other.
As outlined above, the larger magnitude obtained for |χ BS 11 | when using QMHRG2016+ arises from the fact that this list uses some quark model states that are also listed in the PDG tables. As a consequence the HRG model contains too many strange baryons; QMHRG2020 gives better agreement with QCD results below T pc,0 .
The differences between lattice QCD calculations and HRG model calculations using point-like, non-interacting hadronic resonances, indicate that these models do not correctly reflect interactions in a strongly interacting medium. It has been attempted to improve the model calculations by taking into account additional repulsive interactions in so-called excluded volume HRG models (EVHRG) [15,16,18]. Generically the magnitude of sec- ond order cumulants involving net baryon-number fluctuations is decreased in such excluded volume HRG models. Indeed, in the case of χ BQ 11 this may improve the comparison to QCD. However, at the same time it is obvious that EVHRG model calculations for χ BS 11 will spoil the good agreement observed between QCD results and HRG model calculations with point-like, non-interacting resonances observed below T pc,0 . While χ BS 11 favors only small excluded volumes, χ BQ 11 would require large volumes to achieve agreement between HRG model calculations and QCD results.
In fact, the relative change in second order cumulants involving net baryon-number fluctuations, calculated in HRG and EVHRG models, respectively, is identical [18] Here b parametrizes the size of the excluded volume for all baryons and P HRG B (T ) denotes the contribution of baryons and anti-baryons to the pressure in a resonance gas model for non-interacting point-like resonances. It contains all information on details of the hadron spectrum. The excluded volume parameter b is related to the hard-sphere radius (r) of a hadron through b = 16πr 3 /3.
The comparison of EVHRG model calculations and QCD results, shown in Figs. 4 and 5 for χ BS 11 and χ BQ 11 , makes it clear that no unique choice for b exists that would give better agreement between QCD and HRG models for both cumulants simultaneously. In fact, using Eq. 16 we may write Within the errors put on QCD results, this ratio should be consistent with unity for EVHRG models in order to be consistent with QCD results. This puts bounds on the magnitude of the excluded volume parameter b. Using the QCD result for χ BX 11 , X = S, Q, and including the combined statistical and systematic error ∆ X , we find the largest (b + ) and smallest (b − ) excluded volume parameters that would yield EVHRG model results consistent with QCD results, Here we used the relation between the baryonic part of the pressure and the second order cumulant of net baryon-number fluctuations in HRG models with pointlike hadrons, P B /T 4 = χ B 2 . At low temperatures these bounds are not very stringent as the hadronic medium is dilute. However, close to T pc,0 , i.e. for T = (150 − 155) MeV, we find that b ≤ 0.4 fm 3 is needed in order for EVHRG calculations to be consistent with QCD results for χ BS 11 , whereas from χ BQ 11 we find that b should be significantly larger, i.e. 1 fm 3 ≤ b ≤ 2 fm 3 .
In fact, the temperature dependence of χ BQ 11 favors b 2 fm 3 . As can be seen in Fig. 2 (top, left), in QCD calculations it is found that χ BQ 11 rises much more slowly with temperature than in HRG model calculations using QMHRG2020. A value of b, significantly larger than 1 fm 3 would be needed to account for this difference. In Fig. 6 we show results for the temperature derivative of χ BQ 11 shown in Fig. 2 (top, left). In the case of QCD these derivatives are obtained from the parametrization given in Eq. 15, using a bootstrap fit to the continuum extrapolated results shown in Fig. 2. As can be seen b > ∼ 2 fm 3 would be needed in an EVHRG calculation to reproduce the QCD results for the temperature derivative of χ BQ 11 . As excluded volume corrections are identical for all baryon correlations and fluctuations, it is instructive to look at ratios of second order cumulants involving net baryon-number fluctuations. In this case excluded volume corrections drop out. Any difference between HRG and QCD cumulants thus is of different origin. We also note that due to the relation given in Eq. 12, it suffices to analyze one ratio, e.g. χ BQ 11 /χ BS 11 . Other ratios, are then simply related to this ratio, e.g., Differences found between QCD and HRG model calculations of χ BQ 11 /χ BS 11 thus translate into corresponding differences for the other two ratios. Results for χ BQ 11 /χ BS 11 are presented in Fig. 7. This shows that deviations from HRG model calculations, which cannot be accounted for in excluded volume models using a single parameter b, become significant already at temperatures T 145 MeV. This, of course, is a direct consequence of the deviation of HRG model calculations from QCD results setting in for χ BQ 11 already at T 140 MeV, while at the same time results for χ BS 11 obtained in HRG model calculations are in good agreement with QCD up to T pc,0 .
A similar conclusion has been drawn from calculations of the second virial coefficients using the S-matrix approach [51], where it has been pointed out that the influence of repulsive interactions, which motivated an excluded volume ansatz, is subtle and quite different in various quantum number channels contributing in the partial wave analysis of the second virial expansion coefficient.
Although the S-matrix approach, which is based on experimental information on phase shifts contributing to the S-matrix, provides a rigorous formulation of the thermodynamics of an interacting hadron gas in the grand canonical ensemble and, at least in principle, does not require any a-priori information on a particular spectrum of resonances, it generally is difficult to arrive at first principle, quantitative results on properties of second order cumulants. In a systematic virial expansion of the partition function, expressed in terms of the S-matrix, even the calculation of the second virial coefficient often suffers from insufficient experimental information even for (bottom) with HRG model results based on QMHRG2020 (red band) and second order virial expansion results (green band) [11,50].
two-particle interactions. Moreover, at higher densities multi-particle interactions become of importance [8,9], which are poorly known and thus require approximations when comparing a S-matrix based analysis to e.g. first principle lattice QCD calculations [10][11][12].
In Fig. 8 we compare QCD results for second order cumulants involving net baryon-number fluctuations, with corresponding results from calculations of the second virial coefficient. In Ref. [11] the second virial co-efficient for the description of correlations between net baryon-number and net strangeness fluctuation has been obtained in a unitary, multi-channel analysis [52,53]. We show results from this analysis in Fig. 8 (top). As can be seen, agreement between S-matrix calculations and lattice QCD results is within (10-15)% below the pseudo-critical temperature. This result is similar to, but does at present not improve over results that can already be achieved within a HRG model calculation based on QMHRG2020.
A calculation of the second virial coefficient for the description of correlations between net baryon-number and net electric-charge fluctuations is more difficult as less information on the relevant interaction channels is known. In [50] χ BQ 11 has been analyzed in the S-matrix approach taking into account two body interactions arising from elastic π-N scattering, and a small contribution from the inelastic π-N → η-N channel. This partial calculation of the second virial expansion coefficient turns out not to be sufficient to achieve good agreement with QCD results [12,50]. Although it improves the comparison with QCD, deviations at low temperatures are about a factor two larger than in the S-matrix analysis of χ BS 11 . The inclusion of further interaction channels and contributions from higher order corrections clearly is needed.
At present the insufficient knowledge on scattering process contributing to χ BQ 11 can only be overcome by modeling contributions from three body and higher order interactions. This has been attempted in [12] by approximating ππN interactions by a simple, structure-less vertex. The strength of the interaction vertex in this model has been tuned to achieve agreement with lattice QCD results at T 150 MeV. However, comparing the results obtained in [12] with the precise QCD results presented here also at lower temperatures shows that agreement in a wider temperature range cannot be achieved with a temperature independent coupling strength for this interaction vertex. At lower temperatures, T ∼ 140 MeV, deviations from QCD are still about 10%.

11
In Fig. 2 (bottom, right) we have shown results for χ QS 11 which also are found to agree well with HRG model calculations using QMHRG2020. As mentioned in section V A 1 we did not include the scalar kaon resonance K * 0 (700) (previously κ in PDG) [54,55] in the QMHRG2020 list, although it is listed as an established resonance in the PDG tables. How to best accommodate for this resonance in thermodynamic calculations is much discussed [46,47,56].
Kaons give the most significant contribution to the correlations between net electric-charge and strangeness fluctuations. At low temperatures, T ∼ 130 MeV, already the ground state kaon and its P-wave excitation K * (892) contribute more than 80% to the 2 nd order cumulant χ QS 11 . All heavier strange mesons and baryons account for the remaining contribution to χ QS 11 . Adding the K * 0 (700) to the list of strange meson resonances would change the HRG model result for χ QS 11 by almost 10%. However, as has been shown in an analysis of the strangeness fluctuation cumulant χ S 2 [47], the contribution of K * 0 (700) is largely reduced when treating its contribution in a virial expansion that makes use of information on scattering amplitudes in the S-wave K-π channel.
The QCD results for χ QS 11 are shown in Fig. 2 and compared to our list of hadron resonances QMHRG2020, in which we do not include the K * 0 (700) resonance. In Fig. 9 we compare in more detail the HRG model calculations, with and without the K * 0 (700) resonance included, with QCD results. Here we also show the S-matrix calculation taken from [47], which includes interactions in the I = 1/2 and I = 3/2 S-wave channels. We combined the result obtained from the virial expansion with those strange hadron resonances from the QMHRG2020 list that are not taken care of in the S-matrix analysis. The analysis of interactions in the S-matrix formulation of an interacting hadron gas thus motivates that the K * 0 (700) resonance does not contribute to the thermodynamics of such a system and thus should be left out when using a gas of point-like, non-interacting only.

Strangeness neutrality and the strangeness chemical potential
The ratio χ BS 11 /χ QS 11 is directly related to the ratio of baryon and strangeness chemical potentials in a strangeness and electric charge neutral medium. As can be seen in Eq. 9 it controls the value of µ S /µ B in the isospin symmetric case (n Q /n B = 1/2), which also gives the dominant contribution in the case most relevant for comparison with heavy ion experiments, n Q /n B = 0.4. The good agreement found for χ BS 11 /χ QS 11 when calculated in lattice QCD and HRG models using QMHRG2020 (see Fig. 9 (bottom)), suggests that a determination of the ratio of chemical potentials from experimental data, using as an intermediate step HRG model based relations, is appropriate at least for small values of the baryon chemical potentials, for which the leading order Taylor expansion expressions (Eq. 8) are valid. The strong sensitivity of χ BS 11 on the strange baryon sector, on the one hand, and the small sensitivity of χ QS 11 on details of the spectrum, on the other hand, also suggests that µ S /µ B provides information on additional strange baryon resonances contributing to the thermodynamics of strongly interacting matter at the freeze-out temperature.
In Fig. 10 [11]. The upper figure shows the result of a virial expansion, based on the analysis of S-wave scattering contributions in the K-π channel to strangeness fluctuations [47], as discussed in the text.
onances with (QMHRG) and without (PDGHRG) additional strange baryon resonances included. As can be seen the ratio µ S /µ B obtained in QCD calculations when imposing the strangeness neutrality condition n S = 0 and n Q /n B = 0.4 differs by about 15% from HRG model calculations, that only use the PDGHRG states, but is in good agreement with QMHRG model calculations. This shows that the ratio of strangeness and baryon chemical potentials indeed is sensitive to the spectrum of strange hadron resonances in a strongly interacting medium. We will discuss this further in a forthcoming publication, where we show that higher order contributions to µ S /µ B are indeed negligible for µ B /T < ∼ 1, which makes this ratio a good observable for probing thermal conditions in a strongly-interacting medium.
B. Volume dependence of the cumulant of net electric-charge fluctuations, χ Q

2
As can be seen in Fig. 2  for the second order cumulant of net electric-charge fluctuations still differ from HRG model calculations. At such low temperatures electric charge fluctuations are dominated by the contribution from pions. In this case it is well known that finite volume effects in a pion gas with mass m π < ∼ T can lead to significant deviations from results obtained in the thermodynamic limit [57][58][59].
The lattice QCD calculations presented in Fig. 2 have been performed on lattices with a fixed ratio of spatial versus temporal extent, N σ /N τ = 4. I.e. the spatial extent of the physical volume, L = N σ a, changes with temperature, T = 1/N τ a, such that LT = 4 stays constant.
In the temperature range shown in Fig. 2, we have calculated χ Q 2 for a pion gas in a finite volume with LT = 4. This can be done directly using the partition function of a Bose gas in a finite volume [58,59]. In order to mimic periodic boundary conditions and cubic box sizes as they are used in lattice QCD calculations, we instead used a scalar field theory discretized on a lattice as described in [57]. For a non-interacting pion and kaon gas we find that deviations from infinite volume results are well parametrized by a straight line ansatz in the temperature interval of interest, T ∈ [130 MeV : 170 MeV], 126 T /T pc,0 , pion gas 1.002 − 0.032 T /T pc,0 , kaon gas (21) At the pseudo-critical temperature, T pc,0 , the net electric-charge fluctuations in a pion gas in a volume LT = 4 thus is about 12% smaller than in the infinite volume limit. In the HRG model, calculated with the resonance spectrum QMHRG2020, this distortion effect gets reduced by almost a factor two, reflecting the relative contribution of pions to the entire net electriccharge fluctuations. In the temperature interval T ∈ [130 MeV : 180 MeV] we find for a HRG model using First error in QCD is obtained by linearly combining statistical and systematical error of our continuum extrapolations while the second one error reflects the error on the determination of Tpc,0. Similarly the error in the HRG is due to the error of Tpc,0. These results agree well with the continuum extrapolated data given in [32]. For χ Q 2 we give results calculated in a finite volume LT ≡ Nσ/Nτ = 4. Numbers in brackets give the corresponding infinite volume result.
the QMHRG2020 list of hadrons, This amounts to a 6% finite volume correction for χ Q 2 at the pseudo-critical temperature, T pc,0 , which only changes slowly as function of temperature. We show this finite volume correction to the QMHRG2020 result for χ Q 2 in Fig. 2 (top, right). Using the finite volume corrected QMHRG2020 results for a comparison with QCD results in a finite volume, we find that at T pc,0 the latter are still smaller by about 5%. this is consistent with large deviations observed in the charged baryon sector (χ BQ 11 ) which contributes only about 15% to the total electric charge fluctuations. Continuum extrapolated lattice QCD results for all six second order cumulants at T pc,0 and corresponding results from HRG model calculations, using different hadronic resonance spectra, as well as results from Smatrix calculations are summarized in Table V. Here we also give results for two of the three independent ratios of second order cumulants. The third ratio would involve χ Q 2 , for which at present no infinite volume extrapolated result exists.

VI. CONCLUSIONS
We have presented an update of continuum extrapolated results for second order cumulants of conserved charge fluctuations and their cross-correlations. These results are based on high-statistics data obtained in simulations within the HISQ discretization scheme for staggered fermions. Our results are found to be consistent with previous results obtained by using the stout discretization scheme for staggered fermions [32].
We compiled a new list of hadron resonances, QMHRG2020, which differs from QMHRG2016+ in particular in the strange baryon sector. HRG model calculations with QHMHRG2020 provide a good description of strangeness fluctuations and correlations with baryon-number and electric charge fluctuations, respectively. Deviations are found to be less than 10% in the temperature range 135 MeV < T < T pc,0 . This puts stringent bounds on the magnitude of excluded volume parameters for strange baryon interactions in EVHRG models.
The largest differences between QCD results and HRG model calculations have been found for correlations between net baryon-number and electric charge fluctuations. They amount to more than 20% at T pc,0 when comparing QCD with HRG models based on point-like, non-interacting resonances. Modeling these deviations, and in particular the quite different temperature dependence of χ BQ 11 found in the vicinity of T pc,0 , in excluded volume HRG models would require a large excluded volume parameter b > 1 fm 3 .
Calculations based on a virial expansion capture basic features of the interplay between repulsive and attractive interactions in the strongly interacting hadronic medium. They successfully explain why in particular the strange meson resonance K * 0 (700) does not contribute in that medium. They also qualitatively describe the smaller magnitude of χ BQ 11 found in QCD compared to HRG model calculations with point-like, non-interacting resonances. However, in order to achieve quantitative agreement with QCD modeling of interactions that are not captured in calculations of the second virial coefficient is needed.
The results presented here form the basis for a systematic analysis of second order cumulants at non-vanishing values of the chemical potentials. This will be discussed in a forthcoming publication. All data from our calculations, presented in the figures of this paper, can be found in [48].  Continuum extrapolations for the cumulants χ BS 11 , χ BQ 11 , χ QS 11 and χ Q 2 (top to bottom) performed with the temperature scales obtained from fK (left) and r1 (right) using linear and quadratic extrapolations in 1/N 2 τ . Shown are results for four values of the temperature below and close to the pseudo-critical temperature of (2 + 1)-flavor QCD. Also shown are results obtained in HRG model calculations using the QMHRG2020 list of hadrons (crosses). Note that for χ BQ 11 and χ Q 2 these HRG results are not shown for all temperature values as the deviations from the corresponding QCD results are too large. For χ Q