Probing Strangeness Canonical Ensemble with $K^{-}$, $\phi(1020)$ and $\Xi^{-}$ Production in Au+Au Collisions at ${\sqrt{s_{\rm NN}}=\rm{3\,GeV}}$

Relativistic heavy ion physics is aiming at the detailed investigation of phase structures of strongly interacting matter, governed by quantum chromodynamics (QCD), under extreme conditions of high temperature and density Adams et al. (2005); Akiba et al. (2015); Busza et al. (2018). Particle production has been studied to investigate properties of the produced QCD matter in heavy-ion collisions. The strange quark mass is comparable to the QCD scale ($\Lambda_{\rm{QCD}}$ $\sim\textup{200 MeV}$), therefore strange quark plays an important role in studying the QCD phase diagram and the Equation-of-State (EoS), particularly in the high density region Rafelski and Muller (1982); Koch et al. (1986); Aichelin and Ko (1985); Fuchs (2006); Ko (2018); Adamczewski-Musch et al. (2019a).

Statistical thermal models have often been used to characterize thermal properties of the produced media Rafelski and Danos (1980); Cleymans and Satz (1993); Braun-Munzinger et al. (1995); Becattini et al. (1998); Braun-Munzinger et al. (1999); Florkowski et al. (2002); Braun-Munzinger et al. (2001, 2004); Cleymans et al. (2004); Petráň and Rafelski (2010); Andronic et al. (2018). In these models, grand canonical ensemble (GCE) and canonical ensemble (CE) statistical descriptions can be applied to conserve electric charge, baryon number, and strangeness number in order to compute the final state particle yields. Both GCE and CE models are able to describe various particle yields including strange particles produced in heavy-ion collisions at RHIC and the LHC at center-of-mass energy ($\sqrt{s_{\rm NN}}$) greater than 7.7 GeV. It has been argued that at lower energies, strangeness number needs to be conserved locally on an event-by-event basis described by the CE, which leads to a reduction in the yields of hadrons with non-zero strangeness number (“Canonical Suppression”) Rafelski and Danos (1980); Redlich and Tounsi (2002); Rafelski and Letessier (2002), but not for the $\phi(1020)$ meson with zero net strangeness number (S=0). The $\phi/K^{-}$ ratio is expected to increase with decreasing collision energy in models using the CE treatment for strangeness, opposite to the trend in the GCE treatment. The canonical suppression power for $\Xi^{-}$ (S=2) is even larger than for $K^{-}$ (S=1). The $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios offer a unique test to scrutinize thermodynamic properties of strange quarks in the hot and dense QCD environment.

In heavy-ion collisions, the near/sub-threshold production of multi-strange hadrons can be achieved from the multiple collisions of nucleons, produced particles, and short-lived resonances Zeeb et al. (2004). The particle production in heavy-ion collisions below its free nucleon-nucleon (NN) threshold ($\sqrt{s_{\rm NN}}$ $\sim$2.89 GeV for $\phi$ and $\sim$3.25 GeV for $\Xi^{-}$) is expected to be sensitive to the stiffness of the nuclear EoS at high density Yong et al. (2021), as it is for single-strange hadrons Aichelin and Ko (1985); Fuchs (2006). The near/sub-threshold production further provides the possibility to observe exotic states of QCD matter McLerran and Pisarski (2007) and signatures of “soft deconfinement” Fukushima et al. (2020).

Previous measurements show that the $\phi/K^{-}$ ratio in heavy-ion collisions stays remarkably flat ($\sim$0.15) at collision energies ${\sqrt{s_{\rm NN}}>\textup{5 GeV}}$ Back et al. (2004); Alt et al. (2008a); Adam et al. (2020). At collision energies close to or below the $\phi$ and $\Xi$ NN-thresholds, recent measurements of $\phi/K$ and $\phi/\Xi$ ratios from HADES and FOPI have achieved a significance about 2.2-3.8 sigma in central heavy ion collisions, and the results indicate a relative enhancement compared to those at high energies  Agakishiev et al. (2009a); Piasecki et al. (2015); Gasik et al. (2016); Adamczewski-Musch et al. (2018), indicative of the applicability of the CE description for strangeness production at these energies. Measurements from $\pi$ or proton induced nuclear reactions Muto et al. (2007); Adamczewski-Musch et al. (2019b) suggest that absorption in cold nuclear matter may play a role in the $K^{-}$ and $\phi$ production yields in nuclear collision at low energies. In this Letter, we report high precision measurement of $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios in Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{3\,GeV}}$ from the STAR experiment.

