Evidence of isospin-symmetry violation in high-energy collisions of atomic nuclei CERN-EP-2023-283 Strong interactions preserve an approximate isospin symmetry between up ($u$) and down ($d$) quarks, part of the more general flavor symmetry. In the case of $K$ meson production, if this isospin symmetry were exact, it would result in equal numbers of charged ($K^{+}$ and $K^{-}$) and neutral ($K^{0}$ and $\overline{K}^{\,0}$) mesons produced in collisions of isospin-symmetric atomic nuclei. Here, we report results on the relative abundance of charged over neutral $K$ meson production in argon and scandium nuclei collisions at a center-of-mass energy of 11.9 GeV per nucleon pair. We find that the production of $K^{+}$ and $K^{-}$ mesons at mid-rapidity is $(18.4\pm 6.1)\%$ higher than that of the neutral $K$ mesons. Although with large uncertainties, earlier data on nucleus-nucleus collisions in the collision center-of-mass energy range $2.6<\sqrt{s_{NN}}<200$ GeV are consistent with the present result. Using well-established models for hadron production, we demonstrate that known isospin-symmetry breaking effects and the initial nuclei containing more neutrons than protons lead only to a small (few percent) deviation of the charged-to-neutral kaon ratio from unity at high energies. Thus, they cannot explain the measurements. The significance of the flavor-symmetry violation beyond the known effects is 4.7$\sigma$ when the compilation of world data with uncertainties quoted by the experiments is used. New systematic, high-precision measurements and theoretical efforts are needed to establish the origin of the observed large isospin-symmetry breaking.

1 Introduction

One of the main aims of basic research is to understand the fundamental constituents of matter and the interactions between them. Within Quantum Chromodynamics (QCD) [1], the theory of strong interactions, the fundamental particles are quarks and gluons carrying color – the charge of strong interactions. Because of confinement, quarks and gluons are hidden in colorless hadrons, particularly protons and neutrons. The strong force binds them, forming atomic nuclei.

Accelerator-based experiments recording collisions of highly energetic hadrons and nuclei allow for systematic studies of the properties of strong interactions. In these collisions, many new particles are produced. They are predominantly mesons containing one valence quark ($\mathit{q}$) and one valence anti-quark ($\mathit{\overline{q}}$). The most copiously produced are the lightest mesons, pions and kaons, built from up ($\mathit{u}$), down ($\mathit{d}$) and strange ($\mathit{s}$) quarks and the corresponding anti-quarks.

QCD assumes that interactions are independent of quark type (flavor) for equal quark masses and in the absence of other interactions, a feature known as flavor symmetry. When only the light quarks up and down are considered, flavor symmetry reduces to isospin symmetry, historically introduced in the pre-QCD period by Heisenberg to understand the properties of nuclei [2]. The masses of up and down quarks, $m_{\mathit{u}}=2.16\pm 0.07$ MeV and $m_{\mathit{d}}=4.70\pm 0.07$ MeV [3], are not equal, but they are much smaller than the QCD scale, $\Lambda_{QCD}$ [4, 5]. (Note that units in this paper follow the Particle Data Group (PDG) [3] convention: masses and energies are expressed in MeV (or GeV), whereas momenta in $\mbox{Me\kern-1.00006ptV}/c$ (or $\mbox{Ge\kern-1.00006ptV}\!/\!c$). The relative differences are given as the ratio of the difference to the mean.) Hence isospin-symmetry breaking effects are small, as confirmed by the mass ratios of pions and kaons, $\left(m_{\pi^{+}}-m_{\pi^{0}}\right)/\left(m_{\pi^{+}}+m_{\pi^{0}}\right)\simeq 0.017$ and $\left(m_{K^{+}}-m_{K^{0}}\right)/\left(m_{K^{+}}+m_{K^{0}}\right)\simeq-0.004$. Moreover, the elastic cross sections for pion-pion, pion-nucleon, and nucleon-nucleon scattering closely follow the predictions of isospin symmetry [6, 7]. Here, of special interest is a specific isospin transformation, an inversion of the third component of the isospin, called the charge transformation for historical reasons. It is equivalent to swapping $u\leftrightarrow d$ quarks. At the hadronic level, the charge transformation implies swapping $p\leftrightarrow n$, $\pi^{+}\leftrightarrow\pi^{-}$, $K^{+}\leftrightarrow K^{0}$, $\overline{K}^{\,0}\leftrightarrow K^{-}$, etc.

Let us consider nucleus-nucleus ($A$+$A$) collisions and, for simplicity, assume that both nuclei have an equal number of protons and neutrons. Without referring to a detailed mathematical formalism, charge symmetry means that strong interactions are invariant under the charge transformation of every nucleus and hadron of the initial and final states. For an ensemble of initial states being invariant under the charge transformation, the probabilities of having initial states related by this transformation are equal. This is indeed the case of nucleus-nucleus collisions, for which each nucleus has an equal number of protons and neutrons. Then, the invariance under charge transformation also holds for the final state ensemble, implying that the mean multiplicities of charge-transformation related hadrons, such as $K^{+}$ and $K^{0}$ as well as $\overline{K}^{\,0}$ and $K^{-}$, coincide:

$$\langle K^{+}\rangle=\langle K^{0}\rangle\text{ and }\langle K^{-}\rangle=\langle\overline{K}^{\,0}\rangle\text{ . }$$ (1)

Note that since models predict only properties of ensembles of events and not outcomes of single events, we need to consider quantities averaged over the event ensembles. The subject has a vast literature; see, for example Refs. [8, 9, 10, 11, 12, 13]. Consequently, the exact isospin symmetry prediction for the charged-to-neutral kaon ratio in nucleus-nucleus collisions with electric charge to baryon number $Q/B=1/2$ reads

$$R_{K}\equiv\frac{\langle K^{+}\rangle+\langle K^{-}\rangle}{\langle K^{0}\rangle+\langle\overline{K}^{\,0}\rangle}=\frac{\langle K^{+}\rangle+\langle K^{-}\rangle}{2\langle K_{S}^{0}\rangle}=1\leavevmode\nobreak\ .$$ (2)

(Note that the $\mathit{K}^{\,0}$ and $\mathit{\overline{K}}^{\,0}$ states are produced in strong interactions, but they decay through weak interactions. Consequently, in the final state, one observes linear combinations of the latter known as $\mathit{K}^{0}$ short ($\mathit{K}^{0}_{S}$) and $\mathit{K}^{0}$ long ($\mathit{K}^{0}_{L}$), where short and long refer to their weak decay lifetime [3]. By neglecting the small CP violation, the multiplicities corresponding to weak and strong eigenstates are related by $\langle K_{S}^{0}\rangle=\frac{1}{2}\langle K^{0}\rangle+\frac{1}{2}\langle\overline{K}^{\,0}\rangle=\langle K_{L}^{0}\rangle$.) The prediction given by Eq. (2) is a reference for experimental testing of the isospin symmetry in hadron production processes. For a more detailed introduction and didactic derivations see Ref. [14].

Here, we report a measurement of the ratio $R_{K}$ in the 10% most central collisions of argon ($\mathrm{Ar}$) and scandium ($\mathrm{Sc}$) nuclei at center-of-mass energy per nucleon pair equal to $\sqrt{s_{NN}}$ $=$ 11.9 GeV. Further on, we compare the NA61/SHINE result with the world data on charged and neutral kaon production in nucleus-nucleus collisions. The results indicate a significant excess of charged over neutral kaon production. This excess cannot be explained by known effects violating the isospin symmetry. This is discussed and demonstrated by comparing experimental results to well-known theoretical approaches, the statistical Hadron Resonance Gas (HRG) [15] and the dynamical Ultrarelativistic Quantum Molecular Dynamics (UrQMD) [16] models. The predictions of models are calculated for reactions corresponding to experimental data, generally with $Q/B<1/2$. They consider isospin-breaking effects in strong interactions and, importantly, the production and subsequent decays of the $\phi$ mesons.

Summarizing, the NA61/SHINE Collaboration measures a charged-to-neutral kaon ratio $R_{K}=1.184\pm 0.061$ in Ar+Sc collisions at 11.9 GeV per nucleon pair. This value aligns with previous experimental measurements, albeit their uncertainties are larger. The significance of the isospin symmetry violation beyond the known effects amounts to 4.7$\sigma$ when all measurements are considered, and uncertainties quoted by the experiments are used. This is the first evidence of an unexplained isospin symmetry violation in hadron production processes.

