Systematics of yields of strange hadrons from heavy-ion collisions around threshold energies

The parametrizations of experimental yields of K±,0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\pm ,0}$$\end{document}, ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} and Λ+Σ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda +\Sigma ^0$$\end{document} are proposed as function of available energy, sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\textrm{NN}}}$$\end{document} , and number of participants, ⟨Apart⟩b\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle A_{\textrm{part}} \rangle _{\textrm{b}}$$\end{document} , for sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\textrm{NN}}}$$\end{document} from 2.15 to 3 GeV. For all the dataset the ⟨Apart⟩b\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle A_{\textrm{part}} \rangle _{\textrm{b}}$$\end{document} was extracted using the Glauber Monte Carlo method. The α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} exponent of yield dependency on ⟨Apart⟩b\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle A_{\textrm{part}} \rangle _{\textrm{b}}$$\end{document} appears not to change with beam energy and is found to be 1.30 ± 0.02. Our parametrization and the predictions of public versions of RQMD.RMF, SMASH and UrQMD transport models are compared to the HADES experimental data for Ar+KCl at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\textrm{NN}}}$$\end{document} of 2.61 GeV. The phenomenological parametrization currently offers the best overall description of these yields. Predictions are given for yields from Ag + Ag collisions at available energies of 2.41 and 2.55 GeV, analysed by HADES, Au + Au experiment at 2.16 and 2.24 GeV planned by this collaboration, some yields for STAR’s Au + Au collisions at 3 GeV, and for Au + Au collisions planned by CBM, up to 3.85 GeV.


Introduction
Three decades after the first systematics of meson production at beam energies around their thresholds in the free nucleon-nucleon (NN) channel [1], a considerable sample of yields of hadrons with strange valence quarks was obtained with help of the KaoS, FOPI and HADES setups installed at the SIS18 accelerator in the GSI.STAR recently explored this region with measurements at the available energy in the NN channel of √ s NN = 3 GeV [2]1 .For HADES, the analysis of strangeness from Ag+Ag at available energies of 2.41 and 2.55 GeV is ongoing [3] and the new experiment at energies between 1.95 and 2.25 GeV is planned.Also, the CBM being constructed together with SIS100 accelerator is expected to deliver soon new strangeness data from heavy-ion collisions at a few GeV [4].With this in mind, it is a good time to compile the state-of-art systematics of yields (P ) of hadrons with strange quarks at energies around thresholds for strangeness production.In this paper we present such systematics for K +,−,0 , φ and Λ + Σ 0 .However, we omit Ξ − and strange resonances due to currently too limited data to provide the systematics.As for nearly all the heavy-ion experiments the Σ 0 decays into Λ were indistinguishable from Λ produced directly in the collision zone, the yields of both hyperons were presented together.Therefore, in our paper we follow this convention.
The values of threshold available energy needed for the production of the considered hadrons in a free NN collision are: 2.55 GeV for K + , K 0 and Λ (within the NN → NKΛ channel), 2.86 GeV for K − (produced in the NN → NNK + K − process) and 2.90 GeV for φ meson (for the NN → NNφ channel).A sub-threshold production of these hadrons in the heavy-ion collisions is possible due to opening up of many alternatives to direct NN collisions.For example, calculations within the transport models in the considered region of energies point to πY and BY (B = baryon, Y = hyperon) -initiated channels [5,6] as possibly the leading sources of K − mesons.The experiments together with some of these models reveal also a considerable feeddown from φ meson decays [7,8].In case of the latter hadrons, a broad mixture of channels may be responsible for their production, possibly including feeddowns from decays of N * and ∆ resonances [5,7,8].
The obvious quantities on which the yields of considered hadrons should depend are √ s NN and the average number of participant nucleons in the initial collision stage, A part b , the latter being the model-dependent function of the event centrality.Such a systematics could be a handy tool in estimating and/or comparing the strangeness production for any analysed or future heavy-ion experiment at comparable energies.Also, it could be helpful for theorists, as any model postulating a mixture of channels leading to production (and absorption) of strangeness in this domain, should reproduce the measured dependency.In our considerations we shall not go much beyond √ s NN ∼ 3 GeV, as the particle production mechanism is thought to pass from hadronlike to soft string-like nature.Mixing data from both regions to obtain one systematics may lead to biased conclusions.The dependency of the yield on A part b is universally found to have a power character [9][10][11].With an accumulated set of yields it is interesting to check whether or not the extracted exponent is the same for all the investigated hadrons, and whether α changes or not with √ s NN .An attempt to directly set together the A part b values from all the published data encounters the problem of different modelling of the collision centrality by various authors.For all the FOPI and part of the KaoS analyses the Geometrical model [12] was used, wherein the nuclear density profile was coarsely approximated by the sharp sphere.Within some KaoS works, the optical Glauber approach [13] was employed [10], however, with only one parametrization of the nuclear shape ("two-parameter Fermi", 2pF).For the recent analyses of HADES [14] and STAR [2] the Glauber Monte Carlo simulations [15] were performed.Therefore, the first task of our analysis is adjustment of the A part b estimations to an equal footing by consistent application of the Glauber Monte Carlo approach to all the data.While we certainly don't claim that our approach is final, we perceive it as a step forward compared to the geometrical sharp-sphere approach and optical Monte Carlo with just 2pF profile with fixed parameters.
The paper is organized as follows: in Sect. 2 the collected 4π yields are presented with emphasis on the centrality models used in the original papers.We then describe our procedure of A part b determination within the Glauber Monte Carlo approach.In Sect. 3 the results of the best fits of the parametrization curve to the yields as function of √ s NN and A part b are presented.The predictions of the curve for considered hadrons are compared to the range of known yields, and shown for some colliding systems where data is being analysed or planned to be measured.We also discuss the behaviour of α exponent of the dependency of yield on A part b .In Sect. 4 we compare the predictive power of our parametrization with results of some publicly available transport codes for the collisions of Ar+KCl at √ s NN = 2.61 GeV.We also present the predictions of our parametrization together with results of transport codes for central collisions of Ag+Ag at √ s NN of 2.41 and 2.55 GeV.Section 5 summarizes the paper.Some additional information is given in Appendices.