The dataset used in this analysis was collected under the fixed target (FXT) setup Meehan (2016) in the 2018 RHIC run. A single beam was provided by RHIC with total energy equal to 3.85 GeV/nucleon and incident on a gold target of thickness 0.25 mm, corresponding to a 1% interaction probability. The target is installed inside the vacuum pipe, 2 cm below the center of the beam axis, and located 200 cm to the west of the center of the STAR detector. The main detectors used are the Time Projection Chamber (TPC) Anderson et al. (2003); Meehan (2016), the Time of Flight (TOF) detector Llope (2012); Meehan (2016), and the Beam-Beam Counter (BBC) Whitten (2008). The trigger is provided by the signal in the east BBC detector and at least five hits in the TOF detector. To best utilize the detector band-width, the beam-on-target collision rate was tuned to around 1.5 kHz, and the pileup contribution to the triggered event is $<1\%$ Abdallah et al. (2021a). Tracking and particle identification (PID) are done using the TPC and TOF. Both the TPC and TOF detectors have full azimuthal coverage within a pseudorapidity range of 0$<$ $\eta$ $<$ 1.88 for the TPC and 0$<$ $\eta$ $<$ 1.5 for the TOF in FXT mode Anderson et al. (2003); Llope (2012); Meehan (2016). Events are selected with the offline reconstructed collision vertex within 1.5 cm of the target center along the beam direction. Approximately 2.6$\times 10^{8}$ minimum bias (MB) triggered events passed the selection criteria and are used in this analysis.

The centrality class is selected using measured charged particle total multiplicity within the TPC acceptance. A Monte Carlo Glauber model, used in conjunction with a negative binomial distribution to model particle production in hadronic collisions, is optimized in order to best match the data and determine the centrality class Ray and Daugherity (2008); Abdallah et al. (2021a). Due to the trigger inefficiency in the low multiplicity region (corresponding to the most peripheral collisions), we only report the results from the 0–60% centrality class in this paper. In addition, in order to reduce the pile-up contamination, events above the reference multiplicity of 195 are removed from the most central centrality class.

$\phi$ mesons are reconstructed via the decay channel $\phi\rightarrow K^{+}K^{-}$ with a branching ratio (BR) of ($49.2\pm 0.5$)% , while the $\Xi^{-}$ hyperons decay via $\Xi^{-}\rightarrow\Lambda\pi^{-}\rightarrow p\pi^{-}\pi^{-}$ with a BR of ($63.8\pm 0.5$)% Zyla et al. (2020). $\Xi^{-}$ reconstruction is performed using the KFParticle package based on the Kalman Filter method Kisel (2020); Adam et al. (2021). The charged tracks are reconstructed with the TPC in a 0.5 T uniform magnetic field, and are required to consist of at least 20 TPC hits (out of a maximum of 45) and have a ratio between the number of hit points and the maximum possible number of hit points larger than 0.52 to ensure good tracking and avoid track splitting. The TPC tracking performance with these requirements in the FXT data is similar to that in other data taken in the collider mode. Monte Carlo simulations also reproduces the distributions of various tracking variables. The charged tracks are identified via a combination of the ionization energy loss ($dE/dx$) measurement with the TPC and the time-of-flight ($tof$) measurement with the TOF Shao et al. (2006); Xu et al. (2010). The resolution-normalized $dE/dx$ or $\beta$ deviation from the expected values are used for the PID selection. A minimum $p_{T}$ cut of 0.2 GeV/$c$ is required in the analysis. Since the $K^{-}/\pi^{-}$ ratio is much smaller than the $K^{+}/\pi^{+}$ ratio at low energies, to reduce the contamination from $\pi^{-}$ and $e^{-}$ tracks, a strict PID criterion for $K^{-}$ is implemented by requiring the TPC $dE/dx$ and TOF $\beta$ to be within three standard deviations of the expected values. $K^{+}$ tracks used for the $\phi$ analysis are selected with a hybrid algorithm, in which the TPC $dE/dx$ requirement is applied at low momentum $p<0.5$ GeV/$c$ while an additional TOF $\beta$ requirement is imposed at $p>0.5$ GeV/$c$. In the $\Xi$ analysis, proton and $\pi^{-}$ tracks are identified by requiring the TPC $dE/dx$ to be within three standard deviations of the expect values and the TOF $\beta$ requirement is only applied when there is a valid measurement.