2 Results

2.1 Production of $\boldsymbol{K}$ mesons in central Ar+Sc collisions at the CERN SPS

The new experimental results presented here have been obtained by the NA61/SHINE fixed-target experiment at the CERN Super Proton Synchrotron [17]. The measurements of $\mathit{K^{+}}$ and $\mathit{K^{-}}$ production in the 10% most central Ar+Sc reactions at $\sqrt{s_{NN}}=11.9$ GeV have been published elsewhere [18]. The analysis procedure and details of systematic uncertainties are given in Refs. [19, 20]. Here, we present the first measurement of $\mathit{K}^{0}_{S}$ production in nucleus-nucleus collisions from NA61/SHINE. Earlier data from this experiment, on $\mathit{K}^{0}_{S}$ production in $p$+C, $\pi^{-}$+C, $\pi^{+}$+C, $\pi^{+}$+Be, and $p$+$p$ collisions can be found in Refs. [21, 22, 23, 24, 25, 26]. For more details concerning the experimental procedure, see Methods’ "Experimental procedure" subsection.

The comparison of the rapidity distribution of $\mathit{K}^{0}_{S}$ mesons to the average of rapidity distributions for $\mathit{K^{+}}$ and $\mathit{K^{-}}$ mesons is presented in Fig. 1. The rapidity $y$ is a relativistic generalization of the particle velocity along the direction of the incoming nuclei. We calculate rapidity in the nucleon-nucleon collision center-of-mass system, and positive $y$ corresponds to the direction of the Ar nucleus.

In the entire range of rapidity covered by the measurement, the averaged charged $\mathit{K}$ mesons yield prevails significantly over the neutral $\mathit{K}^{0}_{S}$ mesons one. To quantify this effect, Table 1 presents the rapidity densities $\mathrm{d}n/\mathrm{d}{y}$ of $\mathit{K^{+}}$, $\mathit{K^{-}}$ and $\mathit{K}^{0}_{S}$ production measured at mid-rapidity ($y\approx 0$). Here, the relative excess of charged mesons is $\mathrm{(18.4\pm 6.1)}$%. Integration of the two distributions in Fig. 1 over positive rapidity, $y>0$, gives $4.28\pm 0.13$ and $3.22\pm 0.37$ for the production rates per collision of $\mathit{(K^{+}+K^{-})/2}$ and $\mathit{K}^{0}_{S}$, respectively (total uncertainties are given; the quantities provided for charged $\mathit{K}$ mesons are based on Ref. [18]). The resulting difference of $1.06\pm 0.39$ corresponds to a surplus of charged ($\mathit{K^{+}}$ and $\mathit{K^{-}}$) over neutral ($\mathit{K}^{0}$ and $\mathit{\overline{K}}^{\,0}$) states equal to $2.12\pm 0.79$ at positive rapidity. Under the assumption that the charged-to-neutral ratio would be similar also at negative rapidity, the total excess would amount to $4.2\pm 1.6$ additional $K^{+}$ or $K^{-}$ mesons per one central Ar+Sc collision.

Figure 1: Comparison of rapidity spectrum of neutral ($\boldsymbol{K^{0}_{S}}$) with the averaged spectrum of charged ($\boldsymbol{K^{+}}$ and $\boldsymbol{K^{-}}$) mesons in the 10% most central Ar+Sc collisions at $\sqrt{s_{NN}}=11.9$ GeV. Total uncertainties, calculated as the square root of the sum of squared statistical and systematic uncertainties ($\sqrt{\sigma_{stat}^{2}+\sigma_{sys}^{2}}$) are drawn. For charged $K$ mesons, the total uncertainties were calculated separately for positively and negatively charged and then propagated.

			statistical	systematic	total
	$\left(\frac{\mathrm{d}n}{\mathrm{d}y}\right)_{y\hskip 0.85358pt\approx\hskip 0.85358pt0}\mathit{(K^{+})}$	3.732	$\pm$ 0.016	$\pm$ 0.148	$\pm$ 0.149
	$\left(\frac{\mathrm{d}n}{\mathrm{d}y}\right)_{y\hskip 0.85358pt\approx\hskip 0.85358pt0}\mathit{(K^{-})}$	2.029	$\pm$ 0.012	$\pm$ 0.069	$\pm$ 0.070
	$\left(\frac{\mathrm{d}n}{\mathrm{d}y}\right)_{y\hskip 0.85358pt\approx\hskip 0.85358pt0}\mathit{(K}^{0}_{S}{\mathit{)}}$	2.433	$\pm$ 0.027	$\pm$ 0.102	$\pm$ 0.106
	charged-to-neutral $\mathit{K}$ meson ratio:
	${R_{K}}$	1.184	$\pm$ 0.014	$\pm$ 0.060	$\pm$ 0.061

Table 1: Rapidity densities of charged and neutral $\boldsymbol{K}$ mesons produced at mid-rapidity. The measurement was performed in the 10% most central Ar+Sc collisions at $\sqrt{s_{NN}}=11.9$ GeV, as described in Methods’ "Experimental procedure" subsection. The excess of charged over neutral mesons is quantified by the ratio ${R_{K}}$ defined in Eq. (2).

A comparison of distributions of $\mathit{K}^{0}_{S}$ with averaged $\mathit{K^{+}}$ and $\mathit{K^{-}}$ mesons as a function of transverse momentum $p_{T}$ (the momentum component perpendicular to the direction of the incoming nuclei) is shown in Fig. 2. Both distributions are integrated over the rapidity range $0<y<2$. The prevalence of charged over neutral $\mathit{K}$ mesons is again evident. The insert in the figure shows the $p_{T}$-dependence of the ratio ${R_{K}}$. The corresponding excess of $\mathit{K}$ mesons containing $\mathit{u}$, $\mathit{\overline{u}}$ over those containing $\mathit{d}$, $\mathit{\overline{d}}$ quarks and anti-quarks remains in the range 6–33% over the considered range of $p_{T}$.

Figure 2: Comparison of transverse momentum spectrum of neutral ($\boldsymbol{K^{0}_{S}}$) with the averaged spectrum of charged ($\boldsymbol{K^{+}}$ and $\boldsymbol{K^{-}}$) mesons in the 10% most central Ar+Sc collisions at $\sqrt{s_{NN}}=11.9$ GeV. The bottom panel shows the ratio of the two distributions, as defined in Eq. (2). The meaning of the total uncertainties drawn is the same as in Fig. 1.

2.2 Comparison to the world data and models

Figure 3 compares the present measurement of the ratio ${R_{K}}$ at mid-rapidity and the world data compiled by us and detailed in Methods’ "World data" subsection. The experimental results were obtained by CERES [27, 28, 29], STAR BES [30, 31], STAR [32, 33, 34, 35], ALICE [36, 37], NA35 [38, 39], NA49 [40, 41, 42], and HADES [43, 44] experiments. We note that the compilation includes measurements at mid-rapidity and total multiplicities. This may increase the overall spread between the data points. We also note the sizeable uncertainties of the earlier measurements. These probably explain why the aforementioned charged-over-neutral anomaly was never reported as an experimental observation. Despite these uncertainties, a consistent picture emerges in the energy range $2.6<\sqrt{s_{NN}}<200$ GeV. The ratio is above one for all experiments except NA35 and ALICE.

Figure 3: The charged-to-neutral kaon ratio $\boldsymbol{R_{K}}$ as a function of collision energy. The symbols show the experimental world data with total uncertainties; see Methods’ "World data" subsection for details. The black line shows the HRG predictions for $Q/B=0.4$. The black dots indicate the HRG predictions for $Q/B$ values corresponding to the ones in the experiments. For different nuclei, $Q/B$ corresponds to the electric charge over the baryon number of the whole system. The gray squares show UrQMD predictions. See Methods’ "Models" subsection for details on models.

Figure 4: The experimental data for the charged-to-neutral kaon ratio divided by the HRG baseline $\boldsymbol{R_{K}/R_{K}^{HRG}}$ as a function of collision energy. The symbols are explained in Fig. 3. The solid black line shows the weighted average of the experimental data, and the shaded area shows the uncertainty of the weighted average.