Data selection and A part b modelling
For the sample of strange hadron yields subjected to our analysis we accepted all the published 4π yields of K ±,0 and φ mesons, as well as Λ + Σ 0 baryons emitted from the nucleus-nucleus collisions at √ s NN = 2.16 -4.86 GeV, where the centralities or multiplicities of charged particles (MUL) were given.We also extracted the yields of K ± from central collisions of Ni+Ni at 2.67 GeV from their rapidity distributions, as the acceptance was wide and the data profile sufficient enough for the systematic uncertainty of the yield to be rather small [16].We did not include the yield of K + emitted from Ar+Ti at very low available energy of 1.92 MeV due to narrow acceptance, large relative uncertainty and no control of centrality [17].In total we collected 107 data points: 44 for K + , 30 for K − , 11 for K 0 s , 13 for Λ and 9 for φ.They are listed in Table 1 for K + , Table 2 for K − and φ and Table 3 for K 0 s and Λ.One can also find there the published centralities (based on cuts on MUL), the A part b values stated in the respective papers, and the types of models originally used by the Authors.Apparently the Geometrical Model dominates with 61 points (all the FOPI and AGS data, 60% of KaoS points, HADES data for Ar+KCl at √ s NN = 2.61 GeV and some others).The Glauber Optical approach was used 19 times, exclusively by KaoS.The Monte Carlo version of the Glauber model was applied 22 times (Au+Au data from HADES at 2.41 GeV and STAR data for 3 GeV).For the remaining 5 points only the centralities were given, not the A part b .Merging data with A part b obtained by different models weakens control over the field.In the prevailing Geometrical Model the nuclear shape profile was assumed to be of the sharp ball, which cannot be justified in light of measurements of the nuclear density profiles [18] and deformation models [19].A need to correct the situation seems to be evident.Therefore, we performed an unified procedure of extraction of A part b for all the data points, within the Glauber Monte Carlo approach.To achieve this we used the TGlauberMC code [20][21][22].The procedure performed for each colliding system at specific centrality range starts with a generation of 4•10 5 collision events, where for each event a set of positions of nucleons within nuclei is randomly generated according to an assumed density profile.A nucleon from a nucleus A is counted as participant if it overlaps with at least one nucleon from a nucleus B. The overlap is defined such that the centers of these nucleons are closer than σ NN,inel /π, where σ NN,inel is the cross section for the inelastic NN collision.This procedure gives a distribution of number of participants (N part ).For each experimental data point, a given centrality range is mapped into the range of N part .An averaged N part within this group is identified as A part b .Two kinds of input to this procedure are necessary: • value of σ NN,inel at a given √ s NN , • nuclear density profile.
For each experimental data point of strangeness yield the σ NN,inel was obtained as an isospin-averaged combination of the cross sections of pp (= nn) and pn interaction, where p and t subscripts correspond to projectile and target, and Z, N and A are the number of protons, neutrons and nucleons.The excitation functions of both the  3 Yields of K 0 and Λ emitted from nuclei colliding at given available energies (as well as beam kinetic energies) and centrality ranges.Values of Apart b are shown as published originally and according to our Glauber Monte Carlo analysis (see text for details.) pp and pn (np) channels were the result of the σ total − σ elastic subtraction of the data from [30].Results of such a subtraction were already shown in Figs.77 and 78 of [31].However, it has to be noted that the experimentally extracted values of np cross sections differ a bit from those from pn experiments.From this database, the values of elementary cross sections were interpolated at given √ s NN via the polynomial fit to the neighbouring points.We found that a precise form of interpolating function and statistical errors had a minor impact on the uncertainties.The dominant factors were: (1) the experiment-driven differences between σ pn and σ np , and (2) a relative scarcity of σ pn data points at √ s 2.25 GeV, as well as of σ np data points at √ s 2.3 GeV.These scarcities cause the susceptibility of the extracted values of the interpolation to the range of fitting window.To assess this contribution to the systematic errors, the highest and lowest variants of the extracted cross sections were found for each experiment.
An experimental extraction of the nuclear profile usually depends on the assumed functional form.The systematics [18,19] and other analyses [20][21][22][33][34][35][36] provide a broad range of variants of profiles for many nuclei.To assess the systematic error contribution to A part b due to choice of a nuclear shape we performed our Glauber MC calculations for the available functional forms and -for each form -the dominant uncertainties of values of relevant parameters.The details of this procedure are described in Appendix A. The values of A part b obtained within this procedure are shown in Tables 1, 2 and 3.The statistical uncertainties are negligible.The systematic ones stem from previously discussed uncertainties of σ NN,inel and of nuclear density profiles.As an uncertainty of each point we took half of distance between the maximal and minimal result due to variations of these contributions.
To shed light on the extent of changes that our method introduces, the ratios of corrected to originally published A part b values for K + data are shown in Fig. 1, where the numbering of data points is the one from the Table 1.The subsequent panels on this figure correspond to three different models estimating A part b , used in the original works.One observes that corrections are largest for the data obtained originally from the geometrical model.They often exceed 10% and in some cases 25%.
For the data where the Optical Glauber approach was originally used, the deviation is usually smaller, but still in one case reaches 20%.For the data where A part b was originally found using the Glauber MC approach, we are in good agreement with the original results, which supports our calculations.We have also verified that the RMS values of the density profiles obtained in our analysis are in agreement with the experimental ones, reported in [37].

