1 Introduction

High-energy heavy-ion collision experiments at various beam energies \(\sqrt{s_{NN}}\) have enriched our understanding of the properties and phases of strongly interacting matter. The transient creation of a color-deconfined state in the laboratories was one of the major scientific successes in the last two decades. The deconfinement transition is an analytic crossover for vanishing baryon chemical potential \(\mu _B\) [2], where the transition region \(T_c=(154\pm 9)\) MeV [3, 4] is rather broad in temperature T. While these information base on first-principle lattice QCD calculations, another fascinating landmark in the phase diagram, the QCD critical point, has not yet been discovered with this method despite being predicted by various approaches [5, 6].

Information about the properties of hot and dense QCD matter can be inferred indirectly from the measured particle spectra and their event-by-event fluctuations. The success of statistical hadronization models in describing the particle production in heavy-ion collisions ranging from AGS to LHC beam energies [7,8,9,10] led to the conclusion that the produced hadronic matter originates from a source in or near thermal and chemical equilibrium. A common freeze-out curve [7] in the phase diagram could be drawn, highlighting the thermal conditions (\(T, \mu _B\)) at chemical freeze-out where the hadrochemistry is fixed.

Higher-order moments of the event-by-event particle multiplicity distributions provide an additional excellent measure to characterize the matter properties and to reveal the phase structure in QCD. While it is debated whether fluctuations originate from an equilibrated hadronic medium [11, 12] it is to a first extent reasonable to assume, that if this is true for the means also the lowest-order fluctuations, i.e. the variances, should be describable within statistical models. First studies [13, 14], utilizing experimental data on net-proton and net-electric charge fluctuations to determine the freeze-out conditions, found at large \(\sqrt{s_{NN}}\) significantly reduced freeze-out parameters compared to [7]. In the work [13], the analysis was quantitatively driven by the contributions from protons and anti-protons as well as charged pions.

In a recent study [1], the experimental data [15] on lowest-order fluctuations in the net-kaon number \(N_{K^+}-N_{K^-}\) as function of \(\sqrt{s_{NN}}\) were used to extract freeze-out conditions. For a unique determination, information from lattice QCD for the isentropic trajectories [16] running through the freeze-out points of [13] was supplemented. The obtained freeze-out temperatures were reported to be substantially larger than those found in [13]. In our work, we combine an analysis of the net-kaon fluctuations [15] and the yield ratios of strange anti-baryons over baryons [17,18,19] as functions of \(\sqrt{s_{NN}}\). This complementary approach allows us to determine whether strange baryon yields and net-kaon fluctuations can be described with the same T and \(\mu _B\).

Our strategy is to analyze observables that are highly sensitive to variations in the thermal parameters. As we discussed in [20], the ratio of variance over mean of the net-kaon number is such an observable that varies rapidly with T in a statistical model. The same is true for the heavier baryons while we use the lighter baryons as a baryometer in this work. Therefore, we consider the combination of observables studied here as optimized for determining freeze-out conditions. One should keep in mind, nevertheless, that the net-kaon number is not a conserved charge in QCD and as such prone to late stage processes like resonance decays.

2 Theoretical framework

We perform our analysis of net-kaon fluctuations and strange anti-baryon to baryon yield ratios using a Hadron Resonance Gas (HRG) model for a grandcanonical ensemble of non-interacting hadrons and resonances. Such a model, in which the strong interaction between hadrons is effectively accounted for by contributions from resonances [21], was shown to reliably describe basic thermodynamic quantities [22,23,24] as well as susceptibilities and their ratios [25,26,27] from lattice QCD. The hadrons and resonances are considered to be point-like in our work. The pressure P in this framework reads

$$\begin{aligned} P = \sum _i (-1)^{B_i+1} \frac{d_iT}{(2\pi )^3} \!\int \!\text {d}^3k \ln \left[ 1+(-1)^{B_i+1}z_i e^{-\epsilon _i/T}\right] , \end{aligned}$$
(1)

where the sum runs over all particle species included in the framework. In Eq. (1) the particle energy reads \(\epsilon _i=\sqrt{k^2+m_i^2}\) for momentum k and particle mass \(m_i\), \(d_i\) is the degeneracy factor and \(z_i=e^{\mu _i/T}\) is the fugacity. The particle chemical potential \(\mu _i\) is defined as \(\mu _i=B_i\mu _B+S_i\mu _S+Q_i\mu _Q\), where \(\mu _X\) denotes the chemical potential of the conserved charge X and \(X_i=B_i, S_i, Q_i\) represent the quantum numbers of the conserved baryon, strangeness and electric charge, respectively.

