Dynamical Freeze-Out Phenomena: The Case of K ± , φ Transverse Momentum Spectra in Collisions of Au(1.23 A GeV) + Au

: We argue for a continuous (dynamical) kinetic freeze-out of K ± , φ observed at midrapidity in collisions of Au(1.23 A GeV) + Au. The simulations, by means of a transport model of Boltzmann-Ühling-Uhlenbeck (BUU) type, point to time independent transverse momentum slope parameters after 20 fm/c. The complex interplay of expansion dynamics and strangeness production/exchange/absorption as well as elastic scatterings involved in the reaction network does not support the previous interpretation of a late freeze-out of K − due to larger cross sections.


Introduction
The KaoS (Kaon Spectrometer) Collaboration [1] parameterized the transverse momentum spectra of K ± , measured in collisions of Au(1.5 A GeV) + Au at midrapidity, by means of Jüttner type phase space distributions. Ignoring the flow, the distributions become Boltzmann-like, depending essentially on a slope parameter which is often called the "kinetic freeze-out temperature", T. The experimental fact of T K + > T K − has been interpreted in [1] as evidence for an expanding and cooling fireball with antikaons (K − ) decoupling at a later, thus cooler, stage due to their significantly larger cross section. The kaons (K + ) decouple at an earlier, thus hotter, stage.
An alternative interpretation has been put forward by the HADES (High Acceptance Di-Electron Spectrometer) Collaboration [2] for the reaction Au(1.23 A GeV) + Au: K + and K − can decouple under the same circumstances, i.e., at the same fireball temperature, thus T K + = T K − at kinetic decoupling where the elastic interactions cease. However, the unexpectedly large yield of φ changes the late spectra due to the decays φ → K + K − . Since the multiplicities hold N K + N K − , the impact of φ decays on T K + is minor. As observed in simulations [3], T K − f rom φ decays is approximately 2 3 T φ for distributions with T K ± ≈ T φ prior to φ decays. As a result of the 26 ± 8% contribution of K − from φ decays, the observed final slope parameters in fact obey T K + > T K − . In order to cope with the experimental results, a sufficiently large yield of φ is a prerequisite for such an interpretation [2]. (For an early assessment of relating K − and φ yields, see [4], and for the first data-based quantitative analysis of φ → K + K − feed-down, cf. [5]. For a comprehensive survey on strangeness production in near-threshold heavy-ion and proton-nucleus collisions, cf. [6], and for new simulation tools, see [7][8][9][10].) Inspired by the controversy of the KaoS and HADES interpretations, we performed kinetic theory simulations of the collisions Au(1.23 A GeV) + Au and studied the time evolution of the slope parameters T K ± ,φ . Our findings can be summarized as follows: T K ± ,φ (t > 20 fm/c) ≈ const. (It should be stressed that the slope parameters T K ± ,φ refer to transverse momentum spectra at midrapidity of all respective hadron species K ± and φ; they are not related to a local or global fireball temperatures steered overwhelmingly by nucleons and their excitations). We interpret such a behavior as dynamic freeze-out. The energy dependence of the various inelastic and elastic cross sections combined with the proper dilution upon expansion of the fireball give rise to the conspiracy of T K ± ,φ (t > 20 fm/c) ≈ const. There is neither the need nor the possibility to define in an unambiguous manner the freeze-out temperature at a certain instant of time within such a kinetic theory approach.
The dynamical freeze-out phenomenon has been familiar for some time for its role in primordial nucleosynthesis in the early universe, where it refers to the abundances of light isotopes. As a result of the energy dependence of cross sections, which translate into a temperature dependence of reaction rates, combined with the temporal temperature dependence due to cosmic expansion, the isotopic abundances Y i stay constant after about a thousand seconds world age: Y i (t > 10 min) ≈ const. While the starting values are determined by statistical nuclear equilibrium [11], the final values are non-Markovian, i.e., they depend on intermediate stages, and are very specific for cross sections (and some other parameters of the system). In particular, the multitude of final (late) abundances cannot be related to chemical equilibrium values at a certain common temperature.
A prototypical freeze-out model is provided in paragraph 5.2 in [11], see Figure 5.1 therein. On the basis of the momentum-integrated Boltzmann equation, the normalized abundance of a massive-particle species is seen to follow for some time the fiducial equilibrium abundance but levels off at certain point and stays constant afterwards, never reaching the ever-dropping equilibrium value. The leveling off is determined by a combination of the expansion (cooling) rate and the reaction rate based on a thermally averaged cross section. That is, at and after freeze-out, the expansion rate exceeds the reaction rate.
From such a perspective, it is astonishing that the various hadron abundances in heavy-ion collisions can be described by chemical equilibrium values at a common (albeit beam-energy dependent) temperature modulo a normalization volume. At LHC (Large Hadron Collider) energies, no other parameters are needed to uncover the hadron and isotopic (including antinuclei) yields over nine orders of magnitudes [12]; at lower beam energies, the baryo-chemical potential becomes important, which is also beam-energy dependent, to describe a multitude of hadron yields, cf. [13,14] for examples.
Besides the chemical freeze-out, related to abundances, the kinetic freeze-out, related to the momentum distributions of hadrons, is also of interest as a signature of the "hadronic life", e.g., after the hadronization of deconfined strong-interaction matter at LHC energies. This freeze-out dynamics may be flavor dependent and may be different for ground state hadrons and short-living resonances [15]. In the fragmentation region, large net-baryon densities are expected in ultrarelativistic heavy-ion collisions [16]-quite similar to conditions achieved in relativistic heavy-ions. Using rare (i.e., strange and charm) probes of strongly compressed baryon matter is a central part of the research program of the CBM collaboration [17,18]. These facets of ultrarelativistic heavy-ion collisions in turn are linked to medium-energy (relativistic) heavy-ion collisions, where one gains complementary important information on hadronic many-body dynamics.
After this digression on freeze-out phenomena, let us return to the primary goal of our study-the time evolution of K ± , φ slope parameters in a kinetic theory simulation of Boltzmann-Ühling-Uhlenbeck (BUU) type for the reaction Au(1.23 A GeV) + Au. The BUU model is briefly described in Section 2, and its numerical results are presented in Section 3. Since we employ a code version which has been successfully utilized in [19] for the reaction Ar(1.75 A GeV) + KCl, we present analog results of the time evolution of T K ± ,φ in Appendix A for comparison. Section 4 summarizes the paper.