Parametrization of yields of strange hadrons
Our aim is to find the type of curve that fits best to the available experimental yields of K ±,0 , φ and Λ at √ s NN between 2.16 and 3.0 GeV.As we show further, the AGS yields at 4.86 GeV deviate already too far from the rest of considered points.The curve should serve as a convenient tool to predict the yield at given √ s NN and A part b , providing also its uncertainty.To keep the numerical stability of these results, such a curve should have as few parameters, as possible.We selected the function type which factorizes the mass scaling as the power law and assumes the exponential dependence on available energy: We found the β parameters negative for every studied hadron type, therefore the curve has a saturating behaviour at large √ s NN .Despite this feature, the formula appears to fit better to the data than the power function, P ∼ √ s NN β .We found that for K + and K − the fitting procedure with all the four parameters left free generally worked well and the predictions of specific yields delivered rather reasonable results with limited uncertainty, which will be discussed further.However, the meager statistics in cases of the other hadrons caused problems with the covariance matrices.Therefore, for K 0 s , Λ + Σ 0 and φ we limited the parameter space by fixing C first, and adjusting its value manually, following the requirement of best χ 2 /ν of the 3-parameter fit.We checked the effect of this procedure for K + and K − , where both approaches were possible.We found that the values of predicted yields are identical up to third significant digit, but the fit in 3-parameter space generated uncertainties about 3 times lower than the full fit.It suggests that uncertainties for yields of K 0 s , Λ + Σ 0 and φ are probably also underestimated.We admit that it is a temporary solution in light of the limited statistics, and welcome any new data points to elevate the procedure to the full approach.
The parameters of the best fits of Eq. 2 to the experimental yields of considered hadrons are shown in Table 4.The full covariance matrices are presented in Appendix B. The χ 2 /ν values appear to be reasonable.For the case of K + , where χ 2 /ν reaches a relatively high value of 3.6, we note that the removal of point No. 34 in Table 1 results in a considerable drop of χ 2 /ν down to 2.1.However, we cannot pinpoint an experimental defect to substantiate the removal of this point.Interestingly, the α exponent of the yield dependence on A part b appears to be identical within uncertainties for K + and K − .This result confirms and sharpens the finding of a constant K − /K + ratio as a function of A part b , obtained for the smaller dataset of Au+Au at √ s NN = 2.52 GeV and Ni+Ni at 2.52 and 2.67 GeV [10].A relative goodness of fitting the Eq. 2 to the data also suggests, that the constancy of α exponent across the studied range of √ s NN may be a good hypothesis.We shall examine this issue further in the paper.A strength of the β parameter steers the steepness of the curve with √ s NN .We note that for Λ + Σ 0 the best fit value corresponds to an unusually high rise, and causes that for the experimental points at highest energy the curve is nearly saturated.Clearly more data is needed to verify this parametrization, and in particular the extrapolations of this curve seem unreliable, as it will be shown further.
For K + and K − the N , β and C parameters exhibit strong (anti-)correlations (c.f.Tab. 8).As a result, despite the fact that N parameters are burdened with large uncertainties, the non-diagonal terms of the covariance matrix strongly reduce the uncertainties of the estimated yields.
Fig. 2 visualizes the experimental data points and the best fit curves describing the yields as a function of √ s NN and A part b for two cases of, respectively, the richest statistics (K + ) and the poorest one (φ).A rather good factorization into these two quantities can be observed.For orientation, the standard deviations of experimental points from their fit predictions are shown in the lower panels of this figure.We shall now move to the predictions of yields with the found curves.The estimations of uncertainties for a given hadron type can be usually done via the error propagation formula including the covariance terms.We cross-checked this method by the bootstrap-like random pulling of curve parameters from the multivariate Gaussian distribution parametrized by the covariance matrix obtained from the fit.We found that the probability distributions of some of yields had non-Gaussian, asymmetric tails.Therefore, we modified the definition of the yield uncertainty to be: half of distance keeping central 68.3% of the probability distribution generated in such a way.In some cases we found the uncertainties up to 2 times larger than those calculated with the error propagation formula.All the uncertainties reported in this work were obtained using the bootstrap method.
Our first aim is to compare with the experimental data the interpolations at the lowest available energies (2.16 GeV for K + and 2.41 GeV for the other hadrons) and at the highest one (3 GeV) used for fitting.As it is seen in the upper panels of Table 5, the agreement is robust.We also compare the interpolations to the experimental data for 35% most central Ar+KCl collisions at √ s NN = 2.61 GeV, as this is the only experiment in the studied energy range, where the yields of all five strange hadrons are available at the same centrality range.Here, for the most of yields a very good agreement is found, with the exception of K + where the distance is about 3.3 standard  5 Experimental and predicted yields of strange hadrons for selected beam energies and centralities.See text for the choice of these points and discussion.
deviations.We also note that the relative errors of the predictions themselves remain within 4-20% range.