In Eq. (1), we use an updated version of the spectrum of hadrons and resonances in line with the recent listing from the Particle Data Group [28]. In previous studies [29, 30], the advantages of using such an update were discussed in detail. Moreover, in [31] the influence of unconfirmed resonances in the 2016 listing [28] versus only confirmed resonances on determined freeze-out parameters was investigated.

The net-density \(n_X\) of a conserved charge is given by \(n_X=\sum _i X_i n_i\), where the individual particle densities \(n_i\) follow from derivatives, \(n_i=(\partial P/\partial \mu _i)\), at fixed T. With this, physical conditions met in a heavy-ion collision experiment can be implemented into the model [32] by requiring that \(n_S=0\) and \(n_Q=x n_B\). These account for net-strangeness neutrality in the fireball and an initial proton to baryon ratio at mid-rapidity which for Au\(+\)Au and Pb\(+\)Pb collisions is approximately \(x\simeq 0.4\). As a consequence, the chemical potentials \(\mu _S\) and \(\mu _Q\) become functions of T and \(\mu _B\). Due to the lack of stopping at high \(\sqrt{s_{NN}}\), the mid-rapidity region is almost isospin symmetric. As \(\mu _B\), and thus \(n_B\), is small for high beam energies this is also approximately satisfied by the second physical condition.

The experimentally realized phase-space coverage, which can be limited in rapidity y, transverse momentum \(k_T\) and azimuthal angle \(\phi \) due to the detector design and demands from the analysis, can be respected in a straightforward way. Following [33], kinematic acceptance cuts are implemented into the model by restricting the momentum integrals in Eq. (1) accordingly. This requires replacing the integration measure \(\text {d}^3k\) by \(k_T\sqrt{k_T^2+m_i^2}\cosh (y)\,\text {d}k_T\,\text {d}y\,\text {d}\phi \) and \(\epsilon _i\) by \(\cosh (y)\sqrt{k_T^2+m_i^2}\). It should be noted that this procedure does not take into account the elastic scatterings between chemical and kinetic freeze-out, which can transport individual particles in and out of the experimental acceptance. In order to take this correctly into account, a fully dynamical transport approach for the hadronic phase would need to be applied. We think, however, that only a small number of all considered particles are actually affected by this kind of final state effect.

Table 1 Summary of considered rapidity windows in the analysis of \(\bar{B}/B\)-ratios. The beam energy \(\sqrt{s_{NN}}\) is given in GeV
Full size table

The effect of resonance decays can be implemented into the framework in an explicit way which allows us to keep correctly track of the strangeness transfer from mother to daughter including resonances such as e.g. \(\Xi (1690)^-\) or N(1650). This is different from the way resonance decays were implemented in [1]. Our method was developed in [34, 35] and applied to study the impact of resonance decays on net-proton fluctuations without [36] and in the presence of a QCD critical point [37]. The final particle number \(N_j\) of a stable, i.e. with respect to strong and electromagnetic decays, hadron is given by the sum \(N_j=N_j^*+\sum _R \langle N_j^R\rangle _R\) of primordially, i.e. directly, produced hadrons \(N_j^*\) and the contributions stemming from resonance decays. Those produce on average over the decays \(\langle N_j^R\rangle _R=N_R^*\langle n_j\rangle _R\) hadrons of type j associated with the branching ratios \(b_r^R\) of resonance R via \(\langle n_j\rangle _R=\sum _r b_r^R n_{j,r}^R\) for integer \(n_{j,r}^R\). On an event-by-event basis, the numbers \(N_j^*\) and \(N_R^*\) fluctuate thermally. But, in addition, the actual number of hadrons of type j originating from the decay of R follows a multinomial probability distribution [34]. Thus, fluctuations in the contributions from resonance decays are caused both by thermal fluctuations in \(N_R^*\) and by the probabilistic character of the decay process. For the mean, the latter has no consequences and thus the mean of the final particle number after resonance decays is given by

$$\begin{aligned} M_j = \langle N_j^*\rangle _T + \sum _R \langle N_R^*\rangle _T \langle n_j\rangle _R \,. \end{aligned}$$
(2)

Accordingly, for the net-kaon number, the mean is given by \(M_K=M_{K^+}-M_{K^-}\). The variance of the net-kaon number is instead influenced by the probabilistic nature of the decay and follows as [34, 35]