Figure 3 also compares the data with the HRG and UrQMD model predictions. The HRG calculations were performed with $Q/B=0.4$ (solid line) and $Q/B$ values corresponding to collisions studied in the experiments (black dots). The UrQMD results were obtained for central Au+Au collisions. See Methods’ "Models" subsection for details on the models. The predictions of models agree well with each other but are systematically lower than the experimental data. At energies larger than 10 GeV, the mass difference between charged and neutral kaons, leading to isospin-symmetry breaking (mostly via $\phi$-meson decays), increases the $R_{K}$ ratio by about 0.03. Other isospin-breaking effects can be neglected. The ratio $R_{K}$ is reduced for collisions with $Q/B<1/2$, but this is insignificant for energies larger than 10 GeV. At energies lower than 10 GeV, the $R_{K}$ ratio is significantly more sensitive to the known isospin-breaking effects as well as to the $Q/B$ ratio; see Methods’ "Models" subsection for more details.

To quantify the measured isospin-symmetry breaking beyond the known effects, the ratio of the measured $R_{K}$ to the corresponding HRG baseline is shown in Fig. 4. We do not consider the lowest energy data point (HADES) because the known isospin-breaking and $Q/B<1/2$ effects are significant at the low collision energies; thus, predictions may be model-dependent. We also do not consider the NA35 point because, unlike other measurements, charged kaons were identified by reconstructing their decays, leading to large statistical uncertainties and possible biases. Thus, the number of selected measurements at different collision energies (from SPS to LHC) is $13$. Out of them, only one is below unity. No significant dependence of the double ratio on collision energy and nuclear mass number of colliding nuclei is visible.

The weighted average of all double ratios shown in Fig. 4 is $1.129\pm 0.027$, where the uncertainty was calculated using kaon uncertainties reported by experiments. The HRG uncertainties are small and were neglected. The significance of the isospin violation is $4.7\sigma$. The $\chi_{\textit{m}in}^{2}/\text{dof}\approx 0.3$ may indicate either a correlation between results or an overestimation of the uncertainties.

3 Discussion

In the following, we discuss possible effects that may potentially contribute to the violation of isospin symmetry in kaon production.

First, we consider symmetry-breaking effects due to the non-equal bare $u$ and $d$ masses in strong interactions. They are included in the HRG and UrQMD models. Then, we discuss the possible influence related to electromagnetic and weak processes.

1. Mass effects within strong interactions. Within QCD, the isospin symmetry is not exact because $u$ and $d$ quark masses are different, $\approx 2.2$ and $\approx 4.7$ MeV, respectively. The different quark masses lead to different masses of hadrons within the isospin multiplets, particularly different masses of charged and neutral kaons. This modifies the ratio $R_{K}$. We list below the considered effects and quantify their influence on $R_{K}$ using the statistical Hadron Resonance Gas model (HRG) [15]. The results are cross-checked with the microscopic transport model, UrQMD [45, 16, 46]. For details, see Methods’ "Models" subsection.

   1. Smaller masses of charged kaons than neutral ones, $m_{K^{+}}=m_{K^{-}}=493.7$ MeV and $m_{K^{0}}=m_{\overline{K}^{\,0}}=497.6$ MeV, lead to an increase of $R_{K}$ resulting from direct kaon production by about 0.02. This was estimated by removing resonances from the particle list of HRG. We have numerically verified that HRG for $Q/B=1/2$ and with exact isospin symmetry gives $R_{K}=1$, as expected.

   2. A significant fraction of kaons (up to several dozen percent) results from resonance decays [47]. Different kaon masses affect the branching ratios of resonances. The most striking example is $\phi(1020)$ meson, which decays about $1.45$ more frequently into charged kaons than neutral ones. This large difference is because the $\phi(1020)$ mass is just above the kaon-kaon thresholds. Including the kaon production from resonance decays increases $R_{K}$ by about 0.03.

   3. In connection to the previous point, other potentially relevant $K\overline{K}$ decays refer to the resonances $a_{0}(980)$ and $f_{0}(980)$, whose masses are close to the $K\overline{K}$ decay threshold. In the PDG review [3], for both resonances $K\overline{K}$ decays are reported as seen. The HRG model used in this work initially assumes equal branching ratios ($BR$) of $K^{+}K^{-}$ and $K^{0}\overline{K}^{\,0}$ decay channels. Yet, just as for the $\phi$ meson, isospin breaking may be relevant. (Note that isospin breaking for $a_{0}(980)$ and $f_{0}(980)$ is experimentally confirmed by the $a_{0}(980)$ - $f_{0}(980)$ mixing [48, 49] leading to the otherwise forbidden decays $a_{0}^{0}(980)\rightarrow f_{0}(980)\rightarrow\pi\pi$ and $f_{0}(980)\rightarrow a_{0}^{0}(980)\rightarrow\pi\eta$; see the experimental results for these small but non-zero transitions in Ref. [50].) Using the estimate (see Methods’ "Models" subsection) $BR(K^{+}K^{-})$/$BR(K^{0}\overline{K}^{\,0})\approx 1.2$ one gets an increase of $R_{K}$ by about $0.5\%$ . The yield of charged kaons may also be affected by other effects, such as the electromagnetic interaction between $K^{+}$ and $K^{-}$ (see also (B) (ii) below). To calculate the upper limit due to this effect, we assumed $BR(K^{+}K^{-})$ $=$ $BR(K\overline{K})$ (no decays to neutral kaons) for both $a_{0}^{0}(980)$ and $f_{0}(980)$. This leads to the increase of $R_{K}$ by at most $3\%$ at the highest collision energy.

   4. Mass differences of hadrons from other isospin multiplets also break the flavor symmetry and affect $R_{K}$. The largest effect comes from the mass difference between proton and neutron, which reduces $R_{K}$ at the lowest collision energies. At energies larger than 10 GeV, this effect is negligible.

2. Electromagnetic processes. Electromagnetic interactions do not obey isospin symmetry because electric charges differ for the quark flavors $u$ and $d$. The electromagnetic interaction slightly affects the masses of hadrons. For instance, the neutral pion is lighter than the charged ones. The HRG and UrQMD include such effects since the physical masses are used. However, the effects mentioned below are not included in the models.

   1. Electromagnetic decays of hadrons are typically suppressed by a factor $\alpha\simeq 1/137$ compared to strong ones. Consequently, decays that involve the production of virtual photons and their subsequent decay into charged kaons are suppressed by a factor $\alpha^{2}$ and thus negligible. Taking into account the charge of the nuclei $Z_{1}$ and $Z_{2}$, one would expect an effect of the type $Z_{1}Z_{2}\alpha^{2}$, which is not observed in the experimental data – isospin-symmetry breaking for collisions of light and heavy nuclei is similar, see Fig. 4.

   2. There may also be non-perturbative electromagnetic effects at the hadronic level that might affect the kaon multiplicities. One notable example is the case of $K^{+}K^{-}$ pairs with low momenta of kaons $p\lesssim m_{K}\alpha\simeq 3.6$ $\mbox{Me\kern-1.00006ptV}/c$, see e.g. Ref. [51]. This is, in particular, the case for the decays of the resonances $a_{0}(980)$ and $f_{0}(980)$, whose masses are close to the two-kaon decay threshold. It is then possible that this effect will increase the number of produced $K^{+}K^{-}$ pairs. Experimental search for $K^{+}K^{-}$ interaction close to the threshold was conducted at COSY in the study of the reaction $p+p\rightarrow p+p+K^{+}+K^{-}$, see e.g. Ref. [52] and the review [53]. An increase in the cross-section close to the threshold has been measured, where different effects appear to be relevant, with an important role played by the $pK^{-}$ Coulomb attraction. We have quantified the impact of $a_{0}(980)$ and $f_{0}(980)$ decays in point (A) (iii) above and in Methods’ "Models" subsection, showing that they cannot explain the measured value of $R_{K}$.

   3. The $u\overline{u}$ and $d\overline{d}$ pair creation in strong processes may be affected by electromagnetic interactions. They are different for $u\overline{u}$ and $d\overline{d}$ pairs due to different electric charges of up and down quarks. This leads to a different phase space for their production, favoring $u\overline{u}$ pairs and thus charged kaons. In particular, the quark-gluon effective coupling is enhanced by QED effects due to the attraction among quarks, leading to a larger coupling of gluons to $u$-quarks than $d$-quarks [54]. A model of the space-time evolution of the pair creation will be needed to quantify the effect. In addition, the isospin breaking due to the Coulomb potential of highly-charged fireballs formed in heavy-ion collisions is discussed in Ref. [55] within the statistical QGP model. We recall that electromagnetic interactions are expected to modify fusion rates in the Big Bang nucleosynthesis epoch; see, for example, Refs. [56, 57].