Predictions for the upcoming data and planned experiments
The overall positive results of comparisons of the parametrizations to the experimental yields encourages us to present several predictions for the colliding systems for which the experimental data for strangeness are currently analysed or planned to be obtained.All the A part b values were calculated with the method described in Sect. 2.
Let us first look at the collisions of Ag+Ag at the available energies of 2.41 and 2.55 GeV (current analysis by HADES).Here, the uncertainties of the predicted yields range from about 3% for K + down to 22% for φ.The predictions appear to be in good agreement with the experimental data in the neighbouring region of √ s NN and Focusing on collisions of Au+Au at √ s NN = 2.16 and 2.24 GeV, being the topic of the HADES S520 experiment accepted by the GSI, we notice that these energies probe the edge of available data for K + , but fall below the edges for the other hadron types.Nevertheless one can consider the following ratios of yields: K 0 s /K + , K − /K + , φ/K − and (Λ + Σ 0 )/K 0 s .The first ratio reaches on average 1.03±0.30at lower energy and 0.73±0.17at the higher one.For the closest beam energy where the experimental ratio can be constructed, √ s NN of 2.41, the value is within 0.47-0.55depending on centrality, with uncertainty of about 0.07.Therefore, our predictions are 1.7 (0.7) standard deviations above the compared experimental values for the lower (higher) available energy.Since the agreement for the K + yield at 2. 16 GeV is very good, we conclude that the predictions of K 0 s are reasonable.Moving to K − /K + , we first note that for the lower energy the uncertainties are larger than the yields.For the higher energy, the values of this ratio are about (1.2 ± 0.7) • 10 −3 , where the considerable error is the consequence of the uncertainty of the K − yield.The experimental value closest to the discussed case is 7 • 10 −3 for √ s NN of 2.41 GeV [24] and the K − /K + ratio is known to systematically drop with decreasing the available energy (c.f.Fig. 6.2 in [6]).The ratio predicted at √ s NN = 2.24 GeV is generally in line with the extrapolation of the trend on this plot.For the φ/K − ratio, we notice that at lower energy the uncertainty of the prediction of φ yield reaches 100%.For the higher energy it is not so drastic, but the ratio of yields, found to be 0.11, is burdened with 80% uncertainty.It is not in disagreement with the closest experimental value of this ratio, 0.46 ± 0.12 for √ s NN = 2.42 GeV [24], however, the large uncertainty of our prediction makes it of little use.The predicted values of the (Λ + Σ 0 ) yield are below 10 −10 , which is more than 6 orders of magnitude below the yield of K 0 and K + .As the basic expected strangeness production channel is BB → BK +/0 Y, where B ∈ {N, ∆}, the hyperon and kaon yields should be similar.This is the case e.g. for √ s NN = 2.41 GeV, where the ratio is at the level of 3. Therefore the predictions for the hyperon production at the S520 energies are unreliable, and probably are the effect of the very steep β parameter found by the fitting procedure, already discussed in the previous section.To sum up, our predictions at 2.24 GeV seem reasonable for all the kaons, but at 2.16 GeV -only for K + and K 0 s , while the rest of predictions is unreliable.
Let us discuss the predictions for √ s NN = 3 GeV (STAR's lowest beam energy) for the currently unpublished yields of K + , K 0 s and Λ for three centrality classes applied in [2].For all the centralities we found K − /K + ratio to be (4.8 ± 0.7)%.This value is well in line with the AGS results obtained at midrapidity (c.f.Fig. 3 in [45]).Since the result from our parametric curve for K − follows the experiment nearly exactly, we conclude that the K + yield is predicted also well.The predicted K 0 s /K + ratio rises from 0.32(5) to 0.49 (7) while passing from central to more peripheral collisions.The nearest available energy and centrality at which such a ratio can be constructed from the experimental data is √ s NN = 2.67 GeV and A part b ∼ 72.The ratio is within 0.41 and 0.51, depending on the centrality class of two adjacent data points for K + .Apparently, the predicted values fit nicely to the experimental ones at the adjacent energy, hence the prediction for K 0 s is reasonable too.Let us now focus on the (Λ + Σ 0 )/K 0 s ratio.For the considered centrality classes our prediction points to values between 1.13(16) and 0.94 (10).However, at √ s NN of 2.67 GeV the experimental value of this ratio is found to be between 2.6 and 3.3.This discrepancy is probably the consequence of the previously discussed reaching the saturation region by Eq. 2 for the Λ case.Hence, we conclude that this prediction is not reliable.To sum up, among the unpublished yields of K + , K 0 s and Λ + Σ 0 at √ s NN of 3.0 GeV, the predictions for two former hadrons make sense, but not for the latter one.
The CBM setup, currently being constructed at FAIR, GSI is planned to measure heavy-ion collisions at available energies within 2.7 -5 GeV [46].Here, we should check, how far in terms of available energy our parametrization can be useful.The experimental yields of K ± in this energy region are found to be rising sharply, however they were measured only at midrapidity, so this data cannot be applied directly to our analysis [45].The 4π yields of K ± and Λ are available at √ s NN of 4.86 GeV.However, as seen in Table 5, our parametrization strongly underestimates the experimental yields.The region of available energies on which our fits base, 2.16 -3.0 GeV, is understood to be of (nearly) purely hadronic character.At higher available energies the string degrees of freedom open up and the chiral symmetry restoration during the hadronization process is expected to additionally raise the yields (c.f.Sect.5.2 in [47]).On the side of fitting, the curve defined by Eq. 2 inevitably saturates with energy.Therefore, an extrapolation of our fit should not go far.In Table 5 we present the predictions up to √ s NN = 3.85 GeV only.They were obtained at A part b = 100.To assess credibility of these results, let us first examine the K − /K + yield ratio.Our fits predict that its values should be 0.031(3), 0.055 (12) and 0.057 (18), respectively from the lowest energy upwards.These results can be compared to the experimental excitation function of this ratio, shown in Fig. 3 of [45], although measured at midrapidity.Apparently, two first predictions are in line with the experimental data, but the third one falls below.Regarding the φ/K − ratio, our curves give the values: 0.40 (7), 0.18(4) and 0.15(4) from the lowest energy upwards.They seem to nicely follow the experimental trend of this ratio, as shown in Fig. 20 of [8].Passing to the K 0 s /K + yield ratio, our parametrization proposes, 0.45(3), 0.41(8) and 0.42 (11), respectively from the lowest energy upwards.As mentioned before, the nearest available energy where the experimental data is available is 2.67 GeV, and the value of the ratio is located between 0.41 and 0.51.All the predictions fit very well to these values.Regarding Λ + Σ 0 , the predicted yield of these hyperons does not change with the beam energy.It is a consequence of reaching the saturation by this curve, as discussed before.However, the ( √ s NN , A part b ) coordinates for the lowest CBM energy are situated close to two points where the experimental values of (Λ + Σ 0 )/K 0 s ratio are available.Namely, at (2.67 GeV, 93) and at (2.7 GeV, 336) the ratio is found to be 2.9(4) and 3.2(4), respectively.Our prediction of the ratio at (2.7 GeV, 100) gives 2.47 (18).As the standard deviations between our prediction and these experimental points are, respectively, 0.8 and 1.7, we conclude that at the lowest CBM available energy of 2.7 GeV the predicted ratio (and consequently the yield of Λ + Σ 0 ) is still reliable.To sum up, for the lowest energy all the predictions seem reasonable, for the middle one -all except Λ + Σ 0 and for the highest one (at √ s NN = 3.85 GeV) -only φ/K − and K 0 s /K + ratios make sense.Overall, the predictions the K + and K 0 yields seem to be reliable within the available energy range of [2.16 -3.0] GeV, whereas for K − the lower limit should be lifted to 2.24 GeV.For Λ + Σ 0 and φ the predictions appear to make sense from 2.41 GeV, up to 2.7 (3.0) GeV, respectively.In terms of the ratios of yields, the predictions for the K − /K + seem to be reliable within [2.24 -3.32] GeV, K 0 s /K + makes sense for all the examined range of [2.16 -3.85] GeV, and φ/K − fits the experimental data for [2.41 -3.85] GeV.

Discussion on the α exponent
A consequence of adopting the parametrization of Eq. 2 is the assumption that for every hadron type the α exponent of the yield dependency on A part b does not change with beam energy.We tested this assumption by fitting α for groups of data points at fixed available energies (or a narrow range of them), in the region between 2.3 and 3.0 GeV.We noticed that the fit of the set of yields of K + emitted at √ s NN = 2.63 GeV resulted in an unacceptably high χ 2 /ν of 15.Seemingly, the data point for C+C colliding at this energy could not be aligned with the yields for Ni+Ni.Therefore, we disregarded this dataset.We also did not consider the φ meson data, as for the only series of points at the same available energy of 2.67 GeV, the result α = 1.0 ± 0.6 is burdened with too high uncertainty to contribute noticeably to the analysis [32].Out of 11 accepted series of data points, 8 were characterized by the same available energy within ± 0.01 GeV.For K − , K 0 s and Λ emitted from collisions at √ s NN ∈ [2.6 − 2.7] GeV we noticed some variations of the fitted α values depending on the choice of subsample.Therefore, for each of these cases we grouped together the data for the above-mentioned energy range and fitted the Eq. 2 where the available energy dependence was subject to fitting.We then estimated the additional contributions to the systematic errors by testing the susceptibility of resulting value of α to the removal of any single point from these samples.
Figure 3 shows the results as function of √ s NN and hadron type.Per hadron type we do not observe any statistically meaningful linear change of α with √ s NN .A fit of the linear function to all the points taken together provides the gradient value of 0.11 ± 0.16, so the global conclusion is the same.Therefore, we are encouraged to fit a constant to the whole dataset, which gives α = 1.30 ± 0.02.
We can now compare this result to the α values obtained per hadron, shown in Table 4. Apparently, for K + , K − and φ mesons the agreement is within 1 standard deviation (σ), and for Λ + Σ 0 the distance is 1.8 σ.However, the α exponents for K 0 s appears to be lower by 5 σ units.In [11], a global value of α = 1.45±0.06was extracted for five strange hadron species emitted from Au+Au collisions at the available energy of 2.41 GeV.This finding is 2.5 σ away from our result obtained for all the energies, so rather not in disagreement.We also found that the values of α exponent for the AGS data on charged kaon emission at √ s NN = 4.86 GeV are: 1.18 ± 0.05 for K + and 1.22 ± 0.04 for K − [29].Interestingly, despite the fact that this energy seems to be already situated in the region of string degrees of freedom, the found values are  3 Values of α exponent of Eq. 2 obtained for fixed beam energies and given hadron types.To increase readability at clusters of points, some uncertainties were shifted to the edges of markers.All the errors are symmetric.Dashed (full) lines mark the linear (constant) functions of the best fit to the data.also within about 2σ in agreement with the constant average found in our analysis for √ s NN ∈ [2.24, 3.0] GeV.

Comparison to transport models
It is of interest to compare the predictive power of the found parametrization formulae to that of various transport model calculations.Currently only the collisions of 40 Ar+ nat KCl at the available energy of 2.61 GeV offer the set of yields for all five strange hadrons considered in this work at the same centrality range of 0-35%.We shall treat this data as the basis for comparison.The simulations of heavy-ion collisions were performed within the recent publicly available versions of RQMD.RMF implemented in the Monte-Carlo event generator JAM2.1 [48], SMASH ([49], version 2.1) and UrQMD ( [50], version 3.4).To approximate the nat KCl, the 37 Ar was taken as the target nucleus.The class of 35% most central events was selected by the condition of b < 5.0 fm, resulting from the Glauber model calculations.The SMASH simulations were performed for two variants of the Equation of State (EoS), characterized by the incompressibility modulus κ = 240 and 380 MeV.To stabilize the mean field map each event was processed in 40 ensembles.For RQMD.RMF two momentum-dependent EoS variants were selected: MD2 and MD4.The variants MD1 and MD3 give results in-between.As currently none of these variants is pointed out as the more applicable one to the considered range of energies, the chosen scenarios can be treated as an estimate of the systematic error of the model originating from an uncertainty of EoS.Since we are interested only in the 4π yields, the Coulomb potential was neglected both in RQMD.RMF and SMASH.The publicly available version of UrQMD offers either the cascade approach or the "hard" variant of EoS.We chose the latter one due to better relevance to AA collisions at beam energies in the hadronic sector.The yields of K 0 s were obtained by dividing the found K 0 yield by 2. The multiplicity of φ meson was found by counting first the so-called "reconstructable" K + K − pairs, i.e. these ones, whose invariant mass remains in a peak around m φ = 1.0195GeV/c 2 .The yield was then corrected for BR (φ → K + K − ) = 49.1 % [30].The yields of φ mesons from UrQMD 3.4 were not included, since the proposed dominant mechanism of φ production from decays of heavy resonances [7] was introduced after this version of the model.Since the Σ 0 hyperons were observed in the considered experiments only indirectly (as the feed-down to the Λ spectrum with BR ≈ 100%) in our simulations the yields of these two hyperons were summed up.
In the upper panel of Fig. 4 the yields for Ar+KCl found experimentally are set together with those obtained from the phenomenological parametrization using Eq. 2, and the ones predicted by the transport models.To quantify the overall discrepancy between the experimental yields of five considered hadrons and their predictions obtained within a given model, we used the total standard deviation: where h ∈ {K + , K − , K 0 s , Λ + Σ 0 , φ}.In cases of the transport codes, by the uncertainties of yields we mean the statistical errors.For the worst case of the least abundant φ mesons they are at the level of 10% of the respective yield, so their influence is minor compared to the experimental uncertainty.The values of found z Model are listed in Table 6.The phenomenological predictions of Eq. 2 appear to fare better than RQMD.RMF and UrQMD calculations both overall, and for each hadron type separately.For SMASH, the overall discrepancy is still larger than our prediction.However, this model fares better for K − and in the soft EoS case -for K + .
In two lower panels of Fig. 4 we set together the predictions of yields of strange hadrons emitted from 10% most central collisions of Ag+Ag at the available energies of 2.41 and 2.55 GeV.These predictions cover both the phenomenological parametrization of Eq. 2 and the selected transport models.We encourage the researchers to include our parametrization into comparisons of the experimental yields of hadrons emitted from these collisions.GeV within centrality of 0-10%.Data labelled as "Exp.interpolation" are the predictions of Eq. 2 using the best-fit parameters from Table 4. See text for details on the SMASH, RQMD.RMF and UrQMD simulations and discussion.

Summary
We analysed the available experimental yields of five strange hadron types (K ±,0 , φ and Λ with feeddown from Σ 0 decays) emitted from collisions of heavy ions at √ s NN in the range of 2.16 -3.0 GeV.Since in the original works the estimation of A part b was done within three different approaches, we used the stated centralities and extracted the numbers of participants for all the data with help of one common method, the Glauber Monte Carlo.We then postulated the parametrization of yields of each of analysed hadrons as function of √ s NN and A part b , and obtained the fits with χ 2 /ν between 0.15 and 3.6.We demonstrated the good description of the yields at the boundaries of experimental datasets that were used for fitting.Therefore, we felt encouraged to provide the predictions of yields for the currently analysed HADES data from Ag+Ag collisions at the available energies of 2.41 and 2.55 GeV, accepted Au+Au experiment at 2.16 and 2.24 GeV, some currently unpublished yields from STAR's Au+Au interactions at 3.0 GeV, as well as planned CBM data from Au+Au collisions up to 3.85 GeV.Overall, regarding the yields, the predictions of our phenomenological curves for K + and K 0 s seem to work within [2.16 -3.0] GeV, whereas for K − the lower limit should be lifted to 2.24 GeV.Concerning φ, the range of applicability spans within [2.41 -3.0] GeV, and in case of Λ + Σ 0 the upper limit should be dropped to 2.7 GeV.However, in terms of yield ratios, the applicability of K + /K 0 s and φ/K − can be extended up to 3.85 GeV, whereas for K − /K + the fits currently make sense up to 3.32 GeV.
The accumulated data allows us also to find that the α power exponent of the yield dependence on the A part b does not change with available energy within the current uncertainties.Its average value was found to be 1.30±0.02.This value is in general agreement with the values of α obtained per hadron type via fit of Eq. 2, with the exception of K 0 s .We also used the yields from Ar+KCl at 2.61 GeV to benchmark our parametrization and the predictions of recent publicly available versions of RQMD.RMF, SMASH and UrQMD.For this system, we concluded that our approach currently appears to have the best degree of overall agreement with the experimental data.We therefore encourage the researchers to include it into their prediction and comparison tools.8 Covariance matrices of the best fits of Eq. 2 to experimental yields of studied strange hadrons.All the matrices were found to be symmetric.For K 0 s , Λ + Σ 0 and φ the C parameter was fixed during fitting.

Fig. 2 (
Fig. 2 (Top panels) systematics of yields of K + and φ mesons emitted from heavy-ion collisions in the near-threshold range, as function of available energy in the NN channel, √ s NN , and mean number of participants, Apart b .The points depict the experimental data.The curves are the best fits of Eq. 2 to these data.(Bottom panels:) the maps of standard deviations of the data points from the fitted curves.
Fig.3Values of α exponent of Eq. 2 obtained for fixed beam energies and given hadron types.To increase readability at clusters of points, some uncertainties were shifted to the edges of markers.All the errors are symmetric.Dashed (full) lines mark the linear (constant) functions of the best fit to the data.

Fig. 4
Fig. 4 Comparison of experimental, interpolated and simulated production yields of five strange hadrons from Ar+KCl at √ s NN = 2.61 GeV within centrality of 0-35% and Ag+Ag at 2.41 and 2.55

Table 2
Yields of K − and φ mesons emitted from nuclei colliding at given available energies (as well as beam kinetic energies) and centrality ranges.Values of Apart b are shown as published originally and according to our Glauber Monte Carlo analysis (see text for details.) √ s NN T b /A Centrality Ratio of Apart b values for K + data obtained in this work to those published in the original papers.Numbering scheme of data points is as in Table1.Full dots correspond to data where the geometrical model was originally used, squares refer to the Optical Glauber, and triangles -to the Glauber MC.

Table 4
Results of fit of curve 2 to the yields of strange hadrons as function of Apart b and √ s NN .For K ± all the parameters were fitted.For the other hadrons C was fixed during the fit.The data points are listed in Tables 1, 2 and 3.For prediction of yields see the full covariance matrices in Tab. 8 and discussion in text.

Table 6
Standard deviations (c.f.Eq. 3) between the experimental yields and the ones predicted by different approaches, for hadrons emitted from Ar+KCl at √ s NN = 2.61 GeV at centrality range of 0-35%."Exp.interpolation" are the predictions of Eq. 2.Last column is the sum of standard deviations over five hadrons.