BUU Code
Here, we employ the same BUU code as utilized in [19] for collisions Ar(1.75 A GeV) + KCl. Only the beam energy, system size, proton-to-neutron ratio, and impact parameter range are adapted; all other parameters are frozen, see Table 1. (In [6], you can find an extended discussion on scalar and vector contributions to the effective in-medium masses (cf. Equation (60) for the implementation in IQMD (Isospin-dependent Quantum Molecular Dynamics), or Equations (20) and (21) for the derivation from a chiral Lagrangian) and the related issue of spectral functions (cf., e.g., [20,21]). In this respect, the often used parameterization m(ρ) ∝ ρ∆m is rather schematic). The code was compared in [22] with other transport codes and some peculiarities have been identified. Nevertheless, as shown in [19], the BUU code copes successfully with the data from [23]. The decisive difference of the reactions Au(1.23 A GeV) + Au and Ar(1.75 A GeV) + KCl are the significantly lower beam energy and the significantly larger system size of the former. This leads us to expect a larger sensitivity to many-body effects, in particular in-medium effects, since the reactions with strange mesons are deeper below the respective thresholds. Figure 1, left panel, exhibits a survey on the thresholds and experiments performed up to now in the threshold region.

Impact Parameter Dependence
It happened that K ± , φ transverse momenta spectra and rapidity distributions of the Ar(1.75 A GeV) + KCl data [23] could be described very well by the impact parameter b = 3.9 fm, see Appendix A. Following such a strategy for Au(1.23 A GeV) + Au, we see that an optimum description of the data [2] for the centrality class 0-40% is accomplished by b = 9 fm for K + and b = 10 fm for K − , φ, see Figures 2-4. As a compromise, we henceforth use b = 9 fm to avoid subtleties of certain impact parameter averaging procedures according to weighting ∝ b db. In fact, Table II in [26] attributes the 0-40% centrality class to the impact parameter interval b = 0-9.3 fm with a mean of 6.2 fm. The upper panel in Figure 8 of [26] indicates that an impact parameter range 7-11 fm centered at 9 fm corresponds to the centrality class 30-40%. (Table II in [26] quotes a mean impact parameter of 8.71 fm for that centrality class). Therefore, the selection of one "representative impact parameter" must be considered with caution. To get some feeling on the impact parameter dependence, we exhibit in the following, the sequence of b = 1-10 fm in steps of 1 fm with the color code displayed in the right panel of Figure 1.