3. Uncertainties in weak decays. The weak interaction does not obey the isospin symmetry. The mean lifetimes [3] of charged and neutral kaons are $\tau(K^{+})=\tau(K^{-})=(1.2380\pm 0.0020)\cdot 10^{-8}$ s and $\tau(K^{0}_{S})=(8.954\pm 0.004)\cdot 10^{-11}$ s, $\tau(K^{0}_{L})=(5.116\pm 0.021)\cdot 10^{-8}$ s. The charged kaons are typically measured by reconstructing their trajectories in a detector. Due to the large mean lifetime, the corrections for their losses caused by weak decays are small. In contrast, the neutral $K^{0}_{S}$ kaons are measured by reconstructing their decays into two charged pions. Typically, the corrections for the losses caused by weak decays are large. This is because the decay should be far enough from the interaction point to separate the decay point from the background in high-multiplicity $A$+$A$ collisions. Assuming that the $K^{0}_{S}$ meson is measurable when the lifetime of a particle in its rest frame is larger than the mean lifetime (typical for NA61/SHINE), one estimates the maximum relative bias of the mean multiplicity, $\Delta(\langle K_{S}^{0}\rangle)/\langle K_{S}^{0}\rangle$, as:

   $$\frac{\Delta(\langle K_{S}^{0}\rangle)}{\langle K_{S}^{0}\rangle}=\frac{3\cdot\sigma(\tau(K_{S}^{0}))}{\tau(K_{S}^{0})}\approx 0.0013\,,$$ (3)

   where $\sigma(\tau(K_{S}^{0}))$ is the uncertainty of the mean $K^{0}_{S}$ lifetime. Thus, the maximum deviation of $R_{K}$ from unity due to the uncertainty of the mean $K^{0}_{S}$ lifetime is 0.13%.

Finally, we discuss the consequences of having collisions with $Q/B<1/2$, which corresponds to many experimental results presented in Fig. 3. The third component of isospin equals $|I_{z}|=|B/2-Q|$ and therefore the total isospin is limited as $|I_{z}|\leq I\leq B/2$. The compiled experimental results in nucleus-nucleus collisions correspond to the $Q/B$ ranging from about 0.4 (Pb+Pb and Au+Au collisions) to 0.5 (S+S collisions); see Fig. 3. These limits correspond to the normalized per baryon third component and total isospin $|I_{z}|/B=0.1$, $0.1<I/B<1/2$ and $|I_{z}|/B=I/B=0$, respectively. The non-zero $I_{z}$ and $I$ for heavier nuclei can affect the charged-to-neutral kaon ratio in two ways:

1. The larger fraction of neutrons than protons for heavy nuclei having $|I_{z}|/B\approx 0.1$ enhances neutral kaon production compared to charged ones and thus reduces $R_{K}$. This fact is taken into account by the employed theoretical models that use the physical value for $Q/B$. The reduction of $R_{K}$ is significant at low collision energies and is small compared to other effects at energies larger than 10 GeV; see Fig. 3 and Methods’ "Models" subsection.

2. For $|I_{z}|/B\approx 0.1$, the total isospin is limited by $0.1\leq I/B\leq 1/2$. Generally, nuclei in the ground state have the lowest possible value of the total isospin [58]. This rule extends to a state of two identical nuclei in the ground state, which, for the considered case, implies $I/B\approx|I_{z}|/B\approx 0.1$. Thus, a possible dependence of $R_{K}$ on $I$ and $I_{z}$ reduces to the dependence on $I=|I_{z}|$. The latter is discussed above in point (a). The experimental data also allow a rough estimate of the influence of the small but non-zero value of the normalized isospin on the charged-to-neutral kaon ratio. Results for heavy nuclei, ²⁰⁸Pb+²⁰⁸Pb and ¹⁹⁷Au+¹⁹⁷Au ($R_{K}\approx 1.15$ for $I/B=|I_{z}|/B\approx 0.1$), are similar to those for intermediate nuclei, ⁴⁰Ar+⁴⁵Sc ($I/B=|I_{z}|/B\approx 0.04$), and ³²S+³²S ($I/B=|I_{z}|/B=0$), see Fig. 3. This suggests that the sensitivity of $R_{K}$ to the $I/B=|I_{z}|/B$ close to zero ($I/B=|I_{z}|/B\leq 0.1$) is low. Note, however, that there are large uncertainties in the experimental results.

3. In the case of central collisions of heavy nuclei, the charge-to-baryon ratio of the interacting nucleon system (participant nucleons) can be higher than the total proton-to-nucleon ratio in colliding nuclei. This is because protons tend to be distributed closer to the center of a nucleus [59]. However, this effect is expected to be small because the ratio $R_{K}$ seems independent of the colliding nuclei size. Low-mass nuclei can be described as alpha clusters [60] (clusters of two protons and two neutrons). Thus, if significant, the effect should disappear for collisions of low-mass nuclei, particularly at the lowest collision energy ($\lesssim 4$ GeV), where $R_{K}$ is sensitive to the small changes of the electric charge to baryon number ratio, $Q/B$. The data do not support this.

Closing comments on future perspectives of experimental and theoretical efforts are in order here.

1. Concerning measurements, reviewing the validity of the past results and confirming them with new high-precision data is important. Systematic results on the collision energy and nuclear mass dependence of the isospin-breaking effect should help us understand its nature. Measurements of the charged-to-neutral kaon ratio in collisions of an equal number of protons and neutrons would reduce uncertainty in its interpretation. The NA61/SHINE experiment plans to perform such measurements for O+O and Mg+Mg collisions [61]. If needed, the measurements of deuteron-deuteron interactions may be possible in the long future. They will require the production of primary deuteron beams in the CERN accelerator complex and the construction of a liquid deuterium target for NA61/SHINE.

2. In connection with the previous point, an interesting experimental test is also possible by considering $\pi^{-}$+C and $\pi^{+}$+C interactions [62]. While the ensemble of only one of them is not invariant under charge transformation, the ensemble having an equal number of $\pi^{-}$+C and $\pi^{+}$+C interactions is invariant. Thus, for the joint ensemble, the exact charge symmetry predicts $R_{K}=1$. NA61/SHINE recorded data on $\pi^{-}$+C and $\pi^{+}$+C interactions at 158 $\mbox{Ge\kern-1.00006ptV}\!/\!c$ in October 2024 for the test.

3. During the current CERN accelerator complex operation period (Run 3), NA61/SHINE records high-statistics data on inelastic Pb+Pb collisions at 150$A\,\mbox{Ge\kern-1.00006ptV}\!/\!c$. This will allow precision measurements of the $R_{K}$ ratio as a function of collision centrality. The larger number of neutrons than protons in Pb nuclei may reduce the ratio $R_{K}$. This should be taken into account when interpreting the results.

4. The relationship between the production rates of charged and neutral kaons in hadronic collisions was used to estimate the flux of secondary neutral kaons produced in $p$+Be collisions at beam momenta of several 100 $\mbox{Ge\kern-1.00006ptV}\!/\!c$ in neutral kaon and neutrino beams [63]. Recently, the charged-to-neutral kaon ratio was studied in $p$+$p$ interactions at SPS energies, and no excess in comparison to a simple quark counting model emerged [64].

5. Kaons play a special role due to their simple isospin structure and easy measurement. This explains why the first results on a large isospin-symmetry breaking in multi-particle production are reported for kaons. Yet, it is important to perform a similar study for other isospin multiplets in the future. For example, using the same methods, we have checked that the (anti-)proton to (anti-)neutron ratio is even less sensitive to the known isospin-symmetry breaking effects and, thus, is predicted to be almost exactly one in nucleus-nucleus collisions with $Q/B=0.5$ in a broad range of collision energies including the low energies.

6. One can extend the current models by introducing new isospin-breaking processes and fitting their parameters to the data. This can be done either for the quark-gluon processes or the hadron-resonance processes. For example, within the statistical hadronization models, one can introduce the quark fugacity factors for $u$ and $d$ quarks separately [65, 66, 67]. This could allow us to make predictions for other hadron ratios but will not explain the origin of the violation.