$$\begin{aligned} \sigma _K^2 =&\,\langle (\Delta N_{K^+}^*)^2\rangle _T + \langle (\Delta N_{K^-}^*)^2\rangle _T \nonumber \\&+ \sum _R \langle (\Delta N_{R}^*)^2\rangle _T \left( \langle n_{K^+}\rangle _R^2 + \langle n_{K^-}\rangle _R^2 \right) \nonumber \\&- 2 \sum _R \langle (\Delta N_{R}^*)^2\rangle _T \langle n_{K^+}\rangle _R \, \langle n_{K^-}\rangle _R \nonumber \\&+ \sum _R \langle N_{R}^*\rangle _T \left( \langle (\Delta n_{K^+})^2\rangle _R + \langle (\Delta n_{K^-})^2\rangle _R\right) \nonumber \\&- 2 \sum _R \langle N_{R}^*\rangle _T \langle \Delta n_{K^+} \Delta n_{K^-}\rangle _R \,, \end{aligned}$$
(3)

where \(\Delta N_i=N_i-\langle N_i\rangle _T\) and \(\langle \Delta n_i\Delta n_l\rangle _R=\langle n_i n_l\rangle _R - \langle n_i\rangle _R \langle n_l\rangle _R = \sum _r b_r^R n_{i,r}^R n_{l,r}^R - \langle n_i\rangle _R \langle n_l\rangle _R\). As apparent from Eq. (3), resonance decays introduce correlations between the number of \(K^+\) and \(K^-\). Neglecting the probabilistic nature of the decay processes, the variance reads instead

$$\begin{aligned} \sigma _K^2 =&\,\langle (\Delta N_{K^+}^*)^2\rangle _T + \langle (\Delta N_{K^-}^*)^2\rangle _T \nonumber \\&+ \sum _R \langle (\Delta N_{R}^*)^2\rangle _T \left( \langle n_{K^+}\rangle _R^2 + \langle n_{K^-}\rangle _R^2 \right) \nonumber \\&- 2 \sum _R \langle (\Delta N_{R}^*)^2\rangle _T \langle n_{K^+}\rangle _R \, \langle n_{K^-}\rangle _R \,. \end{aligned}$$
(4)

In Eqs. (2)–(4) only the thermal averages \(\langle \cdot \rangle _T\) can be obtained from the HRG model pressure in Eq. (1) as the derivatives

$$\begin{aligned} \langle N_i\rangle _T =&\, VT^3\, \frac{\partial (P/T^4)}{\partial (\mu _i/T)} \equiv V n_i \,, \end{aligned}$$
(5)
$$\begin{aligned} \langle (\Delta N_i)^2\rangle _T =&\, VT^3\, \frac{\partial ^2(P/T^4)}{\partial (\mu _i/T)^2} \,, \end{aligned}$$
(6)

where V is the volume which cancels in ratios.

3 Analysis and results

We determine the conditions for T and \(\mu _X\) at chemical freeze-out by applying the framework outlined above to optimally describe experimental data from RHIC on strangeness observables measured by the STAR Collaboration. We analyze data on yield ratios of strange anti-baryons over baryons, \(\bar{B}/B\), as well as lowest-order net-kaon fluctuations, \(\sigma _K^2/M_K\), simultaneously.

For the \(\bar{B}/B\)-ratios we study the published results in [17, 18] for \(\sqrt{s_{NN}}=200\) and 62.4 GeV, and for smaller \(\sqrt{s_{NN}}\) the preliminary results from the RHIC Beam Energy Scan reported in [19]. The ratios for most central collisions, shown as functions of \(\sqrt{s_{NN}}\) in Fig. 1 (upper panel) for \(\Lambda \) (squares), \(\Xi ^-\) (triangles) and \(\Omega ^-\) (circles), were determined from \(\phi \)- and \(k_T\)-integrated yields in a given rapidity window around mid-rapidity. In line with the experimental set-up, we consider the y-ranges summarized in Table 1.

For the net-kaon fluctuations we take the results recently reported in [15] which are corrected for finite detector efficiency and the centrality bin width effect. The data for most central collisions are shown by the squares in Fig. 1 (lower panel) where we include only the dominating systematic error bars in the plot. For \(\sigma _K^2/M_K\), the experimental kinematic acceptance limitations are 0.2 GeV/c \(\le k_T \le 1.6\) GeV/c and \(-0.5\le y\le 0.5\) with full azimuthal coverage.