Figure 1 (a) shows the invariant mass distribution of $K^{+}K^{-}$ pairs in the transverse momentum ($p_{T}$) region of 0.4–1.6 GeV/$c$ for 0–60% central collisions. The combinatorial background is estimated with the mixed-event (ME) technique in which $K^{+}$ and $K^{-}$ from different events of similar characteristics (centrality, event plane angle) are paired. The mixed-event spectra are normalized to the same-event (SE) distributions in the mass range of 1.04–1.08 GeV/$c^{2}$. After the subtraction of the combinatorial background, the remainder distribution, shown as red solid circles, is fitted with a Breit-Wigner function for the signal plus a linear function which represents the remaining correlated background ($<1\%$) from a partial reconstruction of strange hadrons. The $\phi$ meson raw yields are extracted from the Breit-Wigner function fit within the corresponding $\pm$3$\Gamma$ mass window (mean value $\mu\sim$ 1.0194 GeV/$c^{2}$, full-width-at-half-maximum $\Gamma\sim$ 6.5 MeV/$c^{2}$). The extracted $\phi$ signal shape is consistent with its intrinsic properties convoluted with the detector smearing effect due to finite momentum resolution ($<3\%$ for single track). Note that a Voigt function has been used to extract the signal counts as a cross check, and the extracted yields are consistent with the default value within uncertainties. Figure 1 (b) shows the invariant mass distribution of $\Lambda(p\pi^{-})\pi^{-}$ in the $p_{T}$ region of 0.5–2.0 GeV/$c$ for 0–40% central collisions. The combinatorial background is estimated with the rotating daughter (Rot) method, in which a daughter track of $\Xi^{-}$ is rotated by a random angle between 150 to 210 degrees in the transverse plane. The rotated spectra are normalized to the same-event distributions in the mass ranges of 1.30–1.31 and 1.34–1.35 GeV/$c^{2}$. After the combinatorial background is subtracted, the $\Lambda\pi^{-}$ invariant mass distribution is fitted with a Gaussian for the signal plus a linear function for the remaining correlated background ($<1\%$). The $\Xi^{-}$ raw yields are obtained via histogram bin counting from the invariant mass distributions with all background subtracted within mass windows of $\pm$3$\sigma$ ($\mu\sim$ 1.3222 GeV/$c^{2}$, Gaussian width $\sigma\sim$ 2.0 MeV/$c^{2}$). The signal-to-background (S/B) ratio for $\phi$ and $\Xi^{-}$ within the reconstructed mass windows is about 0.7 and 0.6 respectively. The reconstructed $\phi$ and $\Xi^{-}$ acceptances ($p_{T}$ vs. $y_{cm}$) in the collision center-of-mass frame are shown in Fig. 1 (c) and (d), respectively. The target is located at $y_{cm}=-1.05$, using the convention where the beam travels in the positive direction. The red curve represents the TPC and TOF acceptance edge.