7. The possibility of having a phase of strongly interacting matter with a significant isospin violation was suggested by Pisarski and Wilczek within a linear $\sigma$ model of QCD [68]. They expect masses of $\pi^{0}$ and $\eta$ mesons to decrease if the $U_{A}(1)$ symmetry is effectively restored at temperatures lower than the one of the chiral phase transition. The chiral anomaly [69, 70] could break isospin (especially, it affects the pion isotriplet) but does not affect the charged-to-neutral kaon ratio.

8. Creating Disoriented-Chiral-Condensate (DCC) domains in heavy-ion collisions has been considered for many years [71, 72, 73]. They may be signaled by large fluctuations of the charged-to-neutral pion [74] and kaon ratios [75, 76]. A puzzling result on kaon fluctuations was recently reported by ALICE at LHC [77]. Its possible interpretation by the DCC or disoriented-isospin-condensates formation is discussed in Refs. [78, 79]. The considered models for the charged-symmetric ensemble of collisions predict $R_{K}=1$ [80]. The inclusion of isospin-breaking effects in the extended Linear Sigma model [81] was recently discussed in Ref. [82], where through a fit to available experimental masses and decays of light mesons, it is shown that the relative difference between the $u\overline{u}$ and $d\overline{d}$ chiral condensates amounts to $0.02\%$, implying that only very small deviations from $R_{K}=1$ are expected from this effect.

Thus, the presented results on the charged-to-neutral kaon ratio are the first evidence of an unexplained isospin symmetry violation in hadron production processes. Further studies are needed to understand the underlying physics, particularly reducing the experimental uncertainties and quantifying the role of electromagnetic effects. If these steps do not solve the issue, more speculative explanations shall be investigated.

Methods

A. Experimental procedure

Experimental setup. The SPS Heavy Ion and Neutrino Experiment (SHINE) is a fixed-target detector operating at the CERN Super Proton Synchrotron (SPS). It is a multi-purpose spectrometer optimized to study hadron production in various collisions (hadron-proton, hadron-nucleus, and nucleus-nucleus). The detection setup used for the measurements reported here is described below. Its details and a description of the detector performance can be found in Ref. [17].

The beamline is equipped with an array of beam detectors upstream and downstream of the target, used to identify and measure the trajectory of the beam particles and trigger the spectrometer data acquisition. The tracking devices of the NA61/SHINE spectrometer are Time Projection Chambers (TPCs). Two Vertex TPCs are placed inside a magnetic field. Two large-volume Main TPCs measure the charged particle trajectories downstream of the 4.5 Tm magnetic field. The latter provides the bending power for a precise determination of particle momenta. The information about the energy losses (${\mathrm{d}}E\!/\!{\mathrm{d}}x$) of the charged particles in the TPCs, together with Time-of-Flight (ToF) measurements, allows for particle identification in a wide momentum range. The most downstream detector on the beamline is the Projectile Spectator Detector (PSD). It measures the energy of the spectator remnant of the projectile nucleus, closely related to the collision centrality in nucleus-nucleus reactions.

Physics objects. This article compares the production of charged and neutral $\mathit{K}$ mesons in Ar+Sc collisions at a center-of-mass energy per nucleon pair of 11.9 GeV. The ${}^{40}_{18}\mathrm{Ar}$ beam had a momentum of 75$A\,\mbox{Ge\kern-1.00006ptV}\!/\!c$. The stationary target consisted of six ${}^{45}_{21}\mathrm{Sc}$ plates, with a total thickness of 6 mm.

A detailed account on the extraction of charged ($\mathit{K^{+}}$, $\mathit{K^{-}}$) yields can be found in Ref. [18]. Only the neutral $K^{0}_{S}$ mesons are considered in the present analysis. They can be detected via their weak decay into two charged pions $(K^{0}_{S}\rightarrow\pi^{+}\pi^{-})$. The mean lifetime $(c\tau)$ for this decay is 2.7 cm. A detailed presentation of the $K^{0}_{S}$ analysis procedure and systematic uncertainties can be found in Ref. [83].

Analysis. Before analyzing $K^{0}_{S}$ mesons, the recorded Ar+Sc collision data undergo event and track selection procedures. Event selection uses information from the beam detectors to ensure the quality of the measured beam trajectory. It rejects events with more than one beam-target interaction during the trigger-time window. It also reduces the background from off-target interactions based on information about the quality of the main interaction vertex. Finally, it selects the 10% most central collisions using the information from the PSD. This is realized by selecting the 10% lowest energy deposits from the spectator remnant of the Ar nucleus. The total number of recorded collisions (events) was $2.77\cdot 10^{6}$, from which $1.03\cdot 10^{6}$ (37%) remained after all cuts. For more details, especially the centrality selection, see Ref. [18].

The next step is reconstructing the charged particle tracks in the TPCs. Pattern recognition algorithms combine space points recorded in the TPCs into tracks. Their curvature and the magnetic field are used to compute the momenta of the corresponding particles. The minimum number of reconstructed space points in the VTPCs must be more than 10, and the computed momenta must be larger than 400 $\mbox{Me\kern-1.00006ptV}/c$ (in the laboratory frame). The latter selection excludes a large fraction of low-momentum electrons from the analysis. The known positions of the target and the most probable intersection point of measured tracks define the position of the primary vertex.

$\boldsymbol{K^{0}_{S}}$ reconstruction. Unlike the charged particles, the neutral $\mathit{K}$ mesons do not leave a measurable track in the detectors. They are measured by reconstructing their oppositely charged decay products (daughter particles). The two-body decays of $K^{0}_{S}$ create characteristic $V$-shaped particle pairs originating at the decay vertex. This topology is called $V^{0}$. It is searched with a dedicated $V^{0}$-finder algorithm that looks for track pairs of particles with opposite charges. These track pairs are extrapolated backwards until their mutual distance of the closest approach is reached. If this distance is smaller than a given limit value, the track pair becomes a $V^{0}$ candidate with its origin at the decay vertex.

Two further cuts are placed on the track pairs. The first cut imposes a minimum value on the angle between the direction of the line joining the primary and decay vertices and the direction given by the vector sum of the momenta of the decay daughters. The second condition requires a minimum distance between the primary and decay vertex (a minimum length of the decaying particle). The corresponding cut requirements depend on $K^{0}_{S}$ rapidity and are listed in Methods’ "Extended data" subsection. Starting with the approximate decay point, a $V^{0}$-fitter program optimizes the decay point position and the momenta of the decay daughters. Assuming that the daughter particles are pions, it is straightforward to reconstruct the invariant mass of the decaying particle (the invariant mass is defined as $m_{inv}=\sqrt{(\sum{E_{i}})^{2}-(\sum{\mathbf{p}_{i}})^{2}}$, where $E_{i}$ are the energies of the decay products, $\mathbf{p}_{i}$ are their momenta, and $c\equiv 1$ is assumed).

The invariant mass distribution of $V^{0}$ candidates is populated by $K^{0}_{S}$ and $\Lambda$ decays, photon conversion in some detector material, and spurious particle crossings. A $K^{0}_{S}$ signal will appear as a peak on a slowly varying background. For a double-differential $K^{0}_{S}$ analysis, the momentum space was divided into seven rapidity bins ranging from $-1.5$ to 2 and nine transverse-momentum bins ranging from 0 to 2.7 $\mbox{Ge\kern-1.00006ptV}\!/\!c$. The raw number of $K^{0}_{S}$ in a given kinematic bin is obtained from fits of appropriate signal and background functions to the invariant mass distribution of the corresponding $V^{0}$ candidates. The fitted signal function is taken as a Lorentzian, and the background function is a third-order Chebychev polynomial. The integral of the signal function divided by the bin width is equal to the raw (uncorrected) number of the reconstructed $K^{0}_{S}$ in a given kinematic bin. Two typical invariant mass distributions with signal and background fits are shown in Methods’ "Extended data" subsection Fig. 5.

Corrections. To correct the results for losses due to detection and data processing inefficiencies, detailed Monte Carlo simulations were performed. These simulations comprised Ar+Sc collisions generated by the EPOS model [84], and particles propagated in the NA61/SHINE detector using the GEANT framework [85]. The charged particle tracks were reconstructed and analyzed using the same software as used for the experimental data. The branching ratio of $K^{0}_{S}$ decays was considered in the GEANT framework. The final output of the simulation consisted of reconstructed $K^{0}_{S}$ multiplicities. The ratio of the simulated and reconstructed numbers of $K^{0}_{S}$ was used as a correction factor in each $y$–$p_{T}$ bin.