In the analysis of the data, the contributions from resonance decays play an important role. For example, for the net-kaon fluctuations, significant correlations between \(K^+\) and \(K^-\) are induced by the decay processes. For the ratio \(\bar{\Lambda }/\Lambda \) we include in addition to the final particle numbers of \(\Lambda \) and \(\bar{\Lambda }\) also the contributions from electromagnetic \(\Sigma ^0\) and \(\bar{\Sigma }^0\) decays. Those contributions were not corrected in the experimental analysis [17,18,19].

Fig. 1
figure 1

Upper panel: Ratio \(\bar{B}/B\) of strange anti-baryon over baryon yields for \(\Lambda \), \(\Xi ^-\) and \(\Omega ^-\) as a function of the beam energy \(\sqrt{s_{NN}}\). The published and preliminary STAR data (symbols) are taken from [17,18,19], respectively. The colored bands show the fit results for \(\bar{B}/B\) within the employed HRG model for the freeze-out conditions determined in this work. Lower panel: Lowest-order net-kaon fluctuation measure \(\sigma _K^2/M_K\) of variance \(\sigma _K^2\) over mean \(M_K\) as function of \(\sqrt{s_{NN}}\). The published STAR data (squares) are taken from [15], where the shown error bars indicate only the dominating systematic error. The colored band depicts our fit results for the freeze-out conditions shown in Fig. 2

Full size image
Fig. 2
figure 2

Chemical freeze-out conditions for T and \(\mu _B\) determined for the beam energies \(\sqrt{s_{NN}}=200, 62.4, 39, 27, 19.6, 11.5\) GeV (from left to right). The red, solid squares show the results determined in this work from the combined analysis of lowest-order net-kaon fluctuations, \(\sigma _K^2/M_K\), and \(\bar{B}/B\) yield ratios. The error bars correspond to the bands in the fits shown in Fig. 1. The blue, open squares depict the freeze-out conditions reported in [13] which were determined from a combined analysis of lowest-order net-proton and net-electric charge fluctuations

Full size image

The optimal fits for \(\bar{B}/B\), where the means are obtained via Eq. (2), and \(\sigma _K^2/M_K\) including the reported error bars are shown by the colored bands in Fig. 1. The corresponding freeze-out conditions for T and \(\mu _B\) are shown in Fig. 2 by the red, solid squares. For comparison, we contrast the freeze-out conditions [13] deduced from an analysis of lowest-order net-proton and net-electric charge fluctuations which are shown by the blue, open squares. Similar to [1], we observe in particular for large \(\sqrt{s_{NN}}\) a significant enhancement of the chemical freeze-out temperature compared to [13] while for smaller \(\sqrt{s_{NN}}\) the two different results approach each other. Moreover, as in [38], which bases its analysis entirely on hadronic yields distinguishing strange from non-strange hadrons, we find a non-negligible modification of \(\mu _B\) compared to [13]. Our results and the ones reported in [1] overlap within errors, but our results are systematically smaller in T and larger in \(\mu _B\) compared to [1]. The determination of the freeze-out temperature is sensitively influenced by the net-kaon fluctuation data as was already discussed in [20]. Here, our use of the heavier (anti-)baryons adds sensitivity to T, while the lighter (anti-)baryons influence stronger the determination of \(\mu _B\).

The electric charge chemical potential \(\mu _Q\) is negative and negligibly small compared to \(\mu _B\). In contrast, the strangeness chemical potential \(\mu _S\) constitutes a non-negligible fraction of \(\mu _B\) according to the condition of strangeness neutrality as was stressed in previous lattice QCD studies [25, 39]. In Fig. 3, we depict our results for the ratio \(\mu _S/\mu _B\) (red, open squares) for the freeze-out conditions shown in Fig. 2. The relatively large errors in \(\mu _S/\mu _B\) are a consequence of the errors in \(\mu _B\) and in T and the imposed strangeness neutrality. For comparison, we also show the lattice QCD result \(\mu _S/\mu _B=s_1(T)+s_3(T)\mu _B^2\) (colored band) evaluated for the central values of T and \(\mu _B\) determined in this work using the data on \(s_1(T)\) and \(s_3(T)/s_1(T)\) published in [25]. The bandwidth results from the error bars reported there.

Fig. 3
figure 3