Figure 2. Left panel:
Rapidity spectra of K + mesons in the center-of-mass system for impact parameters 1 ≤ b ≤ 10 fm (top to bottom). Right panel: Transverse mass spectra within the rapidity interval −0.1 ≤ y ≤ 0.1 for the same impact parameter range as in the left panel. The black symbols represent the experimental data from [2] (centrality 20-40%). The color code for the impact parameters is depicted in the right panel of Figure 1.   The transverse momentum spectra unravel some deficits of our model: T K ± are somewhat too low (note that the first data point is at the lower edge of the histogram, while the last data point is at the upper edge of the histogram in Figures 2 and 3), and T φ is much too low as evidenced by the much flatter m t − m 0 distribution in Figure 4). These facts are quantified by the slope parameters, see Figure 5 below. For φ, these deficits cannot be cured by some impact parameter averaging since we fail to meet the experimental values of T φ for all values of b. On the other hand, [2] quotes slope parameters of about 91 MeV (K + , centrality class 30-40%) and 69 ± 7 MeV (K − , centrality class 20-40%), which are not too small in comparison with our results, as quantified in Figure 5 below. As in the experiment in [2], we define the slope parameters by at midrapidity y 0 in the rapidity interval y = y 0 ± 0.1. In addition, a normalization factor C is attributed separately to each species. The transverse mass is defined by m 2 t = p 2 t + m 2 0 with rest masses m 0 of the considered species. Henceforth, we denote T B (y 0 ) = T which is again specific for each species. The experimental data [2] for the centrality class 0-40% are represented as gray bars. "Total K − " means inclusion of K − from φ decays, while "direct K − " is without that contribution. Bottom figure: Corresponding reduced χ 2 values to the fit parameters T and C.

Nucleon Density Evolution
We refrain here completely from determining any local temperature. Instead, to visualize the time evolution of the system by the nucleon density in the reaction plane, we exhibit in Figure 6 several snapshots for the representative impact parameter b = 9 fm. The maximum nucleon density is about 2ρ 0 (or 2.75ρ 0 for b = 1 fm) at t ≈ 12 fm/c in the central cell of volume 1 fm 3 . The high-density stage with ρ ≥ ρ 0 is for t = 7 · · · 20 fm/c, see Figure 7.

Time Evolution of K ± , φ Rates
Production rates of K ± , φ as well as Λ + Σ ±,0 as a function of time (left panel) are displayed in Figure 8. The analog elastic collision rates of K ± with nucleons and the absorption rates are exhibited in Figure 9. Concerning the production rates, one sees some delay of K − , φ relative to K + , Λ + Σ ±,0 .
For brevity, Σ * denotes all Σ ±,0 . Focusing still on the evolution, one cannot see a later decoupling of K − ; instead, the last elastic interactions of K + seem to go on for a somewhat longer time, see Figure 10. These investigations are aimed at elucidating whether there is a clear time ordering of production and freeze-out of K ± and φ.     Since the in-medium masses of K ± are strikingly different (see Table 1), it has been argued that production/rescattering/absorption/last elastic interaction probe different densities. In fact, K − production happens at somewhat larger densities (see right panel of Figure 8); nevertheless, the K ± production rates peak at about 1.6ρ 0 . Elastic rescatterings peak at 1.3ρ 0 (K + ) and 1.4ρ 0 (K − ), respectively, while the absorptions acquire maxima at 1.3ρ 0 (K − ) and 1.7ρ 0 (K + ), see right panel in Figure 9. The maximum of last interaction rates is at 1.2ρ 0 (K ± ), though with more last interactions of K + at lower densities, ρ < ρ 0 , see right panel in Figure 10. This seems to be in contrast to the above-mentioned time and density ordering of K + and K − collisions.