Particle raw yields are calculated in each centrality and $p_{T}$ bin within each rapidity slice. The raw yields are corrected for the TPC acceptance and tracking efficiency, the particle identification efficiency, and the TOF matching and PID efficiency. The TPC acceptance and tracking efficiency is obtained using the standard STAR embedding technique Adamczyk et al. (2017); Adam et al. (2020), in which a small number of MC tracks are processed through the GEANT (v3.21) simulation Brun et al. (1994), then mixed with the real data and reconstructed using the same algorithm as in the real data. The TPC PID, TOF matching and PID efficiencies are obtained from the data-driven method similar as in Ref. Adam et al. (2019). The final average reconstruction (including acceptance etc.) efficiency is $\sim$0.30, 0.04, and 0.02 for $K^{-}$, $\phi$ and $\Xi^{-}$, respectively. MC embedding simulation also reproduces various topological variables used in the $\Xi^{-}$ reconstruction. As a cross-check, we conducted the measurement of $\Xi^{-}$ lifetime from the same data and the result is $164.2\pm 6.6$ (stat.) ps, consistent with the PDG value, $163.9\pm 1.5$ ps. The corrected $p_{T}$ spectra in symmetric rapidity bins (-0.2,0) vs. (0,0.2) are also consistent.

The systematic uncertainty of the raw yield extraction is estimated by changing the histogram fitting method to bin counting method or by changing the fitting ranges. The maximum difference between these scenarios and the default one is considered as one standard deviation. The contribution varies by $p_{T}$, rapidity, and centrality and the overall contribution is less than 5% for the invariant yield. The systematic uncertainty in the TPC acceptance and efficiency correction $\varepsilon_{\rm TPC}$ is estimated by varying the cuts on track selection criteria  Adamczyk et al. (2017) and topological variables (for $\Xi^{-}$ only). The contribution to the total yield is 4-5% for $K^{-}$, 13-16% for $\phi$ and 6-10% for $\Xi^{-}$. This leads to a 10-13% (12-18%) uncertainty in the measured $\phi/K^{-}$ ($\phi/\Xi^{-}$) ratio. The uncertainty of the PID efficiency correction is estimated by varying the PID selection cuts and the contribution is less than 3% to the total yield. For the $p_{T}$ integrated yield, the uncertainty due to the extrapolation to the full $p_{T}$ range is estimated by choosing several fitting functions including Levy, Blast-Wave, $m_{T}$-exponential, $p_{T}$-exponential Abelev et al. (2009a), and the maximum difference between these scenarios and the default one ($m_{T}$-exponential) is quoted as one standard deviation. This contribution is 5-7% for $K^{-}$, 14-17% for $\phi$ and 13-15% for $\Xi^{-}$, respectively. This measurement covers nearly the full rapidity range from $y$=0 to the target region. The systematic uncertainty due to the rapidity coverage extrapolation is negligible compared to other systematic sources. For each individual $\phi$-meson, $K^{-}$ and $\Xi^{-}$ measurement, some of the uncertainties are correlated or partially correlated (e.g. TPC and PID). To avoid the correlation in the ratio measurement, we vary the above cuts simultaneously for $\phi$, $K^{-}$ and $\Xi^{-}$, then quote the difference in the final ratios as the systematic uncertainties.

Figure 2 and  3 show the acceptance $\times$ efficiency corrected $K^{-}$, $\phi$ and $\Xi^{-}$ invariant yields as a function of $m_{T}-m_{0}$ ($m_{T}=\sqrt{m_{0}^{2}+p_{T}^{2}/c^{2}}$, where $m_{0}$ is particle rest mass, and $c$ is the speed of light) for various rapidity ranges in 0–10%, 10–40% and 40–60% Centrality Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{3\,GeV}}$. Dashed and solid lines depict fits to the spectra with the $m_{T}$-exponential function in order to extrapolate to the unmeasured $p_{T}$ ranges ($\sim$20-40% for $K^{-}$ which vary rapidity, $\sim$33-50% for $\phi$ and $\sim$40-60% for $\Xi^{-}$). The fitted inverse slope parameters indicate harder spectra for the $\phi$-mesons compared to the $K^{-}$ and $\Xi^{-}$ within uncertainties. The inverse slope parameters gradually decrease from mid-rapidity to forward/backward rapidity and follow the $T_{\rm eff}/\cosh(y)$ distribution well. The inverse slope parameter at $y=0$, $T_{\rm eff}$, is extracted to be $177\pm 5(stat)\pm 8(sys)$ MeV for $\phi$ meson, $158\pm 3(stat)\pm 3(sys)$ MeV for $K^{-}$ and $156\pm 3(stat)\pm 24(sys)$ MeV for $\Xi^{-}$ in 0–10% central collisions. This agrees with the collision energy dependence trend from other experiments Alt et al. (2008a); Adamczewski-Musch et al. (2018). Table 1 lists the extracted $T_{\rm eff}$ parameter for these particles in different centrality bins from this measurement.