Ratio \(\mu _S/\mu _B\) of strangeness over baryon chemical potential as a function of \(\sqrt{s_{NN}}\) (red, open squares) for the freeze-out conditions shown in Fig. 2. The error bars account for the bands in the fits shown in Fig. 1. Using our determined central values for T and \(\mu _B\) we show, in comparison, the corresponding lattice QCD ratio \(\mu _S/\mu _B=s_1(T)+s_3(T)\mu _B^2\) (colored band), where the bandwidth is deduced from the data on \(s_1(T)\) and \(s_3(T)/s_1(T)\) and their errors reported in [25]

Full size image

Finally, the importance of the fluctuation contributions due to the probabilistic character of resonance decays can be tested explicitly in our framework. Leaving these contributions out, the variance \(\sigma _K^2\) is given by Eq. (4), while the mean remains unaffected. Performing the analysis of the experimental data again with this set-up, we find a reduction of the determined freeze-out temperature by about 5%, while the obtained \(\mu _B\) is up to 18% underestimated. This highlights that, indeed, the probabilistic nature of the decay processes influences our final results in a non-negligible way.

4 Conclusions and outlook

In this work we determined chemical freeze-out conditions by analyzing strangeness observables from RHIC measurements at different beam energies. In our analysis we studied both the lowest-order net-kaon fluctuations in [15] and strange anti-baryon over baryon yield ratios from [17,18,19] within a Hadron Resonance Gas model framework. From the combined optimized fit, the freeze-out temperature as well as the chemical potentials associated with the conserved charges of QCD were inferred. We find that the obtained freeze-out temperature is significantly enhanced at large \(\sqrt{s_{NN}}\) compared to the results of a previous study [13] which used data on net-proton and net-electric charge fluctuations instead. The baryon chemical potential is also visibly enhanced compared to [13] except for the smallest \(\sqrt{s_{NN}}\). For smaller \(\sqrt{s_{NN}}\), the freeze-out conditions from both approaches start to converge. We note that the results for T and \(\mu _B\) determined in this work are compatible within errors with those from the statistical hadronization model fits of hadron yields (or ratios) [7,8,9,10].

The results presented here are in qualitative agreement with the findings of a recent study [1] which used the net-kaon fluctuation data supplemented by information from lattice QCD. Our complementary study shows that also measured yield ratios of strange anti-baryons over baryons can be described by similar freeze-out conditions. For a precise determination of the latter, the correct implementation of resonance decay contributions plays an important role. We find that resonance decays lead, for example, to significant correlations between \(K^+\) and \(K^-\) reducing the lowest-order fluctuation measure \(\sigma _K^2/M_K\) of the net-kaon number by about 15% from its Skellam limit for the same thermal parameters.

The freeze-out conditions presented in Figs. 2 and 3 are, of course, subject to both the quality of the analyzed data and limitations in the applied framework. The data on \(\bar{B}/B\)-ratios from the Beam Energy Scan [19] are still preliminary and therefore quantitative changes in the determined freeze-out conditions at these \(\sqrt{s_{NN}}\) can be expected for future published data, especially in terms of the considered error bars. Nevertheless, we expect these changes to be small in our combined analysis of net-kaon fluctuations and \(\bar{B}/B\)-ratios.

The shown errors in our results base entirely on the errors reported for the analyzed data. Additional uncertainties, stemming from limitations in our framework, have not been included. Different aspects that could lead to a refined analysis are, for example:

(i) If the pion density in the fireball is large and the duration of the late hadronic stage long enough, the regeneration and subsequent decay of \(K^*\) resonances provides an additional source of fluctuations in the net-kaon number. The impact of even only a partial isospin randomization [40, 41] is to bring a distribution closer to the corresponding Skellam limit as we discussed in [36] in the case of net-proton fluctuations. To describe the data [15] on net-kaon fluctuations in this case would imply an enhancement of the chemical freeze-out conditions. Quantitative predictions of this effect would depend on the made assumptions.

(ii) Exact global charge conservation realized on an event-by-event basis can cause large effects on the fluctuations. This effect was discussed, for example, in [42,43,44,45] suggesting a future study in a canonical ensemble formulation.

(iii) Fluctuations in the number of participants for a given centrality class can sizeably influence the measured net-particle number fluctuations [44]. This effect is stronger for higher-order cumulants and smaller beam energies, and is more pronounced for central than very peripheral collisions. Its inclusion into the theoretical framework could potentially affect the results on freeze-out conditions deduced from fluctuation observables in a refined analysis.

Our work provides a baseline which neglects dynamical as well as critical fluctuation effects on the net-kaon number fluctuations. First estimates for the impact of critical fluctuations on these will be reported elsewhere.