Time Evolution of T K ± ,φ
The time evolution of the parameters T K ± ,φ and C K ± ,φ defined in Equation (1) is exhibited in Figure 5 together with the K ± , φ multiplicities. The most striking point is T K ± ,φ ≈ const for t > 20 fm/c and all impact parameters; T K − even increases slightly after 20 fm/c for the most central collisions. This invalidates the expectation that the slope parameter T can be related to a local medium temperature. Rather, the interaction rates (see Figures 9 and 10) drop to zero at t > 40 fm/c, meaning the dynamic freeze-out in the same spirit as in big-bang nucleosynthesis. Prior to 40 fm/c, but after 20 fm/c, the intricate network of production, absorption, and elastic reactions on top of the overall expansion (see Figures 6 and 7) result in T K ± ,φ ≈ const.
A completely different picture emerges when considering the parameters T K ± ,φ in the central cell: after 15 fm/c (K ± ) or 17 fm/c (φ), the values of T K ± ,φ rapidly drop and become < 20 MeV for t > 30 fm/c, see Figure 11. Of course, such a correlation of momentum space and position space is experimentally hardly accessible. We emphasize in this context the importance of cross sections. For instance, scaling up the total K − + anything cross section by a factor of ten increases the maximum of T K − by 15 MeV; switching off the absorption channels lets the maximum of T K − further increase towards 100 MeV, see Figure 12. Switching off the absorption for our standard setting of cross sections (see Appendix A in [27]) also increases the maximum of T K − by about 15 MeV.  A quantity, which can serve as a proxy of the temperature in a local off-equilibrium situation, is the mean kinetic energy. We define E kin = E − m * in the center-of-mass system. The effective in-medium mass is defined by m * = m + ∆mρ with values of ∆mρ 0 = ∆m(ρ 0 ) listed in Table 1. Interestingly, we observe E kin,K − > E kin,K + at t < 18 fm/c, see Figure 13-left panel. The pronounced dropping of E kin,K − at t > 13 fm/c can be attributed to the change in the effective in-medium mass m * with dropping nucleon density. In fact, by switching off the effective in-medium masses, i.e., m * → m 0 , the dropping of E kin,K − ceases, see the dashed green curve in right panel of Figure 13. Since the in-medium modifications of K + and φ are minor, there is no noticeable impact on E kin,K + ,φ vs. t. The solid green curves in Figure 13 include the K − from φ decays. The "cooling" of the K − spectrum is due to the above-mentioned fact that the mean kinetic energies of decay-K − are less than the mean kinetic energies of K − in a medium with a temperature scale O(100) MeV.