The $p_{T}$ integrated rapidity distributions $dN/dy$ are displayed in Fig. 4 for Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{3\,GeV}}$ for three different centralities. Solid curves depict Gaussian function fits to the data points with the centroid parameter fixed to zero. They are used to extrapolate to the unmeasured rapidity region ($\sim$5% for $K^{-}$, $\sim$9% for $\phi$ and $\sim$6% for $\Xi^{-}$) for calculating total multiplicities. The integral yields follow the collision energy trend from other experiments and drop quickly toward the low energies around threshold  Alt et al. (2008a); Agakishiev et al. (2009a); Piasecki et al. (2015); Gasik et al. (2016); Adamczewski-Musch et al. (2018).

The $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios are presented in Fig. 5 as a function of collision energy $\sqrt{s_{\rm NN}}$, including the midrapidity data in central Au+Au or Pb+Pb data from the AGS, SPS and RHIC BES at higher energies and $4\pi$ acceptance data from SIS at lower energies. The black solid circles show our measurements in the 0-10% centrality bin in Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{3\,GeV}}$. The measured $\phi$, $K^{-}$ and $\Xi^{-}$ yields in 4$\pi$ and the $\phi/K^{-}$, $\phi/\Xi^{-}$ ratios in different centrality bins are listed in Tab. 2. The $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios measured at 3 GeV are slightly higher than, or comparable to, the values at high energies for $\sqrt{s_{\rm NN}}\geqslant$ 5 GeV Afanasiev et al. (2002); Back et al. (2004); Alt et al. (2008a, b, c); Abelev et al. (2009b); Agakishiev et al. (2009b); Abelev et al. (2015); Adam et al. (2020) despite the collision energy being very close to the $\phi$ threshold and below the $\Xi^{-}$ threshold in NN collisions. Note that the enhancement of $\phi/K^{-}$ and $\phi/\Xi^{-}$ were also observed at lower collision energies in ${\sqrt{s_{\rm NN}}=\rm{2.4\,GeV}}$ Au+Au Adamczewski-Musch et al. (2018) and ${\sqrt{s_{\rm NN}}=\rm{2.6\,GeV}}$ Ar+KCl collisions Agakishiev et al. (2009a, b), respectively.

Various curves in Fig. 5 represent the predictions of $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios from several model calculations in central A+A collisions. Statistical model calculations, based on the Grand Canonical Ensemble and Canonical Ensemble for strangeness with several different choices of strangeness correlation length ($r_{c}$), were calculated using the THERMUS package Wheaton et al. (2009) with energy dependent freeze-out parameters (chemical freeze-out temperature $T_{\rm ch}$, baryon chemical potentials $\mu_{B}$) taken from Andronic et al. (2018), for instance, $T_{\rm ch}$ = 72.9 MeV and $\mu_{B}$ = 701.4 MeV for $\sqrt{s_{\rm NN}}$ = 3 GeV. We noted that the $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios from GCE depend on strangeness chemical potential, $\mu_{S}$. From the results of the thermal model fit to the STAR BES-I data Adam et al. (2020), there is an empirical relation $\mu_{S}$ = $\mu_{B}/4$ in the collision energy region between 7.7 - 39 GeV. The same relation was assumed and used in the GCE calculation at lower energies presented in Fig. 5. Our measured $\phi/K^{-}$ and $\phi/\Xi^{-}$ ratios are larger than this GCE calculation: $\chi^{2}$/ndf = 26.0/2 ($p$-value $<$ $1e^{-5}$), which indicates the event-by-event strangeness conservation is crucial Braun-Munzinger et al. (2004) in such collisions. The exact GCE calculation depends on the precise determination of $T_{\rm ch}$, $\mu_{B}$, $\mu_{S}$ etc, which can be extracted through a global fit to various other particle yields at 3 GeV. In the canonical approach, the correlation length, $r_{c}$, defines a region of the particle production phase space inside which the production of the strangeness is canonically conserved. Both the $\phi/K^{-}$ and $\phi/\Xi^{-}$ data from our measurement favor the CE thermodynamics for strangeness with a small strangeness correlation length ($r_{c}\sim 2.7$ fm for $\phi/K^{-}$ and $r_{c}\sim 4.2$ fm for $\phi/\Xi^{-}$). It is worthwhile to point out that the CE calculations with the same $r_{c}$ parameter cannot describe our $\phi/K^{-}$ and $\phi/\Xi^{-}$ data simultaneously. The CE calculation with $r_{c}\sim 4.2$ fm describes $\phi/\Xi^{-}$ well while it deviates by about 3.5$\sigma$ for $\phi/K^{-}$. $r_{c}$ is an approximation in the CE for reproducing the strange production in heavy-ion collisions. It is unclear if the same value of $r_{c}$ should fit for both S=1 (e.g. Kaon) and S=2 (e.g. $\Xi^{-}$). On the other hand, transport model calculations Steinheimer and Bleicher (2015); Steinberg et al. (2019) with high mass strange resonances reproduce the data implying that the feed down is relevant.

