Off-diagonal correlators of conserved charges from lattice QCD and how to relate them to experiment

Like fluctuations, nondiagonal correlators of conserved charges provide a tool for the study of chemical freeze-out in heavy ion collisions. They can be calculated in thermal equilibrium using lattice simulations, and be connected to moments of event-by-event net-particle multiplicity distributions. We calculate them from continuum-extrapolated lattice simulations at μB 1⁄4 0, and present a finite-μB extrapolation, comparing two different methods. In order to relate the grand canonical observables to the experimentally available netparticle fluctuations and correlations, we perform a hadron resonance gas model analysis, which allows us to completely break down the contributions from different hadrons. We then construct suitable hadronic proxies for fluctuation ratios, and study their behavior at finite chemical potentials. We also study the effect of introducing acceptance cuts, and argue that the small dependence of certain ratios on the latter allows for a direct comparison with lattice QCD results, provided that the same cuts are applied to all hadronic species. Finally, we perform a comparison for the constructed quantities for experimentally available measurements from the STAR Collaboration. Thus, we estimate the chemical freeze-out temperature to 165 MeV using a strangeness-related proxy. This is a rather high temperature for the use of the hadron resonance gas; thus, further lattice studies are necessary to provide first principle results at intermediate μB.


I. INTRODUCTION
The study of the phase diagram of quantum chromodynamics (QCD) has been the object of intense effort from both theory and experiment in the last decades. Relativistic heavy ion collision experiments both at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) have been able to create the quark gluon plasma (QGP) in the laboratory, and explore the low-to-moderate baryon density region of the QCD phase diagram.
At low baryon density, the transition from a hadron gas to a deconfined QGP was shown by lattice QCD calculations to be a broad crossover [1] at T ≃ 155 MeV [1][2][3][4]. At large baryon densities, the nature of the phase transition is expected to change into first order, thus implying the presence of a critical end point. A strong experimental effort is currently in place through the second Beam Energy Scan (BES-II) program at the RHIC in 2019-2021, with the goal of discovering such a critical point.
The structure of the QCD phase diagram cannot currently be theoretically calculated from first principles, as lattice calculations are hindered by the sign problem at finite density. Several methods have been utilized in order to expand the reach of lattice QCD to finite density, like full reweighting [5], Taylor expansion of the observables around μ B ¼ 0 [6][7][8][9][10], or their analytical continuation from imaginary chemical potential [5,[11][12][13][14].
We remark here that there are alternative approaches to lattice QCD for the thermodynamical description. Specific truncations of the Dyson-Schwinger equations allow the calculation of the crossover line and also to extract baryonic fluctuations [15,16]. Another theoretical result on the baryon-strangeness correlator has been calculated using functional methods from the Polyakov-loop-extended quark meson model in [17,18].
The confined, low-temperature regime of the theory is well described by the hadron resonance gas (HRG) model, which is able to reproduce the vast majority of lattice QCD results in this regime [19][20][21][22][23][24]. Moreover, the HRG model has been extremely successful in reproducing experimental results for particle yields over several orders of magnitude [25][26][27][28][29]. These are usually referred to as thermal fits, since the goal of the procedure is the determination of the temperature and chemical potential at which the particle yields are frozen. This moment in the evolution of a heavy ion collision is called chemical freeze-out, and takes place when inelastic collisions within the hot hadronic medium cease. The underlying assumption is that the system produced in heavy ion collisions eventually reaches thermal equilibrium [30][31][32], and therefore a comparison between thermal models and experiment is possible [33].
Although the net number of individual particles may change after the chemical freeze-out through resonance decays, the net baryon number, strangeness and electric charge are conserved. Their event-by-event fluctuations are expected to correspond to a grand canonical ensemble. In general, when dealing with fluctuations in QCD, and in particular in relation to heavy ion collisions, it is important to relate fluctuations of such conserved charges and the event-by-event fluctuations of observed (hadronic) species. The former have been extensively studied with lattice simulations [34][35][36][37][38][39][40][41][42][43], and are essential to the study of the QCD phase diagram for multiple reasons. First, they are directly related to the Taylor coefficients in the expansion of the pressure to finite chemical potential and have been utilized to reconstruct the equation of state of QCD at finite density, both in the case of sole baryon number conservation [41,44,45], and with the inclusion of all conserved charges [46]. Second, higher order fluctuations are expected to diverge as powers of the correlation length in the vicinity of the critical point, and have thus been proposed as natural signatures for its experimental search [47,48]. On the other hand, fluctuations of observable particles can be measured in experiments, and are very closely related to conserved charge fluctuations. With some caveats [49][50][51][52], comparisons between the two can be made, provided that certain effects are taken into account.
Previous studies found that, for certain particle species, fluctuations are more sensitive to the freeze-out parameters than yields [53]. In recent years, the STAR Collaboration has published results for the fluctuations of net proton [54], net charge [55], net kaon [56], and more recently net Λ [57] and for correlators between different hadronic species [57,58]. From the analysis of net-proton and net-charge fluctuations in the HRG model, it was found that the obtained freeze-out temperatures are lower than the corresponding ones from fits of the yields [52]. More recent analyses [59,60] of the moments of net-kaon distributions showed that it is not possible to reproduce the experimental results for net-kaon fluctuations with the same freeze-out parameters obtained from the analysis of net proton and net charge. In particular, the obtained freeze-out temperature is consistently higher, with a separation that increases with the collision energy. In [59], predictions for the moments of net-Λ distributions were provided, calculated at the freezeout of net kaons and net proton/net charge.
Correlations between different conserved charges in QCD provide yet another possibility for the comparison of theory and experiment. They will likely receive further contribution from measurements in the future, with new species being analyzed and increased statistics allowing for better determination of moments of event-by-event distributions [58].
In this manuscript, we present continuum-extrapolated lattice QCD results for all second-order nondiagonal correlators of conserved charges. We then identify the contribution of the single particle species to these correlators, distinguishing between measured and nonmeasured species. Finally, we identify a set of observables, which can serve as proxies to measure the conserved charge correlators. The manuscript is organized as follows. In Sec. II we present the continuum extrapolated lattice results for second-order nondiagonal correlators of conserved charges and discuss the extrapolation to finite μ B in Sec. III. In Sec. IV we show the comparison with HRG model calculations, and describe the breakdown of the different contributions to the observables shown in the previous section. In Sec. V we propose new observables which can serve as proxies to directly study the correlation of conserved charges. In Sec. VI we analyze the behavior of the constructed proxies at finite chemical potential, and study the effect of acceptance cuts in the HRG calculations. We argue that the small dependence on experimental effects allows for a direct comparison with lattice QCD results. We also perform a comparison to experimental results for selected observables in Sec. VII. Finally, in Sec. VIII we present our conclusions.

II. LATTICE QCD AND THE GRAND CANONICAL ENSEMBLE
The lattice formulation of quantum chromodynamics opens a nonperturbative approach to the underlying quantum field theory in equilibrium. Its partition function belongs to a grand canonical ensemble, parametrized by the baryochemical potential μ B , the strangeness chemical potential μ S and the temperature T. Additional parameters include the volume L 3 , which is assumed to be large enough to have negligible volume effects, and the quark masses. The latter control the pion and kaon masses, and are set to reproduce their physical values. At the level of accuracy of this study we can assume the degeneracy of the light quarks m u ¼ m d and neglect the effects coming from quantum electrodynamics.
There is a conserved charge corresponding to each flavor of QCD. The grand canonical partition function can be then written in terms of quark number chemical potentials (μ u , μ d , μ s ). The derivatives of the grand potential with respect to these chemical potentials are the susceptibilities of quark flavors, defined as withμ q ¼ μ q =T. These derivatives are normalized to be dimensionless and finite in the complete temperature range. For the purpose of phenomenology we introduce for the B (baryon number), Q (electric charge) and S (strangeness) a chemical potential μ B , μ Q and μ S , respectively. The basis of μ u , μ d , μ s can be transformed into a basis of μ B , μ Q , μ S using the B, Q and S charges of the individual quarks, Susceptibilities are then defined as It is straightforward to express the derivatives of p=T 4 with respect to μ B , μ Q and μ S in terms of the coefficients in Eq. (1) [36,39,61]. For the cross-correlators we have Such derivatives play an important role in experiment. In an ideal setup the mean of a conserved charge i can be expressed as the first derivative with respect to the chemical potential, while fluctuations and cross-correlators (say between charges i and j) are second derivatives, In these formulas N i indicates the net number of charge carriers. For example, for the baryon number B, it corresponds to the number of baryons minus the number of antibaryons. The procedure to define the chemical potential on the lattice [62] and to extract the derivatives in Eq. (1) from simulations that run at μ u ¼ μ d ¼ μ s ¼ 0 has been worked out long ago [6] and has been the basis of many studies ever since [7,61,63,64]. Since the derivatives with respect to the chemical potential require no renormalization, a continuum limit could be computed as soon as results on sufficiently fine lattices emerged [3,36]. Later the temperature range and the accuracy of these extrapolations were extended in [39,40].
In this work, we extend our previous results [39] to nondiagonal correlators and calculate a specific ratio that will be later compared to experiment.
These expectation values are naturally volume dependent. Their leading volume dependence can, however, be canceled by forming ratios. In [37,38,65] such ratios were formed between various moments of electric charge fluctuations, and also for baryon fluctuations. For the same ratios the STAR experiment has provided proxies as part of the Beam Energy Scan I program [54,55].
The gauge action is defined by the tree-level Symanzik improvement, and the fermion action is a one link staggered with four levels of stout smearing. The parameters of the discretization as well as the bare couplings and quark masses are given in [39].
The charm quark is also included in our simulations, in order to account for its partial pressure at temperatures above 200 MeV, where it is no longer negligible [66]. In the range of the expected chemical freeze-out temperature between 135 and 165 MeV the effect of the charm quark is not noticeable on the lighter flavors [39].
In this work we use the lattice sizes of 32 3 × 8, 40 3 × 10, 48 3 × 12, 64 3 × 16 as well as 80 3 × 20. Thus, the physical volume L 3 is given in terms of the temperature as LT ¼ 4 throughout this paper. The coarsest lattice was never in the scaling region. The finest lattice lacks the precision of the others, and we only use it when the coarser lattices, e.g., 40 3 × 10, are not well in the scaling region and if, for a particular observable, the 80 3 × 20 has small enough error bars. If this was used, the data set is also shown in the plots.
We show here the continuum extrapolated cross-correlators at zero chemical potential. In Fig. 1 we show an example of continuum extrapolation for the three crosscorrelators, with T ¼ 150 MeV and the w 0 -based scale setting [67]. Figure 2 shows χ BQ 11 ðTÞ for the four different lattices, as well as the continuum extrapolation. Although our simulation contains a dynamical charm quark, we did not account for its baryon charge. Thus, the Stefan-Boltzmann limit of this quantity is 0. This limit is reached when the mass difference between the strange and light quarks becomes negligible in comparison to the temperature. The peak is seen at a higher temperature than T c ≈ 155 MeV, while in the transition region there is an inflection point. Below T c this correlator is dominated by protons and charged hyperons. In Sec. IV we account in detail for various hadronic contributions in the confined phase.  The χ BS 11 ðTÞ correlator is shown in Fig. 3. Unlike the BQ correlator, we have now a monotonic function with a high temperature limit of −1=3, which is the baryon number of the strange quark. The transition has a remarkably small effect on this quantity. At low temperatures this correlator is basically the hyperon free energy.
The χ QS 11 ðTÞ correlator in Fig. 4 is also monotonic, converging to 1=3 at high T, which is the electric charge of the strange quark. At low temperatures this quantity is dominated by the charged kaons, which were the focus of recent experimental investigations [56].
The continuum-extrapolated results at μ B ¼ 0 for χ BQ 11 , χ BS 11 , χ QS 11 and the ratio χ BS 11 =χ S 2 we discuss below are provided in Table I of the Appendix.

III. RESULTS AT FINITE DENSITY
Since lattice QCD can be defined at finite values of the B, Q and S chemical potentials, and is capable of calculating derivatives of the free energy as a function of these chemical potentials one could expect that the extension of the simulations to finite density is a mere technical detail. Unfortunately, at any finite real value of the quark chemical potential μ q the fermionic contribution to the action becomes complex and most simulation algorithms break down.
There are several options to extract physics at finite densities, nevertheless. It seems natural to use algorithms that were designed to work on complex actions-both the complex Langevin equation [68] and the Lefschetz thimble approach [69] have shown promising results recently-yet their direct application to phenomenology requires further research.
Instead, we use here the parameter domain that is available for mainstream lattice simulations. In fact, besides zero chemical potential, simulations at imaginary μ B are also possible, and have been exploited in the past to extrapolate the transition temperature [70][71][72][73], fluctuations of conserved charges [42,43] and the equation of state [41]. In all these works it was assumed that the thermodynamical observables are all analytical functions ofμ 2 B . A conceptually very similar method, the Taylor method, provides the extrapolation in terms of calculating higher derivatives with respect toμ B . The series is truncated at a certain order, which is typically limited by the statistics of the lattice simulation.
Since we relate the baryon-strangeness correlator to experimental observables later on we use χ BS 11 as an example for the Taylor expansion, where the terms proportional toμ B andμ S vanish since they contain odd derivatives, which are forbidden by the C-symmetry of QCD. (The charge chemical potential is omitted for simplicity).
In most phenomenological lattice studies the chemical potentials are selected such that the strangeness vanishes for each set of ðμ B ; μ Q ; μ S Þ. More precisely, for each T and μ B we selectμ Q andμ S values such that The factor 0.4 is the typical Z=A ratio for the projectiles in the heavy ion collision setup, and the value we use in the HRG model calculations below. In our lattice study, however, we use 0.5. This introduces a small effect compared to the statistical and systematic errors of the extrapolation, and results in substantial simplification of the formalism: μ Q can be chosen to be 0. The would-be μ Q value is about one tenth of μ S in the transition region [37,38].
We checked the impact of our simplification on the results we present here, with the Taylor expansion method. Utilizing our simulations data from ensembles at μ B ¼ 0, we calculate the correction to the ratio χ BS 11 =χ S 2 . By construction, this correction vanishes at μ B ¼ 0, and we find that it grows to at most 1% at μ B =T ¼ 1, and at most 1.5% at μ B =T ¼ 2. These systematic errors are considerably smaller than the uncertainties we have on our results, as can be seen in Fig. 5.

A. Taylor method
The Taylor coefficients for correlators can be easily obtained by considering the higher derivatives with respect to μ B . For later reference we select the quantity χ BS 11 =χ S 2 for closer inspection, up to Oðμ 4 B Þ corrections, with The derivatives on the right-hand side are all taken at Whether we extract the required derivatives from a single simulation (one per temperature) at μ B ¼ μ S ¼ 0 or we determine the fourth order derivatives numerically by the (imaginary) μ B -dependence of second-order derivatives is a result of a cost benefit analysis. The equivalence of these two choices has been shown on simulation data (of the chemical potential dependence of the transition temperature) in [73].
In [71] we calculated direct derivatives; however, we obtained smaller errors by using imaginary μ B simulations in [43]. Thus, we take the Taylor coefficients from the latter analysis, now extended to the new observable. The results for several fixed temperatures are shown in Fig. 5. As a first observable we show the μ S =μ B ratio that realizes strangeness neutrality. Then we show the χ BS 11 =χ S 2 ratio as a function of positiveμ 2 B . In the plots we show results from a specific lattice 48 3 × 12, which has the highest statistics, so that it provides the best ground to compare different extrapolation strategies.

B. Sector method
In Fig. 5 we compare two extrapolation strategies; here we describe the second approach, the sector method.
We are building on our earlier work in [74], where we have written the pressure of QCD as a sum of the sectors These sectors were also studied in [75,76]. Obviously, QCD receives contributions from sectors with higher quantum numbers as well. The sectors in Eq. (17) are the only ones receiving contributions from the ideal hadron resonance gas model in the Boltzmann approximation. (The dependence on the electric charge chemical potential is not considered now, since we selected The partitioning of the QCD pressure in sectors is very natural in the space of imaginary chemical It is expected that higher sectors are increasingly relevant as T c is approached from below. A study using Wuppertal-Budapest simulation data has shown that below T ≈ 165 MeV the sectors jBj ¼ 0, 1, 2 give a reasonable description, e.g., by calculating χ B 4 from the sector coefficients and comparing to direct results.
Thus, for this work we considered the next-to-leading order of the sector expansion, including the LO∶ P BS 01 ; P BS 11 ; P BS 12 ; P BS 13 ; NLO∶ P BS 02 ; P BS 1;−1 ; P BS 21 ; P BS 22 : It is somewhat ambiguous how the NLO is to be defined. One option would be to include the next-higher jBj quantum number, making our approach second order in this expansion. We do include P BS 21 and P BS 22 ; however adding further higher strangeness sectors, e.g., P BS 23 , did not improve the agreement with our data, so it was not included. Moreover, we included the multistrange sector P BS 02 and the exotic sector P BS 1;−1 , which is the coefficient of the term proportional to coshðμ B þμ S Þ. On the other hand, removing the sectors included in the NLO never resulted in higher temperatures in a smaller χ 2 for the fit; e.g., at T ¼ 160 MeV removing the terms P BS 0;2 , P BS 1;−1 P BS 2;1 or P BS 2;0 from the fit (removing only one term at a time) resulted in a χ 2 =N dof of 72.6=53, 181.1=53, 72.4=53, or 137.2=53, respectively, while including all gave 72.3=52.
In Fig. 6 the results for the sectors in the LO and NLO are shown at different temperatures. The results for the LO sectors shown in the top panel (figure from Ref. [74]) are continuum extrapolated-except for the jBj ¼ 0, jSj ¼ 1 sector. In the lower panel new results for the sectors in the NLO are shown, generated from a 48 3 × 12 lattice.
The alert reader may ask why we do not include the P BS 10 sector, accounting e.g., for protons. In fact, the sectors with jSj ¼ 0 do not contribute to the observables χ BS 11 , χ S 2 or χ S 1 either at 0, or at any real or imaginary chemical potential. In the analysis we included results for χ S 1 , χ S 2 , χ BS 11 from various data sets at various imaginary chemical potentials: at [43], the strangeness data set with μ I B > 0 of [41], and finally the set (only including χ S 1 and χ S 2 ) with μ B ¼ 0, μ I S > 0 from [74]. For the lower temperatures the model defined with the coefficients in Eq. (19) resulted in good fits (Q values ranging from 0 to 1)-the worst fit was at T ¼ 165 MeV with Q ≈ 0.05. This is the temperature where the model is expected to break down. Now we can compare the results to the Taylor expansion. In Fig. 5 we show the sector results with error bars, while the bands refer to the Taylor method. At low temperatures we see good agreement even for large values of the chemical potential; near the transition, however, the two approaches deviate already in the experimentally relevant region. It is obvious that the sector method breaks down above T c . Its systematic improvement to higher jBj quantum numbers requires much higher statistics (the same is true for the Taylor coefficients). Each further order enables the extrapolation to somewhat larger chemical potential, and in the case of the jBj sectors, to a somewhat larger temperature. Let us note that the Taylor method has limitations as well, slightly above T c , because the subsequent orders are not getting smaller [43]. The reason for this behavior is the fact that, between T ¼ 160 MeV and T ¼ 180 MeV, there is a crossover transition in the imaginary domain ofμ B ; then higher Taylor coefficients facilitate an extrapolation through that crossover.
In conclusion, we consider only the chemical potential range where our two methods agree in the extrapolation. At present, our lattice data allow a continuum extrapolation from the sector method only, which we do using 40 3 × 10, 48 3 × 12 and 64 3 × 16 lattices in the temperature range 135-165 MeV for a selection of fixed realμ B values. The method for the continuum extrapolation is the same we also used in Sec. II. We show the result in Fig. 7.
The large error bars in comparison to the μ B ¼ 0 results and the limited range inμ B indicate that the extraction of finite density physics from μ B ¼ 0 or imaginary μ B simulations is a highly nontrivial task. Still, both the Taylor and the sector methods can be systematically improved to cover more of the range of interest for the Beam Energy Scan II program. Given the high chemical freeze-out temperatures for jSj ¼ 1 particles-see Sec. VII-that emerge from STAR data ( [56] and preliminary [57]), the use of continuum extrapolated lattice simulations to  [74]. In the second panel we show the nonstandard sectors in a linear scale that we use for the μ B extrapolation in this work. Note that results from the same lattice for the jBj-only sectors are shown in [77], and can be taken as a comparison. calculate the grand canonical features of QCD is highly motivated.

IV. CORRELATORS IN THE HRG MODEL
The HRG model is based on the idea that a gas of interacting hadrons in their ground state can bewell described by a gas of noninteracting hadrons and resonances. The partition function of the model can thus be written as a sum of ideal gas contributions of all known hadronic resonances R, where every quantity with a subscript R depends on the specific particle in the sum. The relativistic energy is and the conserved charges B R , Q R and S R are the baryon number, electric charge and strangeness respectively. Moreover, d R is the spin degeneracy, m R is the mass, and the factor η R ¼ ð−1Þ 1þB R is 1 for (anti-)baryons and −1 for mesons.
The temperature and the three chemical potentials are not independent, as the conditions in Eqs. (12) are imposed on the baryon, electric charge and strangeness densities. We use these constraints to set both μ Q ðT; μ B Þ and μ S ðT; μ B Þ in our HRG model calculations.
In this work we utilize the hadron list PDG2016+ from [74], which was constructed with all the hadronic states (with the exclusion of charm and bottom quarks) listed by the Particle Data Group (PDG), including the less-established states labeled by *, ** [78]. The decay properties of the states in the list, when not available (or complete) from the PDG, were completed with a procedure explained in [79], and then utilized in [59,79].
In the HRG model the χ BQS ijk susceptibilities of Eq. (5) can be expressed as where B R , Q R , S R are the baryon number, electric charge and strangeness of the species R and the phase space integral at order i þ j þ k reads (note that it is completely symmetric in all indices, hence i þ j þ k ¼ l), The HRG model has the advantage, when comparing to experiment, of allowing for the inclusion of acceptance cuts and resonance decay feed-down, which cannot be taken into account in lattice QCD calculations.
The acceptance cuts on transverse momentum and rapidity (or pseudorapidity) can be easily taken into account in the phase space integrations via the change(s) of variables, in the case of rapidity and pseudorapidity respectively, where in all cases the trivial angular integrals were carried out [52].
A. Correlators of measured particle species The rich information contained in the system created in a heavy ion collision about the correlations between conserved charges is eventually carried over to the final stages through hadronic species correlations and self-correlations. It is convenient, in the framework of the HRG model, to consider the hadronic species which are stable under strong interactions, as these are the observable states accessible to experiment. However, due to experimental limitations, charged particles and lighter particles are easier to measure, and so we cannot access every relevant hadron related to conserved charges. Thus, historically protons have served as a proxy for baryon number, kaons have served as a proxy for strangeness, and net electric charge is measured through p, π, and K.
In our framework, we consider the following species, stable under strong interactions: π 0 , π AE , K AE , K 0 ,K 0 , p,p, n,n, Λ,Λ, Σ þ ,Σ − , Σ − ,Σ þ , Ξ 0 ,Ξ 0 , Ξ − ,Ξ þ , Ω − ,Ω þ . Of these, the commonly measured ones are the following: π AE ; K AE ; pðpÞ; ΛðΛÞ; Ξ − ðΞ þ Þ; Ω − ðΩ þ Þ: A few remarks are in order here. First of all, we refer to the listed species as commonly measured because, although some others are potentially measurable (especially the charged Σ baryons), results for their yields or fluctuations are not routinely performed both at the RHIC and the LHC. In the following, we keep our nomenclature of measured and "nonmeasured" in accordance to the separation we adopt here. Obviously, neutral pions can be measured with the process π 0 → γγ, but they are not included here as they do not carry any of the conserved charges of strong interactions. An additional note is necessary for K 0 S : although the measurement of K 0 S is extremely common in experiments, it is not of use for the treatment we carry on in this work. This is because, from K 0 S only, it is not possible to construct a net-particle quantity (it is its own antiparticle), and additionally part of the information on the mixing between K 0 andK 0 is lost because K 0 L cannot be measured. For this reason, in the following we consider K 0 andK 0 instead, and treat them as "not measured." Finally, we note that, since the decay Σ 0 → Λ þ γ has a branching ratio of ∼100%, effectively what we indicate with Λ contains the entire Σ 0 contribution as well; this well reproduces the experimental situation, where Λ and Σ 0 are treated as the same state.
It is straightforward to adapt the HRG model so that it is expressed in terms of stable hadronic states only. The sum over the whole hadronic spectrum is converted into a sum over both the whole hadronic spectrum and the list of states which are stable under strong interactions, with l þ m þ n ¼ p, and where the first sum only runs over the particles which are stable under strong interactions, and the sum P R→i ¼ P α N α R→i n R i;α gives the average number of particle i produced by each particle R after the whole decay chain. The sum runs over particle R decay modes, where N α R→i is the branching ratio of the mode α, and n R i;α is the number of particles i produced by a particle R in the channel α.
In light of the above considerations, it is useful to define the contribution to the conserved charges from final state stable hadrons. In the following, we adopt the convention where the net number of particles of species A (i.e., the number of particles A minus the number of antiparticlesĀ) With this definition, we can express conserved charges as net-B∶p þñ þΛ þΣ þ þΣ − þΞ 0 þΞ − þΩ − ; Using this decomposition, we can write as an example the BQ correlator where P R→net−B ¼P R→p þP R→ñ þP R→Λ þP R→Σ þ þP R→Σ − þ P R→Ξ 0 þP R→Ξ − þP R→Ω − , and e.g., P R→p ¼ P R→p − P R→p . Analogous expressions apply to net-Q and net-S.
The result of this decomposition is that each of the correlators one can build between conserved charges is formed from the sum of many different particle-particle correlations. In particular, the sum of those correlators which entirely consist of observable species yields the measured part of a certain correlator, while its nonmeasured part consists of all other terms, which include at least one nonobservable species. In Fig. 8 the nondiagonal correlators are shown as a function of the temperature at vanishing chemical potential. The measured and nonmeasured contributions are shown with blue, dashed-dotted and red, dashed lines, respectively, while the full contribution is shown with a solid, thicker black line. Alongside the HRG model results, continuum extrapolated lattice results are shown as magenta points as introduced in Sec. II.
We notice that both the BQ and QS correlators are largely reproduced by the measured contribution (for the BQ correlator, the measured portion even exceeds the full one, as the nonmeasured contribution is negative), while the BS correlator is roughly split in half between measured and nonmeasured terms. This is because the former are unsurprisingly dominated by the net-proton and net-kaon contributions, respectively, which in this temperature regime form the bulk of particle production, together with the pions. The BS correlator conversely receives its main contributions from strange baryons, which are almost equally split between measured and nonmeasured.

B. Breakdown of the measured and nonmeasured contributions
The decomposition in Eq. (26) allows one to break down the different contributions to any cross correlator, as well as the diagonal ones, entirely. In Figs. 9 and 10, we show the breakdown of the measured portion of the single final state hadronic (self-) correlations to the nondiagonal and diagonal correlators, respectively. Let us start from the nondiagonal case.
A few features can be readily noticed. First, in all cases only a handful of the most sizable contributions account for  the measured portion of the corresponding observable. As stated above, the BQ and QS correlators are expected to be dominated by the contribution from net-proton and netkaon self-correlations, respectively: indeed, in both cases the measured part almost entirely consists of these major contributions. Second, it is worth noticing how, with the only exception of the proton-pion correlator within χ BQ 11 , all correlators between different species yield a very modest contribution. This is the case for the proton-kaon, kaonpion, Lambda-pion and Lambda-kaon correlators in χ BQ 11 and χ QS 11 , as well as theproton-kaon, Lambda-kaon and Lambda-proton correlators in χ BS 11 . In our setup, correlations between different particle species can only arise from the decay of heavier resonances. Whenever a resonance R has a nonzero probability to decay, after the whole decay cascade, both into stable species A and B, then a correlation arises between A and B. It can be seen from Eq. (27) that only when both probabilities in parentheses are nonzero can a nonzero correlation arise. For the same reason, correlations between different baryons arise, although no single decay mode with more that one baryon (or antibaryon) is present in our decay list. In fact, if a state exists which has a finite probability to produce-after the whole decay cascade-both baryon A and baryon B, then a correlation between A and B is generated through Eq. (27). Finally, since both Ξ − and Ω − carry all three conserved charges, they contribute to all three correlators through their self-correlations, and their contribution is not negligible in all cases.
The case of χ BS 11 is slightly different, as the measured part is smaller than in the cases of χ QS 11 and χ BQ 11 . This is due to the fact that there is no significant separation in mass between the lightest observable particle carrying both baryon number and strangeness-the Λ baryon-and the lightest of the nonmeasured ones-the Σ AE baryons. In fact, the contribution from both charged Σ baryons is comparable to the one from the Λ, which thus cannot play as big of a role as the proton and kaon in the other two correlators, as well as because of the previously mentioned fact that correlators of different species do not contribute significantly.
In the diagonal case, a similar picture appears. The χ Q 2 correlator is almost identical to its measured portion, dominated by the self-correlations of pions, kaons, and protons. The other two correlators have a similar situation to that of χ BS 11 , with the measured part roughly amounting to half of the total. Again, the only non-negligible correlator between different species is the proton-pion correlator in χ Q 2 . We notice that in general, the leading single contribution is not as close to the whole measured portion, as it was in the case of the cross-correlators. This aspect is important in the following, where we move to the analysis of ratios of correlators, and look for suitable proxies.
In Figs. 11 and 12 we show the breakdown of the nonmeasured portion of the final state hadronic (self-) correlations, analogously to what we showed in Figs. 9 and 10 for the measured portion. The situation in this case is slightly different from the previous one: it is generally more difficult to identify a leading contribution, with multiple terms yielding comparable results. In the case of BQ and QS, leading terms come with opposite signs, which further complicates the picture. The reasons for these features come from the fact that (i) the number of single contributions that are not measured is much larger than that of the measured ones, hence it is less probable that few terms dramatically dominate; and (ii) in general, nonmeasured species are heavier than the measured ones, hence single contributions tend to be smaller. Obviously, exceptions to this are the neutron and K 0 . In fact, the diagonal correlators  χ B 2 and χ S 2 show a sizable input from σ 2 n (the variance of the neutron distribution), and both σ 2 K 0 and σ K 0 K , respectively. The case of χ Q 2 is peculiar since, as evident from Fig. 8, the nonmeasured contribution is almost negligible, when compared to the measured one.

C. Isospin randomization
Another important effect we have not addressed yet, which is present in experiment, is the isospin randomization [50,51]. This effect is caused by reactions that take place in the hadronic phase between nucleons and pions, and consist of the generation and decay of Δ resonances [Δð1232Þ prominently], through processes like with analogous ones for the antibaryons. For collision energies ffiffi ffi s p ≳ 10 GeV, the lifetime of the fireball is long enough to allow several of such cycles to take place, resulting in a complete randomization of the isospin of the nucleons [50,51]. This expectation has been confronted with data in [80] confirming a complete randomization with the exception of the highest energy: ffiffi ffi s p ¼ 200 GeV. For this paper, though, we assume complete randomization throughout. The distributions of protons and neutrons then factorize and the correlation between the two is erased. The average number of protons and neutrons, as well as antiprotons and antineutrons, and consequently the average net-proton and net-neutron number, are left unchanged by such reactions, but fluctuations are not. In particular, this results in an enhancement of both the net-proton and net-neutron variance, at the expense of the correlation between the two (note that the variance of net nucleon σ 2 N ¼ σ 2 p þ 2σ pn þ σ 2 n cannot be changed by these reactions). Similarly, charge conservation ensures that the sumQ ¼p þπ is conserved in the reactions in Eq. (28). It can be shown that this results in the net-pion variance σ 2 π being increased by the same amount as the net-proton variance σ 2 p . Since the sum σ 2 pþπ ¼ σ 2 p þ 2σ pπ þ σ 2 π must also be left unchanged, we have that σ pπ is decreased by the same amount again: σ pn .
Thus there is information lost through the process of randomization, and the original σ pn or σ 2 p cannot be reproduced from data individually. The experimental access to those correlators where either of these plays an important role is very difficult. This is the case e.g., for χ BQ 11 , which is completely dominated by σ 2 p .

V. PROXIES
The issue of having particles that cannot be detected poses the problem of a loss of conserved charges. Historically, the proxies for baryon number, electric charge, and strangeness have been the protons, the p, π, K combination, and the kaons themselves, respectively.
We have seen in Figs. 9 and 10 how the single hadronic, measured (self-) correlators relate to the fluctuations of conserved charges. We can then find a correspondence between fluctuations of conserved charges and measurable (and calculable) hadronic fluctuations.
Both in theory and experiment, it is customary to consider ratios of fluctuations, in order to eliminate, at least at leading order, the dependence on the system volume. For this reason, we focus on the ratios χ BS 11 =χ S 2 and χ QS 11 =χ S 2 , for which we construct proxies using solely fluctuations of (measured) hadrons. The underlying assumption when considering ratios is that the freeze-out   of all species involved occurs at the same time in the evolution of the system, hence at the same volume. Let us start considering the χ BS 11 correlator. One could expect that, having both kaons and protons in the bulk of particle production, their correlator σ pK would be a good proxy. However, as we can see in Fig. 9, this is clearly not the case, as the proton-kaon correlator gives a negligible contribution to χ BS 11 . On the contrary, the variance of the net-Lambda distribution σ 2 Λ represents a much more sizable contribution to the total correlator.
In the upper panel of Fig. 13 we show the HRG model results for the ratio χ BS 11 =χ S 2 at μ B ¼ 0 (black, thicker line). As already mentioned, from Figs. 9 and 10 we see how the leading contributions to the two correlators come from σ 2 Λ and σ 2 K , respectively. We can then construct a tentative proxy asC which is shown as a green, dashed line. We see that, although this quantity reproduces very well the full result at low temperatures-where the kaons dominate-it overshoots at higher temperatures, and in particular around the QCD transition and chemical freeze-out temperatures, which are obviously the interesting regime. It is worth noticing that, in order to construct a good proxy for a ratio of conserved charges fluctuations, it is not sufficient to choose the best proxy for both the numerator and the denominator. In fact, a good proxy for the ratio is obtained when the proxy in the numerator and the denominator are equally good. Some guidance in this construction is then provided by Fig. 10, where the extent to which a hadronic correlator reproduces the corresponding BQS fluctuation is most evident. For this reason, we consider adding the contribution from the net-Λ fluctuations to χ S 2 too, and definẽ which is shown as a blue, dotted line. We see how this second proxy is much better at reproducing the full result, as it is very close to it at all temperatures, including in the vicinity of the QCD transition. In addition, again referring to Figs. 9 and 10, it is interesting to try and include the contributions from multistrange hadrons, both in the numerator and denominator. With these, one has which is shown as the orange, dashed-dotted line, and also reproduces very well the behavior of the full ratio, although not really improving the situation over the previous one. As a final check, one can build a proxy from the σ pK correlator asC pK;ΛK which is shown as the yellow, dashed-double-dotted line. Not unexpectedly, this combination is not able to serve as a good proxy. The case of χ QS 11 =χ Q 2 follows directly from the previous one, and is shown in the lower panel of Fig. 13. In fact, in a system with 2 þ 1 quarks (with no isospin symmetry breaking) the following relation applies: from which one can derive that Thus, exploiting this relation and the good proxyC Λ;ΛK BS;SS we have defined for χ BS 11 =χ S 2 , we can define and have a proxy for χ QS 11 =χ S 2 for free, which indeed works very well over the whole temperature range. Now, we wish to consider the correlator χ BQ , which is the only one of the three nondiagonal correlators to be influenced by the isospin randomization discussed in the previous section. By looking at Figs. 9 and 10, it is natural to construct the proxỹ where the net charge is typically defined asQ ¼p þ π þK. In Fig. 14 the total contribution and this proxy are shown with a solid, thick black line and a dashed green line, respectively, and the agreement is extremely good. When including the effect of isospin randomization-which does not affect the denominator-the situation is radically different, and the corresponding curve is shown with a dotted blue line. The increase in the net-proton variance spoils the effectiveness of this proxy. A quantity which is not affected by this effect is the sum σ 2 p þ σ pπ , from which we can define the ratiõ This quantity is also shown in Fig. 14 as a dashed-dotted orange line, and clearly cannot serve as a good proxy. It is interesting to notice how the increase that σ 2 p receives by the isospin randomization is almost exactly equal to σ pπ . As we discussed in Sec. IV C the effect of isospin randomization on σ 2 p amounts to σ pn . The presently used HRGbased approach introduces a p − π (σ pπ ) and p − n (σ pn ) correlation through the decay of the same resonances (Δ).
From Fig. 14 we see that, because of this effect, it is not possible to build a suitable proxy for χ BQ 11 =χ Q 2 . For analogous reasons, it is not possible to create a good proxy for the ratio χ BQ 11 =χ B 2 . Having discussed all three combinations of the offdiagonal cross-correlators we are lacking a good proxy for a correlator ratio involving only the light quarks. As a detour from the main line of the discussion we show that this is also a difficult task in the case of the diagonal correlators. Consider the ratio χ B 2 =χ Q 2 . In Fig. 15 we see the temperature dependence of this ratio at μ B ¼ 0, and the behavior of some tentative proxies alongside it. We start by considering the quantitỹ where we take advantage of the fact that, after the isospin randomization, one has σ 2 Net-N ¼ 2σ 2 p . This quantity is shown in Fig. 15 as a green dashed line, and we see that its contribution is not sufficient. We then consider adding the contribution from Λ baryons, and show as a dotted blue line the quantitỹ which improves on the previous one, but is still not satisfactory. We then try removing the contribution from net kaons at the denominator-which we can do regardless of isospin randomization, C Net-NΛ;pπ and finally include the contribution from multistrange baryons in the numerator, These last two proxies are also shown in Fig 15 as an orange dashed-dotted and as a yellow-double-dotted line, respectively. We see that both compare relatively well with the total contribution, with the latter being the better one. It is quite interesting how difficult it was to construct a suitable proxy for light-quark-dominated observables, in comparison to the previous cases of χ BS 11 and χ QS 11 . This is mainly due to the fact that (i) net charge is such a good proxy for χ Q 2 that is hard to match for other correlators, and (ii) isospin randomization prevents one from building proxies with fluctuations of net proton only.
We have seen in this section how to construct good proxies for ratios including both diagonal and crosscorrelators of conserved charges. The proxies including strangeness make use of only a couple of hadronic observables, namely the variances σ 2 K and σ 2 Λ , more precisely, only their ratio. It is also remarkable how the addition of multistrange baryons to the proxy for χ BS 11 =χ S 2 is not necessary, as it does not improve the existing agreement. We also saw that for light-quark-dominated observables, isospin randomization modifies the correlators of net proton, net pion and net neutron, preventing the construction of useful proxies for such observables.

VI. FINITE CHEMICAL POTENTIAL AND KINEMATIC CUTS
Since experimental measurements for moments of netparticle distributions are currently available both from the LHC and RHIC, it is interesting to analyze the behavior of the quantities we are studying also at finite values of the baryon chemical potential. In the left panel of Fig. 16, we show the behavior of the proxies along parametrized chemical freeze-out lines-shifted in T from the parametrization in [81]-with T intersects at T 0 ¼ 145; 165 MeV, so to "bracket" the crossover region of QCD. The ratios χ BS 11 =χ S 2 and χ QS 11 =χ S 2 are shown in the first and second row, respectively. We see that for these ratios, the agreement with the considered proxies does not worsen with the increase in the chemical potential, and the curves remain very close for a broad range of collision energies. This means that the scope of the proxies we have constructed to reproduce the behavior of fluctuations of conserved charges is not limited to small μ B , but can be extended to the study in the Beam Energy Scan as well.
In Sec. IV we have mentioned that one of the strengths of the HRG model is the possibility it offers to include effects that are present in the experimental situation, like the use of cuts on the kinematics. In the central panel of Fig. 16 we show the same scenario as in the left panel, but with the inclusion of exemplary mock cuts: 0.2 ≤ p T ≤ 2.0 GeV, jyj ≤ 1.0. These cuts do not correspond to any past or ongoing measurement at the LHC or RHIC, but are constructed such as to be reproducible in the experiment, and still give a hint of the effect of including the cuts at all. For a systematic treatment of the dependence of fluctuations on the kinematic cuts-which is beyond the scope of this work-see [82], where it is studied in a thermal model with an older hadron list and without the inclusion of resonance decays. In our example, the same cuts are applied to all particle species. We see that for all the observables considered the agreement between net-charge fluctuation ratios and proxies remains the same as in the case without cuts, for both freeze-out lines.
Finally, in the right panel of Fig. 16 we show the selected proxies, for both freeze-out curves, comparing the cases with and without the cuts. We see that the effect is very minimal for the two ratios χ BS 11 =χ S 2 and χ QS 11 =χ S 2 . This is obviously of key importance in light of a potential direct comparison to results from lattice QCD calculations, as the one discussed in Sec. III for χ BS 11 =χ S 2 . This is one of the main reasons these proxies were built in the first place.
The third row in Fig. 16 shows the behavior of the ratio χ B 2 =χ Q 2 when acceptance cuts are introduced. As opposed to the discussed off-diagonal ratios it shows a large dependence on the cuts. Thus, even though this ratio does not suffer from the effect of isospin randomization, a comparison to lattice simulations can be problematic. Thus, we focus on the strangeness related off-diagonal correlators in the next section, where we compare to experimental data.

VII. COMPARISON TO EXPERIMENTAL RESULTS
In the previous section we have considered the impact of including kinematic cuts on the proxies we have defined previously, by considering some exemplary cuts which were chosen to be the same for all particle species. However, experimental measurements exist for different species, and it is possible to test how the proxies we constructed compare to the experimental results, this time including the corresponding cuts on a species-by-species (or measurement-by-measurement) basis.
In Fig. 17 we show the behavior of the proxiesC Λ;ΛK BS;SS andC K;ΛK QS;SS from Eqs. (30) and (35), along the same freezeout lines used in Fig. 16, and compare them to available experimental results from the STAR Collaboration [56,57]. The important difference is that now the experimental cuts are the ones taken from the actual measurements, and namely they are not the same for the different species.
We see that the proxy (it is only one independent quantity as discussed above) works well also in comparison with available experimental data, when the considered freeze-out line is the one with a temperature Tðμ B ¼ 0Þ ¼ 165 MeV. This is in line with results from other analyses, which indicate that strange particles seem to prefer a higher chemical freeze-out temperature [59].
One more remark is in order: by comparing, e.g., the curves in Figs a certain proxy. In fact, the same ratioC Λ;ΛK BS;SS is shown, with the difference that in Fig. 16 the same cuts are applied to both Λ and K, while in Fig. 17 the cuts utilized are those from the experimental analyses, namely 0.9 < p T < 2.0 GeV, jyj < 0.5 for net-Λ [57] and 0.4 < p T < 1.6 GeV, jyj < 0.5 for net kaon [56]. Due to this difference in the applied kinematic cuts, more than a factor 2 separates the two curves. For this reason, a direct comparison to lattice QCD would be premature.

VIII. CONCLUSIONS
In this work, we first presented new continuumextrapolated lattice QCD results for second-order nondiagonal correlators of conserved charges. While the continuum extrapolation is a straightforward task at μ B ¼ 0, results need to be extrapolated to the real μ B regime, which cannot be simulated directly. This is always ambiguous, so we compared two different schemes for the ratio χ BS 11 =χ S 2 and performed a continuum extrapolation in the regime where the two approaches agree.
We performed an HRG-model-based study on the secondorder correlators, both diagonal and nondiagonal. At μ B ¼ 0 we found agreement with lattice. Then we showed how they relate to fluctuations of those hadronic species which can be measured in heavy ion collision experiments. What percentage of these correlators is accounted for by particles that can actually be detected in the experiment varies quite considerably from observable to observable.
In order to compare either to lattice QCD results or experimental measurements, we focused on ratios of fluctuations, whose behavior can be reproduced through commonly measured hadronic observables, i.e., proxies.
In the following we summarize the findings for a ratio with each of the three possible cross-correlators of baryon (B), electric charge (Q) and strangeness (S).
The BQ correlator in equilibrium is dominated by proton fluctuations, with the other contributions-most notably the proton-pion correlation and hyperons self-correlationsalmost perfectly canceling each other. Nonetheless, the information loss caused by isospin randomization prevents one from constructing successful proxies for ratios including χ BQ 11 . Luckily, neither the isospin randomization nor the introduction of cuts on the kinematics had a significant effect on either χ BS 11 or χ QS 11 . Because of this, we were able to construct proxies for the ratios χ BS 11 =χ S 2 and χ QS 11 =χ S 2 that are within 10% of the grand canonical prediction. These two ratios are not independent, since in the isospin symmetric case they are related by the Gell-Mann-Nishijima formula. It is striking that only two measured quantities, namely the variances of net-kaon and net-Lambda distributions, were sufficient to build the proxies for χ BS 11 =χ S 2 and χ QS 11 =χ S 2 . We showed that the inclusion of multistrange hyperons does not improve the quality of the proxy. Moreover, although these are cross-correlators of conserved charges, particle species cross-correlators do not contribute significantly. In fact, none of the particle cross-correlators contributes to any of the charge fluctuations or crosscorrelators, with the exception of the proton-pion one, but the latter is largely affected by isospin randomization.
Thus, we have a ratio at hand that is available both from lattice simulations and for experimental measurement. The ratio χ BS 11 =χ S 2 behaves as a strangeness-related thermometer for chemical freeze-out. We provided continuum extrapolated results at 0 and finite chemical potential for this quantity.
Finally, we compare our results to experiment. A direct use of lattice data in the experimental context would require the use of the same kinematic cuts for Λ and K. The STAR Collaboration has published results for fluctuations of K and preliminary results for Λ fluctuations, though with different kinematic cuts. To test our proxy, we recalculated its HRG model prediction with the actual cuts used in experiment. We saw that the σ 2 Λ =ðσ 2 K þ σ 2 Λ Þ ratio quite evidently favors the higher freeze-out temperature, in line with what was already shown by other analyses [59,60].
These high temperatures for the chemical freeze-out motivates the use of lattice QCD in future studies since they fall at the limit of the validity of the HRG model.