Discussion and Summary
Contrary to the interpretations [1,2] of T K + > T K − in relativistic heavy-ion collisions at subthreshold beam energy, the use of a BUU transport simulation points to a more complex picture. Figure 5 suggests that from the very beginning of the fireball expansion, the relation T K + > T K − is established and does not change in the subsequent evolution. That is, the complicated interplay of various reaction channels involved in strangeness production, exchange, and absorption up to elastic scatterings cause flatter K + transverse momentum spectra at midrapidity. We emphasize that the slope parameters T K ± ,φ are used to parameterize these transverse momentum spectra. Enforcing a correlation to position space, e.g., by considering only mesons in the central cell, substantially changes the behavior of T K ± ,φ . Instead of staying approximately constant (as for all mesons in the fireball), it rapidly drops after achieving a maximum. The picture becomes more obscured by analyzing the mean kinetic energies of K ± , φ. The coupling of the effective in-medium mass to the nucleon density causes a significant reduction in the K − kinetic energy below the K + kinetic energy, while without in-medium effects, the K − kinetic energy is below the K + kinetic energy, both being roughly constant in time.
Irrespectively of the temporal aspects, we also tried to check the hypothesis that at the production instant, the K − have lower kinetic energies due to peculiarities of the production channels. However, subsequent absorption and scattering/transfer reactions make such a hypothesis unconvincing.
Our interpretation is hampered in some details since we fail to accurately reproduce the φ multiplicity and the K ± , φ slope parameters in the analyzed reaction Au(1.23 A GeV) + Au. The same code, however, describes well the available data of the reaction Ar(1.756 A GeV) + KCl (see Appendix A). It happens that the weights of various channels are fairly different, see the panels in Figure 14 (the ρN → φN channel is special according to [28] (Figure 2): it increases towards small energies in the entrance channel. As a consequence, the relative weight of this channel is larger at smaller beam energies. Analogous, at small energies in the entrance channel, the cross sections obey σ ∆∆→φX > σ N∆→φX > σ NN→φX , according to [28], while at larger entrance energies σ NN→φX dominates. These trends manifest themselves in the smaller contribution of the NN channel at the lower beam energy. Furthermore, the decays of heavy baryon resonances with a φ in the exit channel (cf. [10,29]) are not included).
This means that one needs to investigate different system sizes and different beam energies for benchmarking both the many microscopic input data and their processing in simulation codes by a limited number of observables. The herein presented results can serve as a useful reference of analyses of strangeness dynamics of the planned experiments of CBM (Compressed Baryonic Matter).  In summary, we put forward arguments in favor of a dynamical (i.e., continuous) kinetic freeze-out of strangeness carrying mesons in subthreshold heavy-ion collisions. The in-depth exploration of strangeness dynamics in near-threshold heavy-ion collisions at SIS18 (SchwerIonenSynchrotron) is a valuable prerequisite to the investigation of multi-strange hadron phenomena at higher beam energies, e.g., at SIS100, NICA (Nuclotron-based Ion Collider fAcility), JPARC (Japan Proton Accelerator Research Complex), etc. The insights gained in the strangeness sector pave the way for future analog studies of charm degrees of freedom as probes of compressed baryon matter.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Ar(1.756 A GeV) + KCl
To enable an easy one-to-one comparison of the results presented in section III, we recollect some analog results for collisions Ar(1.756 A GeV) + KCl in Figures A1-A11. Further details can be found in [19].
Instead of scanning through impact parameters, as done in Figures 2-4, we compare in Figures A1 and A2 the rapidity and transverse momentum spectra with data for one optimized impact parameter, b = 3.9 fm. The latter value corresponds to the mean impact parameter enforced by the trigger setting LVL1 in the experiment [23]. Note that the LVL1 trigger condition enforces a mean impact parameter b 4 fm, cf. [19]. In general, the db b weighting attributes a higher representativity to larger impact parameters. Moreover, one could argue on a b /R 2 scaling which can be disturbed by the centrality class selection in Au(1.23 A GeV) + Au vs. LVL1 selection in Ar(1.756 A GeV) + KCl (here we use the estimate R(A) ≈ 1.25 fm A 1/3 for A = 197 and A = 40 yielding the radii 7 fm and 4.1 fm). We emphasize the quite accurate agreement of data and simulations for K ± . The φ slope parameter is acceptable (see the right panel in Figure A3), but the yield, to be extracted from the rapidity distribution (the left panel in Figure A3), is only marginally consistent with the data.  Figure A1. Left panel: Rapidity spectrum of K + mesons in the center-of-mass system for impact parameter b = 3.9 fm. Right panel: Transverse mass spectra for b = 3.9 fm for the six rapidity intervals 0.1 ≤ y lab ≤ 0.2 to 0.6 ≤ y lab ≤ 0.7 (bottom to top) with respective scaling factors 10 0 to 10 5 . The black symbols represent the corresponding experimental data from [23] under the LVL1 trigger setting.  Figure A2. The same as Figure A1, but for K − mesons.  Figure A3. Left panel: The same as Figure A1, left panel, but for φ mesons. Right panel: Transverse mass spectrum for b = 3.9 fm for the rapidity interval 0.2 ≤ y lab ≤ 0.6 and corresponding experimental data from [23] as black symbols.
The high-density stage is much shorter for Ar(1.756 A GeV) + KCl (see Figures A4 and A5) than for Au(1.23 A GeV) + Au (see Figures 6 and 7). Correspondingly, the various rates of production/absorption/elastic scattering/last elastic scattering of K ± , φ are concentrated in shorter time intervals, see the left panels of Figures A6-A8 and compare them with Figures 8-10. The same rates, however, as a function of local density, look very similar for the production channels (compare the right panels of Figure A6 with Figure 8), while the elastic and absorption rates of K ± peak at somewhat higher local densities, with the exception of K + absorption (compare the right panels in Figure A7 with Figure 9). The last elastic K ± , φ scatterings are fairly smoothly distributed over all densities (see the right panel in Figure A8), while for Au(1.23 A GeV) + Au, an apparent peaking at ρ = (1 · · · 1.5)ρ 0 occurs (see the right panel in Figure 10). Such considerations have the motivation to elucidate whether the strangeness-carrying probes are specific for certain density ranges. x (fm) ρ/ρ 0 z (fm) Figure A4. The same as Figure 6, but for Ar(1.756 A GeV) + KCl and for the impact parameter b = 3.9 fm.     Figure A7. The same as Figure 9, but for Ar(1.756 A GeV) + KCl and for b = 3.9 fm.  Figure A8. The same as Figure 10, but for Ar(1.756 A GeV) + KCl and for b = 3.9 fm. Note our restriction to the last elastic interaction. In contrast, Figure 10 in [19] is based on a counting scheme of all last interactions. As a result of the partially perturbative treatment [19] of K − absorption in dilute nuclear matter, the apparently late K − interaction is overestimated.
The patterns of the time dependence of slope parameters and normalizations of Ar(1.756 A GeV)+ KCl look very similar to the case of Au(1.23 A GeV) + Au, where the experimental values of the former are nicely reproduced, see Figure A9 and compare it with Figure 5. An analog statement holds for the slope parameters and normalizations in the central cell-again within shorter time intervals for Ar(1.756 A GeV) + KCl (see Figure A9 and compare it with Figure 11).  Figure A9. The same as Figure 5, but for Ar(1.756 A GeV) + KCl and for impact parameters 1 ≤ b ≤ 6 fm. Top figure: Slope parameter T parameterizing the transverse momentum spectra (top row), normalization C (middle row), and multiplicity N (bottom row) of K + (left column), K − (middle column) and φ mesons (right column) as functions of time for impact parameters 1 ≤ b ≤ 6 fm. The experimental data [23] are represented as gray bars. "Total K − " means inclusion of K − from φ decays, while "direct K − " is without that contribution. Bottom figure: Corresponding reduced χ 2 values to the fit parameters T and C.
Comparing the time evolution of the mean kinetic energies of K ± , φ with and without KN potentials, one observes analog patterns in both collision systems, see Figures 13 and A11. However, the local maximum for K − and the onset of the flat sections for K + and φ are achieved earlier for Ar(1.756 A GeV) + KCl, and the different beam energies reflect themselves in higher values. The sharp drops in the K − mean kinetic momentum (green curves) at t = 60 fm/c are due to the K − component from φ decay added at the end of the simulation run time. This is the "cooling of the K − spectra" by φ decay [2].  Figure A11. The same as Figure 13, but for Ar(1.756 A GeV) + KCl and for b = 3.9 fm.