Previous measurements from smaller collision systems (Ar+KCl and Al+Al collisions) show comparable or higher $\phi/K^{-}$ and/or $\phi/\Xi^{-}$ ratios at energies below 3 GeV Agakishiev et al. (2009a, b); Piasecki et al. (2015); Gasik et al. (2016). The exclusive measurement in $p$+$p$ collisions at 2.7 GeV shows a much larger $\phi/K^{-}$ ratio ($1.04\pm 0.23$) Maeda et al. (2008), while the measured ratio at 17.3 GeV ($0.11\pm 0.01$) Aduszkiewicz et al. (2017, 2020a) is comparable to that in central Au+Au/Pb+Pb collisions at similar energies Alt et al. (2008a); Adam et al. (2020). The $\phi/\Xi^{-}$ ratio in $p$+$p$ collisions at 17.3 GeV Aduszkiewicz et al. (2020a, b), $5.09\pm 0.36$, is also significantly larger than that in central Au+Au/Pb+Pb collisions  Alt et al. (2008a, c); Adam et al. (2020). In our measurement at 3 GeV, there is no obvious difference in the $\phi/K^{-}$ ratio between the 0–10% and 10–40% central bins, while the result in the most peripheral 40–60% central bin shows a hint of a larger value, as shown in Tab. 2. Similarly, the $\phi/\Xi^{-}$ ratio in mid-central collisions seems to be larger than that in central collisions. Overall, these observations are qualitatively consistent with the expectation that a smaller canonical volume in the smaller system leads to a higher observed $\phi/K^{-}$ and/or $\phi/\Xi^{-}$ ratio.