Systematic uncertainties of the measured data points were estimated by comparing the results of the entire analysis (including Monte Carlo simulations and corrections) obtained with varying cut values. The reliability of the $V^{0}$ reconstruction and $K^{0}_{S}$ fitting procedures can be scrutinized by studying the $K^{0}_{S}$ lifetime. Methods’ "Extended data" subsection Fig. 6 shows the computed mean lifetime of $K^{0}_{S}$ in seven rapidity bins. Good agreement with the average value provided by the PDG [3] is observed.

Transverse momentum distributions. The distributions shown in Methods’ "Extended data" subsection Fig. 7 represent the final results of the $K^{0}_{S}$ analysis. The $K^{0}_{S}$ yields are shown as a function of transverse momentum in seven bins of rapidity. The data points are fitted with the function:

$$f(p_{T})=A\cdot p_{T}\cdot\exp\left(-\leavevmode\nobreak\ \frac{\sqrt{p_{T}^{2}+m_{0}^{2}}}{T}\right),$$ (4)

in which $A$ is a normalisation factor, $T$ is the inverse slope parameter, and $m_{0}$ is the $K^{0}_{S}$ mass taken from Ref. [3]. The formula assumes $c\equiv 1$ for simplicity. The fit functions are plotted as red curves, and the inverse slope parameters obtained from the fits are reported in the figure legends.

The transverse momentum distributions of charged and neutral $K$ mesons drawn in Fig. 2 (of the main text) are also fitted with the function defined by Eq. (4). The bottom panel of the figure presents the ratio of the two fitted curves, with its uncertainty band obtained by the propagation of the uncertainties of the fitted parameters.

Rapidity distribution. The final $K^{0}_{S}$ yields in each bin of rapidity were obtained as the integrals of the curves fitted to the respective transverse momentum spectra, Eq. (4), including extrapolations to unmeasured regions. A comparison to the alternative method of replacing integrals in the measured regions by sums of data points only brought a negligible contribution to the systematic uncertainty. Methods’ "Extended data" subsection Fig. 7 showed that extrapolations were needed only in the first and last rapidity bin. They amount to 88% and 6.2%, respectively. The large extrapolation in the first bin of rapidity increases the total uncertainty of the corresponding data point shown in Fig. 1 (of the main text). In this figure, the obtained rapidity distribution of the $K^{0}_{S}$ has been fitted with a function consisting of two Gaussians with centers displaced by a value of $\pm\Delta y$ with respect to $y=0$. These Gaussians have the same widths but may have different amplitudes. The resulting small asymmetry of the fitted rapidity distribution originates from a combined effect of the mass asymmetry of the colliding target and projectile nuclei ($A_{target}=45$ and $A_{projectile}=40$) and the selection of central collisions by the energy measured in the kinematic region of the projectile spectator remnants. The former favors backward and the latter forward rapidities. The yields of charged $K$ mesons at mid-rapidity listed in Table 1 (of the main text) were taken from Ref. [18]. They were determined in the interval $0.0<y<0.2$ as discussed therein. The yield of neutral $K^{0}_{S}$ mesons at mid-rapidity was determined at $y=0$ from the aforementioned fit. Its systematic uncertainty was estimated the same way as for the data points (see above), and its statistical uncertainty was obtained by propagation of the statistical uncertainties of the fit. Both statistical and systematic uncertainties of charged and neutral $K$ yields were propagated into the ratio $R_{K}$. The additional uncertainty of $R_{K}$ resulting from the difference in the mid-rapidity definition for charged and neutral mesons was estimated to be 0.5%, about 10% of the total systematic uncertainty.

B. World data

This section presents the yields of charged and neutral kaons measured by various experiments across different collision systems and energies and within specified centrality and rapidity regions. The results, presented in Table LABEL:expConditions, are sourced directly from the original experimental publications without any modifications to ensure consistency of the quantities reported. The exceptions are HADES $K^{+}$ and $K^{-}$ yields, where two sources of systematic uncertainties were reported [43]. In our analysis, they were added in quadrature, and the square root of such a sum is shown in Table LABEL:expConditions as the final systematic uncertainty (although in further calculations we used more precise values than 0.0014 and 0.000032 displayed in the table). In the NA49 experiment, the $K^{+}$ yield in Pb+Pb at $\sqrt{s_{NN}}$ $=$ 7.6 GeV [40] was reported with asymmetric systematic uncertainty; in this case, the upper limit was taken as $\sigma_{sys}$. For $K^{+}$ and $K^{-}$ yields in NA35 S+S collisions at $\sqrt{s_{NN}}$ $=$ 19.4 GeV only statistical uncertainties were reported in the form of numerical values [38]. We took the NA35 estimate of systematic uncertainty as $\thicksim$3% [38], and the resulting numerical values are presented in the table. Finally, for the STAR experiment at $\sqrt{s_{NN}}$ $=$ 130 GeV [32], two types of uncertainties were reported: uncorrelated errors (first) and correlated systematic errors (second); see Ref. [32] for details.

For all kaon yields reported with statistical and systematic uncertainties separately, we calculated the total uncertainties as $\sigma_{total}=\sqrt{\sigma_{stat}^{2}+\sigma_{sys}^{2}}$ (for STAR at $\sqrt{s_{NN}}$ $=$ 130 GeV $\sigma_{total}$ was taken as $\sqrt{\sigma_{uncorr}^{2}+\sigma_{corr}^{2}}$). Their rounded (to two significant digits) values are displayed in the third column of Table LABEL:expConditions however, more precise values were used when propagating them to $\sigma_{total}$ of $R_{K}$ presented in Table LABEL:Rkvalues.

The notation "Yield (4$\pi$)" refers to particle mean multiplicity in full phase space. The "Yield ($y\approx 0$)" corresponds to mid-rapidity production, in most cases expressed as rapidity density $\mathrm{d}n/\mathrm{d}y$ measured in the region specified in Table LABEL:expConditions as "${y}$ range" (for CERES results and NA61/SHINE $K^{0}_{S}$ mesons the fits at mid-rapidity were used). In some cases, different intervals were used for charged and neutral kaons. When calculating the charged-to-neutral kaon ratio, we used the originally published results.

The HADES data for Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 2.4 GeV [86, 87], the FOPI data for Al+Al collisions at $\sqrt{s_{NN}}$ $=$ 2.7 GeV [88, 89], and the NA49 data for Pb+Pb collisions at $\sqrt{s_{NN}}$ $=$ 17.3 GeV [41, 42] are excluded from this paper, as the charged and neutral kaons were measured in significantly different centrality intervals. Normalizing these results by the number of participants would introduce model dependence of the $R_{K}$ ratio. Moreover, we also omit kaon yields evaluated by the Authors of Ref. [90] based on rapidity spectra measured by AGS experiments in Si+Al/Si collisions at $\sqrt{s_{NN}}$ $=$ 5.4 GeV. The spectra of charged and neutral kaons were measured for different centralities [90], and the type of presented uncertainties is not clear. Finally, we also exclude NA35 kaon yields from S+Ag collisions at $\sqrt{s_{NN}}$ $=$ 19.4 GeV [90, 39]. The type of uncertainties for charged kaon yields [90] is not specified, and the charged and neutral kaons might have been measured for different centralities [90, 91, 39].

Table 2: The compilation of world data on charged and neutral kaon yields in nucleus-nucleus collisions. The yields labeled $y\approx 0$ and $4\pi$ correspond to the rapidity density at mid-rapidity and mean multiplicity in full phase space. The uncertainty fields are left empty in case they are not published. The systematic uncertainties labeled by ^(∗) were estimated for this analysis based on the information given in the original papers (see the text). Only results corresponding to the same centrality for charged and neutral kaons are compiled. The "${y}$ range" specifies the rapidity range used to obtain a given kaon yield.