Hadronic transport models are widely used in the high baryon density region to study the properties of the produced dense matter Bass et al. (1998); Bleicher et al. (1999); Hartnack et al. (2012); Steinheimer and Bleicher (2015); Steinberg et al. (2019); Song et al. (2021). In the modified version of the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model Steinheimer and Bleicher (2015), $\textup{UrQMD}^{2}$, new decay channels from high mass baryon resonances to $\phi$ and $\Xi^{-}$ are deployed. The relevant decay branching fraction was determined by fitting the experimental data from $p$+$p$ collisions Maeda et al. (2008). From the comparison shown in Fig. 5, the modified $\textup{UrQMD}^{2}$ calculation for central ($\rm{b}<5\,\rm{fm}$) Au+Au collisions agrees with the data points at low ${\sqrt{s_{\rm NN}}}$, including our new measurement for $\phi/K^{-}$. However calculations from the public $\textup{UrQMD}^{1}$ model Bass et al. (1998); Bleicher et al. (1999) underestimate our measurements for both $\phi/K^{-}$ and $\phi/\Xi^{-}$. The SMASH (Simulating Many Accelerated Strongly-interacting Hadrons) model Steinberg et al. (2019) attempts to incorporate the newest available experimental data to constrain the resonance branching ratios. These data include both elementary hadronic cross sections and dilepton invariant mass spectra. The $\phi/K^{-}$ ratio is reasonably reproduced using SMASH in the smaller system and ${\sqrt{s_{\rm NN}}}$ below 3 GeV, despite the overestimation of each individual ($\phi$, $K^{-}$) transverse mass spectrum measured, e.g. in Au+Au 0-40% system by HADES Adamczewski-Musch et al. (2018); Steinberg et al. (2019). The predicted $\phi/K^{-}$ ratio from the same model is about 2.5$\sigma$ higher than central Au+Au 0–10% collisions at 3 GeV. This indicates that some important in-medium mechanism for strangeness production and propagation may be missing for the large system in SMASH. Both UrQMD and SMASH calculations reproduced the measured strangeness data highlighting the importance of the contributions of the resonances in the low energies. Furthermore, the $\phi$-meson scattering with the baryonic medium remains an open question from recent measurements of $\pi$ induced nucleus reactions and the $p$-$\phi$ femtoscopy Adamczewski-Musch et al. (2019b); Acharya et al. (2021). More detailed investigations are needed in order to understand the dynamics of strange and multi-strange hadrons at low energy nuclear collisions.

Our measurement of $K^{-}$, $\phi$ and $\Xi$ production yields in 3 GeV Au+Au collisions demonstrates the necessity of the Canonical Ensemble for strangeness at low energy heavy-ion collisions. In the meantime, hadronic transport model calculations (UrQMD and SMASH) including resonance contributions reproduce the data. These observations suggest a change of the medium properties at 3 GeV compared to those from higher energy collisions. Similar conclusions have been reached from the measurements of collectivity Abdallah et al. (2021b) and high moment of protons Abdallah et al. (2021a) in 3 GeV Au+Au collisions.

In summary, we report the systematic measurements of $K^{-}$, $\phi(1020)$ and $\Xi^{-}$ production yields and the $\phi/K^{-}$, $\phi/\Xi^{-}$ ratios in Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{3\,GeV}}$ with the STAR experiment at RHIC. The measured $\phi/K^{-}$ ratio is significantly larger than the statistical model prediction based on Grand Canonical Ensemble in the 0–10% central collisions. Both the results of $\phi/K^{-}$ ($r_{c}\sim 2.7$ fm) and $\phi/\Xi^{-}$ ($r_{c}\sim 4.2$ fm) ratios favor the Canonical Ensemble model for strangeness production in such collisions. Transport models, including the resonance decays, could reasonably describe our measured $\phi/K^{-}$ ratio at 3 GeV and the increasing trend of $\phi/\Xi^{-}$ at lower energies. The new results from this paper suggest a significant change in the strangeness production for ${\sqrt{s_{\rm NN}}<\rm{5\,GeV}}$, providing new insights towards the understanding of the QCD medium properties at high baryon density.

We would like to thank K. Redlich and J. Steinheimer for fruitful discussions. We thank the RHIC Operations Group and RCF at BNL, the NERSC Center at LBNL, and the Open Science Grid consortium for providing resources and support. This work was supported in part by the Office of Nuclear Physics within the U.S. DOE Office of Science, the U.S. National Science Foundation, National Natural Science Foundation of China, Chinese Academy of Science, the Ministry of Science and Technology of China and the Chinese Ministry of Education, the Higher Education Sprout Project by Ministry of Education at NCKU, the National Research Foundation of Korea, Czech Science Foundation and Ministry of Education, Youth and Sports of the Czech Republic, Hungarian National Research, Development and Innovation Office, New National Excellency Programme of the Hungarian Ministry of Human Capacities, Department of Atomic Energy and Department of Science and Technology of the Government of India, the National Science Centre of Poland, the Ministry of Science, Education and Sports of the Republic of Croatia, German Bundesministerium für Bildung, Wissenschaft, Forschung and Technologie (BMBF), Helmholtz Association, Ministry of Education, Culture, Sports, Science, and Technology (MEXT) and Japan Society for the Promotion of Science (JSPS).