NA61/SHINE experiment
Ar+Sc collisions at $\sqrt{s_{NN}}$ $=$ 11.9 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	3.732 $\pm$ 0.016 $\pm$ 0.148	0.15	0–10%	0.0 < $y$ < 0.2	[18]
$K^{-}$	2.029 $\pm$ 0.012 $\pm$ 0.069	0.070	0–10%	0.0 < $y$ < 0.2	[18]
$K_{S}^{0}$	2.433 $\pm$ 0.027 $\pm$ 0.102	0.11	0–10%	${y}$ $=$ 0	this analysis
HADES experiment
Ar+KCl collisions at $\sqrt{s_{NN}}$ $=$ 2.6 GeV
hadron	Yields ($4\pi$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	0.028 $\pm$ 0.002 $\pm$ 0.0014 (∗)	0.0024	0–35%	extrapolated to $4\pi$	[43]
$K^{-}$	0.00071 $\pm$ 0.00015 $\pm$ 0.000032 (∗)	0.00015	0–35%	extrapolated to $4\pi$	[43]
$K_{S}^{0}$	0.0115 $\pm$ 0.0005 $\pm$ 0.0009	0.0010	0–35%	extrapolated to $4\pi$	[44]
STAR (BES I) experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 7.7 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	20.8	1.7	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K^{-}$	7.7	0.6	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K_{S}^{0}$	12.67$\pm$ 0.12 $\pm$ 0.44	0.46	0–5%	$-$0.5 < $y$ < 0.5	[31]
STAR (BES I) experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 11.5 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	25.0	2.5	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K^{-}$	12.3	1.2	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K_{S}^{0}$	15.93 $\pm$ 0.12 $\pm$ 0.58	0.59	0–5%	$-$0.5 < $y$ < 0.5	[31]
STAR (BES I) experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 19.6 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	29.6	2.9	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K^{-}$	18.8	1.9	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K_{S}^{0}$	20.89 $\pm$ 0.08 $\pm$ 0.67	0.67	0–5%	$-$0.5 < $y$ < 0.5	[31]
STAR (BES I) experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 27 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	31.1	2.8	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K^{-}$	22.6	2.0	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K_{S}^{0}$	23.24 $\pm$ 0.09 $\pm$ 0.70	0.71	0–5%	$-$0.5 < $y$ < 0.5	[31]
STAR (BES I) experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 39 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	32.0	2.9	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K^{-}$	25.0	2.3	0–5%	$-$0.1 < $y$ < 0.1	[30]
$K_{S}^{0}$	24.9 $\pm$ 0.1 $\pm$ 1.7	1.7	0–5%	$-$0.5 < $y$ < 0.5	[31]
NA49 experiment
Pb+Pb collisions at $\sqrt{s_{NN}}$ $=$ 7.6 GeV
hadron	Yields ($4\pi$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	52.9 $\pm$ 0.9 $\pm$ 3.5 (∗)	3.6	0–7.2%	extrapolated to $4\pi$	[40]
$K^{-}$	16.0 $\pm$ 0.2 $\pm$ 0.4	0.45	0–7.2%	extrapolated to $4\pi$	[40]
$K_{S}^{0}$	29.3 $\pm$ 0.3 $\pm$ 2.9	2.9	0–7.2%	extrapolated to $4\pi$	[42]
NA49 experiment
Pb+Pb collisions at $\sqrt{s_{NN}}$ $=$ 8.7 GeV
hadron	Yields ($4\pi$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	59.1 $\pm$ 1.9 $\pm$ 3	3.6	0–7.2%	extrapolated to $4\pi$	[41]
$K^{-}$	19.2 $\pm$ 0.5 $\pm$ 1.0	1.1	0–7.2%	extrapolated to $4\pi$	[41]
$K_{S}^{0}$	34.2 $\pm$ 0.2 $\pm$ 3.4	3.4	0–7.2%	extrapolated to $4\pi$	[42]
CERES experiment
Pb+Au collisions at $\sqrt{s_{NN}}$ $=$ 17.3 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	31.8 $\pm$ 0.6 $\pm$ 2.5	2.6	0–7%	$y$ $=$ 0	[27]
$K^{-}$	19.3 $\pm$ 0.4 $\pm$ 2.0	2.0	0–7%	$y$ $=$ 0	[27]
$K_{S}^{0}$	21.2 $\pm$ 0.9 $\pm$ 1.7	1.9	0–7%	$y$ $=$ 0	[28, 29]
NA35 experiment
S+S collisions at $\sqrt{s_{NN}}$ $=$ 19.4 GeV
hadron	Yields ($4\pi$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	12.5 $\pm$ 0.4 $\pm$ 0.375 (∗)	0.55	0–2%	extrapolated to $4\pi$	[38]
$K^{-}$	6.9 $\pm$ 0.4 $\pm$ 0.207 (∗)	0.45	0–2%	extrapolated to $4\pi$	[38]
$K_{S}^{0}$	10.5	1.7	0–2%	extrapolated to $4\pi$	[39]
STAR experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 62.4 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	37.6	2.7	0–5%	$-$0.1 < $y$ < 0.1	[33]
$K^{-}$	32.4	2.3	0–5%	$-$0.1 < $y$ < 0.1	[33]
$K_{S}^{0}$	27.4 $\pm$ 0.6 $\pm$ 2.9	3.0	0–5%	$-$1 < $y$ < 1	[34]
STAR experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 130 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{uncorr}$ $\pm$ $\sigma_{corr}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	46.2 $\pm$ 0.6 $\pm$ 6.0	6.0	0–6%	${}^{-0.1}_{-0.5}$ <$y$ < ${}^{0.1}_{0.5}$	[32]
$K^{-}$	41.9 $\pm$ 0.6 $\pm$ 5.4	5.4	0–6%	${}^{-0.1}_{-0.5}$ < $y$ < ${}^{0.1}_{0.5}$	[32]
$K_{S}^{0}$	33.9 $\pm$ 1.1 $\pm$ 5.1	5.2	0–6%	$-$0.5 < $y$ < 0.5	[32]
STAR experiment
Au+Au collisions at $\sqrt{s_{NN}}$ $=$ 200 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	51.3	6.5	0–5%	$-$0.1 < $y$ < 0.1	[33]
$K^{-}$	49.5	6.2	0–5%	$-$0.1 < $y$ < 0.1	[33]
$K_{S}^{0}$	43.5	2.4	0–5%	$-$0.5 < $y$ < 0.5	[35]
ALICE experiment
Pb+Pb collisions at $\sqrt{s_{NN}}$ $=$ 2760 GeV
hadron	Yields ($y\approx 0$) $\pm$ $\sigma_{stat}$ $\pm$ $\sigma_{sys}$	$\sigma_{total}$	Centrality	${y}$ ranges	Ref.
$K^{+}$	109	9	0–5%	$-$0.5 < $y$ < 0.5	[36]
$K^{-}$	109	9	0–5%	$-$0.5 < $y$ < 0.5	[36]
$K_{S}^{0}$	110	10	0–5%	$-$0.5 < $y$ < 0.5	[37]

Table LABEL:Rkvalues presents the ratios of charged-to-neutral kaons from various experiments, with estimated statistical and total uncertainties where available. Taking into account the possibility of using our compilation in future analyses, the numerical values in Table LABEL:Rkvalues are presented with unusually high precision.

Table 3: Ratios of charged kaons to neutral kaons in different experiments.

Experiment	Collision system	$\sqrt{s_{NN}}$ (GeV)	$R_{K}$	$\sigma_{stat}$	$\sigma_{total}$
NA61/SHINE	Ar+Sc	11.9	1.1839	0.0138	0.0615
HADES	Ar+KCl	2.6	1.2483	0.1027	0.1545
STAR (BES I)	Au+Au	7.7	1.1247	-	0.0819
STAR (BES I)	Au+Au	11.5	1.1707	-	0.0973
STAR (BES I)	Au+Au	19.6	1.1584	-	0.0910
STAR (BES I)	Au+Au	27	1.1553	-	0.0819
STAR (BES I)	Au+Au	39	1.1446	-	0.1079
NA49	Pb+Pb	7.6	1.1758	0.0198	0.1325
NA49	Pb+Pb	8.7	1.1447	0.0295	0.1263
CERES	Pb+Au	17.3	1.2052	0.0539	0.1340
NA35	S+S	19.4	0.9238	-	0.1533
STAR	Au+Au	62.4	1.2774	-	0.1525
STAR	Au+Au	130	1.2994	-	0.2331
STAR	Au+Au	200	1.1586	-	0.1214
ALICE	Pb+Pb	2760	0.9909	-	0.1071

C. Models

Hadron Resonance Gas model. We use the Hadron Resonance Gas model implementation from Ref. [15] to quantify the isospin-breaking effects and their interplay. HRG includes all hadrons and resonances with confirmed status in the PDG tables [92]. The PDG-listed masses, charges, lifetimes, and decay modes are used. Thus, HRG includes the isospin-symmetry violation due to masses and branching ratios of hadrons and resonances.

In HRG calculations, the exact net strangeness conservation is enforced, i.e., the calculations are done within the strangeness canonical ensemble (SCE) [93, 94]. The model parameters are baryo-chemical potential, $\mu_{B}$, temperature, $T$, volume of the system $V$, and the strangeness under-saturation parameter $\gamma_{S}$ (see Ref. [95]). We adopt the simple parametrization of $\mu_{B}$ and $T$ as a function of collision energy introduced in Ref. [96]. We have checked that $R_{K}$ is weakly sensitive to the strangeness suppression effect introduced by having the parameter $\gamma_{S}<1$. Therefore, calculations are done for $\gamma_{S}=1$ for simplicity. The energy-dependent Breit-Wigner spectra [97] model the resonance widths.

Figure 3 (black line) shows the HRG predictions for $R_{K}$ as a function of collision energy $\sqrt{s_{NN}}$ for $Q/B=0.4$. At low collision energies, $R_{K}<1$ due to the enhancement of neutral kaon production caused by a larger number of neutrons than protons, $Q/B<1/2$. On the other hand, at high collision energies, HRG predicts $R_{K}\approx 1.018$ due to the mass difference between charged and neutral kaons produced directly. Finally, $R_{K}\approx 1.032$ when kaon production via resonance decays is included. This is mostly due to $\phi$ decays, which strongly prefer the decay into charged kaons over the one into neutral kaons.

We have checked that the system electric to baryon charge ratio $Q/B$ significantly affects $R_{K}$ up to $\sqrt{s_{NN}}\approx 10$ GeV. At higher energies, pions dominate, and total electric and baryon charges are significantly larger than the corresponding net charges. Thus, $R_{K}$ becomes increasingly less sensitive to $Q/B$ with an increase in $\sqrt{s_{NN}}$.

We have checked that the strangeness grand-canonical ensemble and other popular parametrizations [98, 47] of the model parameters as a function of collision energy lead to quantitatively similar results for $R_{K}$, for $\sqrt{s_{NN}}\gtrsim 4$ GeV. The uncertainties of $R_{K}$ estimated by changing the parametrization of model parameters [96] are less than $1$% for $\sqrt{s_{NN}}>3$ GeV. The effect of including light nuclei in the particle list is negligible.

A dedicated discussion for resonances $a_{0}(980)$ and $f_{0}(980)$ is needed. For $a_{0}(980)$ three strong decay channels are reported as seen in PDG [3]: $\pi\eta$, $K\overline{K}$, and $\pi\eta^{\prime}$, where the latter is phase-space suppressed because $m_{\pi}+m_{\eta^{\prime}}=1.093\text{ GeV}>m_{a_{0}}=0.98\pm 0.02\text{ GeV}$. The PDG average ratio of $a_{0}(980)$ decay widths $\Gamma_{K\overline{K}}/\Gamma_{\pi\eta}=0.172\pm 0.019$ [3] implies that the $K\overline{K}$ branching ratio amounts to $\left(15\pm 1\right)\%$:

$$BR(K\overline{K})\approx\frac{\Gamma_{K\overline{K}}}{\Gamma_{K\overline{K}}+\Gamma_{\pi\eta}+\Gamma_{\pi\eta^{\prime}}}\approx\frac{\Gamma_{K\overline{K}}}{\Gamma_{K\overline{K}}+\Gamma_{\pi\eta}}=\frac{1}{1+\frac{\Gamma_{\pi\eta}}{\Gamma_{K\overline{K}}}}\approx 0.15\pm 0.01\text{ .}$$ (5)

For $f_{0}(980)$ the strong decay channels $\pi\pi$ and $K\overline{K}$ have been seen [3]. A list of measurements for the $\pi\pi$ branching ratios is reported, each of them larger than $50\%$, but no PDG average is provided. Building an average for the $K\overline{K}$ branching ratio leads to $\left(19\pm 2\right)\%$, but to enhance reliability we vary the $K\overline{K}$ range between 10–40%. As a starting point, the HRG approach used in this work assumes $BR(K^{+}K^{-})$ $=$ $BR(K^{0}\overline{K}^{\,0})$ for both resonances. Taking into account the difference of the charged and neutral kaon masses within Flatté-like distributions [99, 100, 101], the ratio of charged and neutral decay rates is about $1.1$ and remains smaller than $1.2$ when the distribution parameters are varied within reasonable ranges (see, e.g. the compilation in Ref. [49]). We therefore recalculated the HRG predictions assuming $20\%$ more charged than neutral kaons produced by decays of $a_{0}^{0}(980)$ and $f_{0}(980)$. The ratio $R_{K}$ increases by less than $0.5\%$. Under the extreme assumption $BR(K^{+}K^{-})$ $=$ $BR(K\overline{K})$ (no decays to neutral kaons) for both resonances, $a_{0}(980)$ and $f_{0}(980)$, and taking the upper limit for the $K\overline{K}$ branching ratio of $f_{0}(980)$ equal to $40\%$, $R_{K}$ increases by up 3% at high collision energies.

UrQMD model. The UrQMD transport model [45, 16, 46] describes $A$+$A$ collisions by explicitly propagating hadrons in phase space. During the propagation, rescattering among hadrons takes place. The particle production in this model happens via resonance decay or string excitation and fragmentation following the LUND model [102].

The gray squares in Fig. 3 indicate the UrQMD model predictions. Here, we have considered central Au+Au collisions ($A$ $=$ 197, $Z$ $=$ 79, $Q/B\approx 0.4$). The predictions are shown within the $\sqrt{s_{NN}}$ range of 2.4 to 20 GeV. At each energy, $10^{4}$ events are used for the analysis.

One sees that for $\sqrt{s_{NN}}\lesssim 7$ GeV, the predictions of UrQMD and HRG are similar. At higher energies, the $R_{K}$ ratio in HRG is systematically higher than the one predicted by UrQMD. This is likely caused by UrQMD assuming $\phi$-meson decays to be exactly isospin symmetric instead of taking the branching ratios from PDG. This is the reason for showing the UrQMD predictions only up to 20 GeV.

D. Extended data

In this part, we present values of track pair cuts used in the analysis (Table 4), examples of fitted invariant mass distributions (Fig. 5), mean lifetime (Fig. 6) and transverse momentum distributions (Fig. 7) of $K^{0}_{S}$ in rapidity bins.

rapidity bin		$(-1.5,-1)$	$(-1,-0.5)$	$(-0.5,0)$	$(0,0.5)$	$(0.5,1)$	$(1,1.5)$	$(1.5,2)$
cut value	cosine of angle	$>$0.999	$>$0.9995	$>$0.9995	$>$0.9995	$>$0.9995	$>$0.9999	$>$0.9999
	distance	$>$5 cm	$>$5 cm	$>$7.5 cm	$>$12.5 cm	$>$12.5 cm	$>$15 cm	$>$12.5 cm

Table 4: Track pair cuts. Values of the cuts on (top row) the cosine of the angle between the line joining the primary and decay vertex and the direction of the vector sum of decay daughter momenta, and (bottom row) the distance between the primary and decay vertex.

Figure 5: Examples of fitted invariant mass distributions. Two studied bins in rapidity $y$ and transverse momentum $p_{T}$ of the $K^{0}_{S}$ are presented, left: $y\in(-1.0,-0.5)$, $p_{T}\in(1.2,1.5)$ $\mbox{Ge\kern-0.90005ptV}\!/\!c$, right: $y\in(0.5,1.0)$, $p_{T}\in(1.2,1.5)$ $\mbox{Ge\kern-0.90005ptV}\!/\!c$. Only statistical uncertainties are presented. The bottom panels show the difference between the experimental data and the fitted (Signal+Background) distribution, divided by the experimental uncertainty.

Figure 6: Mean lifetime of K${}^{\mathbf{0}}_{S}$ mesons as a function of rapidity. The values obtained by NA61/SHINE are divided by the PDG value [3]. Statistical uncertainties are shown by vertical bars and systematic ones by shaded boxes.

Figure 7: K${}^{\mathbf{0}}_{S}$ transverse momentum spectra in rapidity bins. Statistical uncertainties are shown by vertical bars and systematic ones by shaded boxes. Red curves represent fits of the data with the function defined in Eq. (4). The inverse slope parameters ($T$), with their statistical uncertainties resulting from the fits, are also displayed inside the panels.

Data availability
All data shown in plots are publicly available on the HEPdata repository (https://hepdata.net).

Code availability
The authors can provide the NA61/SHINE source code used upon reasonable request.