Constraints on long-lived light scalars with flavor-changing couplings and the KOTO anomaly

Recently, the KOTO experiment at J-PARC has observed three anomalous events in the flavor-changing rare decayKL → π0νν̄, which indicates that the corresponding branching ratio is almost two orders of magnitude larger than the Standard Model (SM) prediction. Taking this intriguing result at face value, we explore model implications of its viable explanation by a long-lived light SM-singlet scalar (S) emission, i.e. KL → π0S, with S decaying outside the KOTO detector. We derive constraints on the parameter space of such a light scalar in the context of three simple models: (i) a real singlet scalar extension of the SM; (ii) a B−L extension where neutrino masses arise via type-I seesaw mechanism from B−L breaking; and (iii) a TeV-scale left-right symmetric model. The flavor-changing couplings needed to explain the KOTO excess in models (i) and (ii) originate from treelevel mixing of the scalar with SM Higgs field (h), and in model (iii), from the mixing of S and h with the neutral component of the heavy bidoublet Higgs field. After taking into account the stringent constraints from high-precision searches for flavor-changing charged and neutral kaon decays at NA62, E949, KOTO and CHARM experiments, as well as the astrophysical and cosmological constraints on a light scalar, such as those from supernova energy loss, big bang nucleosynthesis and relativistic degrees of freedom, we find that the light scalar interpretation of the KOTO excess is excluded in models (i) and (ii) by the supernova limit. However, there is still a narrow parameter space available in the TeVscale left-right symmetric model to explain the KOTO excess, which can be tested in future NA62 and DUNE experiments.


Introduction
In the Standard Model (SM) of particle physics, flavor-changing neutral currents (FC-NCs) are absent at tree-level and are predicted to be small at loop level, suppressed by the Glashow-Iliopoulos-Maiani (GIM) mechanism and the small off-diagonal Cabibbo-Kobayashi-Maskawa (CKM) matrix elements in the quark sector (or by the tiny neutrino masses in the lepton sector, if we allow nonzero neutrino masses to be part of the 'new' SM) [1]. Any observation of FCNC above the SM prediction would therefore be a clear signature of beyond the SM (BSM) physics. Very recently, the KOTO experiment at J-PARC has observed four candidate events in the signal region of one such rare flavor-changing decay K L → π 0 νν [2]. While one of the events is suspected to have originated from SM activity upstream from the detector and can be vetoed away, the remaining three events cannot be explained by currently known backgrounds, with the SM expectation of only 0.05 ± 0.02 events. This corresponds to a decay branching ratio (BR) of [2] BR(K L → π 0 νν) KOTO19 = 2.1 +2.0(+4.1) This result is consistent with their previously reported 90% CL upper bound of [3] BR(K L → π 0 νν) KOTO18 < 3.0 × 10 −9 . (1. 2) The central value in Eq. (1.1) is almost two orders of magnitude larger than the SM prediction of [4] BR(K L → π 0 νν) SM = 3.4 ± 0.6 × 10 −11 . (1. 3) Needless to say, more experimental information on the source of these intriguing events, as well as a careful reevaluation of background estimations, is needed to confirm whether the signal is indeed due to some BSM physics. But given the far-reaching consequences, we take the KOTO result (1.1) at face value and explore possible implications for some simple BSM scenarios that can be independently tested in other ongoing or future experiments. At the phenomenological level, the KOTO signal can be interpreted as the emission of a new light, long-lived scalar particle S in the two-body kaon decay K L → π 0 S, which subsequently decays outside the KOTO detector, thus mimicking the invisible νν final states in K L → π 0 νν [5][6][7]. 1 In this paper, we consider possible ultraviolet (UV)-complete model frameworks for such a new light scalar particle S with m S < m K − m π 0 and with a flavor-changing effective coupling of the form KπS so that it can be emitted in kaon decay. Being light, there are stringent constraints on this particle from laboratory searches for FCNCs in the K, D and B meson decays. In particular, the scalar emission through the effective KπS coupling contributes to both neutral and charged kaon decays, i.e. K L → π 0 S and K + → π + S, whose branching ratios (for an invisible S) are correlated by the Grossman-Nir bound [9] BR(K L → π 0 νν) ≤ 4.3 BR(K + → π + νν) .
( 1.4) No excess has been reported in the charged kaon decay mode K + → π + νν, whose branching ratio is currently constrained by the NA62 experiment [10] (and also by E949 experiment [11]) to be BR(K + → π + νν) NA62 < 2.44 × 10 −10 (1.5) at 95% CL, which is consistent with the SM prediction of [4] BR(K + → π + νν) SM = (8.4 ± 1.0) × 10 −11 . (1.6) There also exist stringent constraints on a light, long-lived scalar decaying into charged lepton or photon pairs from the searches for + − and γγ in rare kaon decays at proton beam-dump experiments, such as CHARM [12]. In addition, a light S particle will be constrained by astrophysical and cosmological observations, such as those from supernova energy loss, and effective relativistic degrees of freedom (∆N eff ) and/or additional energy injection at the big bang nucleosynthesis (BBN) epoch. To see whether the KOTO excess is consistent with all these constraints, it is convenient to work within specific models so that the new scalar interactions with SM particles have a definite profile. Due to the suppressed nature of these interactions, only models where the particle S has no tree-level direct coupling to SM quarks need to be considered. We find the following three BSM scenarios which fall into this category: (i) Scalar singlet model: Here the FCNC couplings of S arise from its mixing with the SM Higgs field h, which has loop-induced FCNC couplings with SM quarks [13][14][15]. The new scalar S could be long-lived if the S − h mixing angle, θ is suitably 1 Other interpretations in terms of either a heavy mediator or a new light particle produced at fixed target and decaying off-axis to two photons (e.g. an axion-like particle) have also been discussed [6]. Similarly, Ref. [8] has considered the possibility of KL → π 0 QQ, where Q is a dark fermion of the dark sector.
small, while avoiding all existing laboratory constraints [16][17][18][19][20][21]. This is a simple, two-parameter model with only m S and θ as the unknown parameters. We call this the SM+S model.
(ii) U (1) B−L model: A class of UV-complete models where such a light scalar without tree-level coupling to SM fermions emerges naturally is based on the gauge group [22][23][24][25]. In this case, three right-handed neutrinos (RHNs) N a (1, 1/2, −1) (with a = 1, 2, 3) are introduced for the purpose of anomaly cancellation. The light scalar S can be identified as the real part of a complex (B − L)-charged scalar ∆(1, −1, 2) that breaks the U (1) I 3R × U (1) B−L symmetry to U (1) Y of the SM and gives mass to the RHNs to implement the type-I seesaw mechanism [26][27][28][29][30]. The long-lived property and the FCNC constraints of this model on S have already been studied in great detail in Ref. [31] (where S was denoted by H 3 ), which will be relied upon here. This model has some new, suppressed decay modes such as S → N N, Z Z (where Z is the massive gauge boson associated with the U (1) B−L breaking) which are absent in the SM+S model. However, as long as all the three RHNs and the Z boson are much heavier than the light scalar S, they will not have any effect on the lifetime of S, and the KOTO phenomenology of light S in the U (1) B−L extension will be the same as in the SM+S model. In what follows, we assume this to be the case and therefore do not separately discuss the U (1) B−L scenario for the KOTO explanation, except for the complementary collider signatures, which are different in the U (1) B−L case due to the additional gauge-portal production.
(iii) Left-right symmetric model: The last class of models studied here is the left-right symmetric model (LRSM) based on the gauge group SU (2) L ×SU (2) R ×U (1) B−L [32][33][34]. Here the light scalar S (denoted by H 3 in Refs. [31,35]) can be identified as the real part of the neutral component of the (B − L)-charged, SU (2) R -triplet field ∆ R (1, 3, 2), which can be light and does not couple directly to SM quark fields prior to symmetry breaking [35][36][37]. It is therefore similar to the SM+S model in many respects and can play a role in resolving the KOTO anomaly. The field ∆ R is responsible for the SU (2) R × U (1) B−L symmetry breaking and the model, like the U (1) B−L model above, has the extra motivation of being connected to neutrino mass generation via type-I seesaw [26][27][28][29][30]. In contrast with the previous two models, the FCNC couplings of S in this case arise at tree-level, due to its mixing with the heavy scalar H 1 from the bidoublet Φ (and the SM Higgs). Another special feature of the light scalar S in the LRSM is that even for small mixing angles , it can still decay into two photons through the W R loop and the heavy charged scalar loops. This makes the FCNC limits, as well as the supernova and BBN limits, on light S in the LRSM very different from the other two models discussed above.
As we show below, it turns out the parameter space preferred by the KOTO excess in the SM+S model, as well as in the U (1) B−L model, is excluded by the constraints from supernova energy loss criteria, whereas the LRSM could still explain the KOTO excess, while being consistent with all laboratory and cosmological/astrophysical constraints. Our results can be tested in the future high-precision intensity frontier experiments, such as NA62 and DUNE. The rest of the paper is organized as follows: In Sec. 2, we discuss the simplest real scalar extension of the SM in light of the KOTO excess vis-á-vis other laboratory and astrophysical/cosmological constraints. Most of this discussion is also applicable to the U (1) B−L case. In Sec. 3, we repeat the same exercise for the LRSM. Our conclusions are given in Sec. 4.

Singlet model
The singlet scalar extension of the SM is one of the simplest and well-motivated BSM scenarios [16]. The most general renormalizable scalar potential of the SM Higgs doublet H and a real singlet scalar S can be written as with µ 2 1,2 > 0 being the mass parameters and λ 1,2,3 being the quartic couplings. We impose a Z 2 -symmetry under which S → −S to prevent the S 3 trilinear term. After spontaneous symmetry breaking, the H and S fields obtain non-vanishing VEVs, with H = (0, v EW ) T with v EW 174 GeV the electroweak (EW) VEV and S = v S . The h − S mixing (where h is the physical SM Higgs field, obtained by expanding the H-field around its VEV, i.e. H = (0, v EW + h) T ) is determined by the quartic coupling λ 2 . In the small mixing limit, the mass of the real component of S is m 2 S 4λ 3 v 2 S to the leading order. For sufficiently small λ 3 and v S , the scalar S could be very light, even down to a few MeV scale. 2 In the U (1) I 3R × U (1) B−L extension discussed in Sec. 1, the S-field can be identified as the real part of a (B − L)-charged scalar, whose VEV breaks the U (1) I 3R × U (1) B−L gauge symmetry down to U (1) Y [31]. We assume the RHNs in this case are all heavier than S, so that the decays of S are identical to those of the SM+S case, being governed only by two parameters, namely, the scalar mass m S and the h − S mixing angle sin θ.

Fitting the KOTO anomaly
The couplings of S to the SM fermions arise from its mixing with the SM Higgs h and are thus flavor-conserving at the tree level. However, FCNCs are generated at one-loop level through Penguin diagrams involving the W −top loop and CKM quark mixings. The effective Lagrangian relevant for FCNC kaon decay is given by [17] where y sd is the effective loop-level coupling in the SM, G F is the Fermi constant, m s, t the strange and top quark masses, and V td, ts the CKM matrix elements. As a result of the CP phase in the CKM matrix, the coupling y sd is complex. If kinematically allowed, this will induce the flavor-changing decays K ± → π ± S and K L → π 0 S, with the partial widths with the kinematic function (2.5) Note that the partial decay widths in Eqs. (2.3) and (2.4) are almost identical, except for the crucial difference that the decay K L → π 0 S depends only on the real part of the coupling y sd . The same h − S mixing is also responsible for S decays into the SM quarks (u, d, s) and charged leptons at tree-level, and gluons and photons at one-loop level, just like the SM Higgs boson decays. In the generic singlet model all the decay modes of S are universally proportional to the mixing angle sin θ, and therefore, the branching ratios depend only on the S mass but not on sin θ. As detailed in Ref. [31], if S is light, say below the GeV-scale, it tends to be long-lived for a wide range of sin θ. If its average decay length is larger than the KOTO detector size, the process of interest will be with S decaying into anything outside the detector. This has the same final state as the decay K L → π 0 νν, i.e. two photons from π 0 → γγ and significant missing energy. In this case, the effective branching ratio is given by 3 where BR(K L → π 0 S) = Γ(K L → π 0 S)/Γ total K L , L = 3 m for the KOTO detector, and b = E S /m S the Lorentz boost factor with energy E S 1.5 GeV. For the total decay width of K L , we use Γ total Using Eq. (2.7), we calculate the preferred region in the (m S , sin θ) parameter space that explains the KOTO excess given by Eq. (1.1) at 95% CL. Our result is shown by the green shaded region in Fig. 1, with the green dashed line corresponding to the KOTO central value in Eq. (1.1). The region with m S > 180 MeV is not included in this fit, because it does not overlap with the KOTO signal region [6]. . For comparison, we show the exclusion regions from a previous KOTO search for K L → π 0 νν (brown shaded) [3], NA62 search for K + → π + νν (blue shaded) [10], E949 search for K + → π + X (magenta shaded) [11], and beam-dump experiment at CHARM (orange shaded) [12]; cf. Sec. 2.2. The pink shaded region is excluded by the supernova constraints (cf. Sec. 2.3).

Laboratory constraints
As shown by the Grossman-Nir bound [cf. Eq. (1.4)], the FCNC decays of charged and neutral kaons are correlated. This relation has to do with general isospin symmetry arguments, which relate the decay amplitudes of K ± to those of K 0 and K 0 , and holds even for the 2-body decays K L → π 0 S and K + → π + S [38]. For our singlet scalar case, this can be explicitly seen from Eqs. (2.3) and (2.4). As a result, the current and future high-precision measurements of the charged and neutral kaon rare decays can be used to set limits on the scalar mass m S and mixing angle sin θ, as discussed below.
In the charged kaon sector, the most stringent limits come from searches of K + → π + νν at NA62 [10] and of K + → π + X (with X being a long-lived particle) at E949 [11]. The NA62 experiment has put a 95% CL upper limit on BR(K + → π + νν) [cf. Eq. (1.5)] which can be translated into an exclusion region in the (m S , sin θ) plane, as shown by the blue shaded region in Fig. 1. Here we have constructed an effective BR, similar to Eq. (2.7), replacing neutral mesons by charged mesons and modifying the experimental parameters to L = 150 m and E S ∼ 37 GeV for NA62. Again we have neglected the O(1) kinematical difference between the 3-body SM decay and the 2-body decay in our scalar case. We see from Fig. 1 that there is a gap in the NA62 excluded region around the pion mass. This is because of the fact that if the scalar mass m S is close to π 0 mass, we will have a large pion background from the SM process K + → π + π 0 with π 0 → νν. With the current limit of BR(π 0 → νν) < 2.7 × 10 −7 [1] and the SM K + → π + π 0 branching ratio of 20.6% [1], the NA62 limit on K + → π + S turns out to be two orders of magnitude weaker in this region, as shown by the gap in Fig. 1.
The E949 experiment has reported 95% CL bounds on BR(K + → π + X), where X is a long-lived particle, as a function of the X mass [11]. The most stringent limit on the branching ratio BR(K + → π + X) could reach up to 5.4 × 10 −11 if the new particle X is stable. Using the same procedure as above for NA62, we evaluate the effective branching ratio following Eq. (2.7), with decay length L = 4 m and energy E S 710 MeV for E949. The corresponding exclusion region is shown by the magenta shaded region in Fig. 1, which is up to a factor of few stronger than the NA62 exclusion region in the low-mass range, but is weaker in the high-mass range and not applicable for m S > 2m π 0 because in this case the S tends to decay quickly, compared to the E949 detector size of 4 m. Like the NA62 limit, there is also a gap for the E949 constraint when m S ∼ m π .
Similarly, a previous KOTO search has reported 90% CL upper limits on the 2-body decay BR(K L → π 0 X), where X is an invisible boson, as a function of the X mass [3]. We can directly use this bound for our scalar case, and following Eq. (2.7) with L = 3 m and E S 1.5 GeV for KOTO, translate it into an exclusion region in the (m S , sin θ) plane, as shown in Fig. 1 by the brown shaded region. Note that there is no gap in the KOTO limit for m S ∼ m π , because the 2-body decay of K 0 L → π 0 π 0 is CP -violating and CKM-suppressed in the SM, with a branching ratio of 8.6 × 10 −4 [1].
Further limits on the light scalar can be derived from the e + e − , µ + µ − and γγ decay products of S produced in neutral and charged kaon decays at proton beam-dump experiment such as CHARM [12]. The production cross section of S at CHARM is given by [18,39,40] where σ pp is the proton-proton cross section, M pp = 11 the average hadron multiplicity, and χ s = 1/7 is the fraction of strange pair-production rate. In Eq. (2.8), the factor D K K /b K cτ K (with K standing for both K ± and K L ) takes into account the re-absorption of Kaons before decaying [7], with K = 15.3 cm the absorption length, b K = E K /m K the Lorentz boost factor with E K 25 GeV, and τ K the total Kaon decay width. It turns out that the re-absorption factors are respectively 8.1 × 10 −4 and 2.0 × 10 −4 for K ± and K L . Normalized to the neutral pion yield σ π 0 σ pp M pp /3, we can predict the total number of S particles produced: N S 2.9 × 10 17 σ S /σ π 0 . Then the number of events collected by the detector would be where L = 480 m is the CHARM beam dump baseline, ∆L = 35 m is the detector fiducial length, and b = E S /m S is the boost factor with E S E K /2 [12]. Due to the huge number of events N S , the mixing angle sin θ is expected to be severely constrained, and the most stringent limits are from the ee and µµ channels, as the γγ channel is comparatively suppressed by the loop factor. Given that no signal event was found at CHARM, an upper limit of N event < 2.3 at the 90% CL on the contribution from BSM physics was set. We use this limit to derive the corresponding exclusion region in the (m S , sin θ) plane, as shown by the orange shaded region in Fig. 1. For lighter S, the boost factor b becomes larger, and fewer S decays happen inside the detector, thereby weakening the constraints. As can be seen from Fig. 1, even if all the laboratory constraints are taken into consideration, the KOTO excess is still allowed in the singlet model. However, as we will see in Sec. 2.3, it is the supernova limit that deals a killer blow to this case. For the sake of completeness, we list here also other limits from the high-precision quark flavor data [31] that are either not applicable or weaker than those shown in Fig. 1.
• The lifetimes of K ± and K 0 are both precisely measured up to the level of 10 −3 , although the absolute theoretical values are subject to a large uncertainty of strange quark mass, up to the order of 10% [1]. Therefore, for sufficiently large mixing angles the contribution of K → πS to the total kaon decay widths will suppress the current uncertainties. However these limits from the total decay widths are comparatively much weaker than those from the rare decays discussed above.
• The limits from the beam-dump experiment NuCal apply only for a light scalar with mass m S 80 MeV [47], and they are not relevant to the KOTO anomaly here.
• In both the singlet scalar and U (1) B−L models, the loop-level FCNC structure is fixed by the SM quark masses and the CKM matrix [see Eq. (2. 2)], thus we can also use the flavor-changing data in the B meson sector to set limits on the S mass and mixing angle sin θ. However, as the B meson is much heavier than the K meson, the production rate of B mesons is much smaller, and as a result the limits from the rare decays of B → Kνν at BaBar and Belle are comparatively much weaker than those from the K meson decays, being respectively 3.2 × 10 −5 [48] and 1.6 × 10 −5 [49] (see Fig. 19 of Ref. [31]). The future prospects at Belle II [50] are also not comparable to those from the Kaon sector. Furthermore, the searches in the visible channels B → K + − ( = e, µ) at BaBar [51], Belle [52] and LHCb [53] are not applicable to the long-lived S case here.
• As detailed in Refs. [31,35], there are also some limits from the measurements of neutral K and B meson oscillations [1], from BR(B s → µ + µ − ) by LHCb [54], BR(B d → γγ) by BaBar [55] and BR(B s → γγ) by Belle [56]. However, these limits are much weaker and are not relevant to the KOTO anomaly.
If the light scalar S is long-lived, it can also be searched for at the LHC [57] and/or the dedicated long-lived particle (LLP) detectors such as MATHUSLA [58]. At the highenergy colliders, S can be produced from the loop-level gluon fusion process gg → gS via mixing with the SM Higgs, and the cross section can go up to (25 pb) × sin 2 θ [31]. The LLP searches at the HL-LHC and FCC-hh could probe a large parameter space of m S and sin θ (see Fig. 20 in Ref. [31]); however, they do not cover the KOTO parameter space of interest in Fig. 1, and hence, are not shown.

Astrophysical and cosmological constraints
A sufficiently light S can be produced in significant amount in the supernova core via the nuclear bremsstrahlung process with N = p, n collectively standing for protons and neutrons. This process is induced by the mixing of S with the SM Higgs field and the effective couplings of the SM Higgs to nucleons. Through the couplings to quarks inside nucleons, the effective couplings of the SM Higgs to nucleons are of order ∼ 10 −3 [59]. Following the calculations in Ref. [60], we estimate the total production rate of S in the supernova core. After being produced, S might decay inside or outside the supernova core, depending on its mass m S and the mixing angle sin θ. If S decays outside the core, then the mass and mixing angle will be constrained by the energy loss due to S emission and must be less than that inferred from SN1987A, i.e. about 3 × 10 53 erg [61]. Setting the core size to be R c = 10 km, we find that the pink shaded region in Fig. 1 is excluded from the supernova energy loss criteria. In the small m S limit, the supernova constraint could go down to sin θ 10 −7 . For larger mixing angles, the emitted S particles from the core get trapped inside the supernova. That happens when the mean path length of S is less than the supernova size. Note also that when m S 2m π , the hadronic decay modes of S will be kinematically open and S will decay much faster. Thus the supernova limit beyond m S 2m π shrinks very quickly. The limits nonetheless extend to m S values much larger than the supernova core temperature T SN 30 MeV, because of the super-high nuclear densities which make the bremsstrahlung process (2.10) very efficient. It is clear from Fig. 1 that the supernova limit excludes all the parameter space favored by the KOTO anomaly. In other words, a long-lived scalar in these two simple singlet scalar models (SM+S and the U (1) B−L model) can not accommodate the KOTO anomaly. If the mixing angle sin θ is very small, say 10 −5 , the lifetime of a light S might be longer than one second, which would potentially affect BBN in the early universe. However such small mixing angles are not relevant to the KOTO anomaly here, thus the BBN limit is not shown in Fig. 1. A detailed analysis of the BBN constraint in the context of the singlet scalar model can be found e.g. in Ref. [62].
A light S might also contribute to the relativistic degrees of freedom N eff in the early universe, thus getting constrained by the current precision Planck data [63]. However, if the mixing angle sin θ is too small, S cannot be kept in equilibrium with the SM photon.
In particular, if the mixing angle sin θ 0.01 and the scalar mass m S 1 MeV, the decay rate Γ(S → γγ) exp[−m S /T * ] will be Boltzmann suppressed and is significantly smaller than the Hubble expansion rate H 10T 2 * /M Pl , with T * ∼ T BBN ∼ MeV and M Pl the Planck mass. As the ∆N eff limit is very weak and not relevant to the KOTO anomaly, we do not show it in Fig. 1.

Left-right symmetric model
In the minimal version of the LRSM based on the gauge group SU (2) L ×SU (2) R ×U (1) B−L , there is one bidoublet Φ(2, 2, 0) and one right-handed triplet ∆ R (1, 3, 2) in the scalar sector: The SU (2) R × U (1) B−L symmetry is broken down to the SM U (1) Y gauge group, once the triplet develops a non-vanishing VEV ∆ 0 , is responsible for breaking the SM gauge group SU (2) L × U (1) Y down to U (1) em and for the generation of SM quark and charged lepton masses as well as the Dirac mass matrix for the type-I seesaw.
In the bidoublet sector, the SM Higgs h is predominantly from the real component of the neutral scalar φ 0 1 . There is a heavy CP-even scalar H 1 from the real component of φ 0 2 , which couples to the SM quarks through the couplings with q L,R = (u, d) T L,R the left-and right-handed quark doublets, Φ = iσ 2 Φ * (with σ 2 being the second Pauli matrix), and h q andh q the quark coupling matrices. After symmetry breaking, the tree-level couplings of H 1 to the SM quarks are flavor-changing, which are governed by the quark masses and the left-and right-handed quark mixing matrices V L, R in the form of with ξ = κ /κ the VEV ratio in the bidoublet sector, i, j the quark generation indices, and Y U,D diagonal Yukawa coupling matrices for the SM up-and down-type quarks. The tree-level FCNC couplings of H 1 contribute significantly to the neutral K and B meson oscillations, and thus H 1 is required to be superheavy, roughly above 15 TeV [36,64,65]. The CP-even neutral component S from the triplet ∆ R couples predominantly to the BSM scalars, heavy W R and Z R gauge bosons and the heavy RHNs in the LRSM, and all the couplings of S to the SM particles are from its mixings with the SM Higgs h and heavy H 1 [36,37]. Therefore in some region of the parameter space, even if the radiative corrections to S mass are taken into consideration, S can be very light, e.g. in the sub-GeV-scale [31,35]. Thus it might be a good candidate to explain the KOTO anomaly.

Fitting the KOTO anomaly
In the LRSM, the FCNC couplings of S are from mixing with the SM Higgs h and the heavy scalar H 1 from the bidoublet. Denoting these mixing angles respectively by sin θ 1, 2 , the FCNC couplings of S to s and d quarks will be proportional to the factor of (ξ sin θ 1 + sin Note that as a result of the CP phase in the V L, R matrices, this coupling is complex. The partial widths for the charged and neutral K meson decays are given by [40,66,67] with the kinematic function β 2 (M, m 1 , m 2 ) defined in Eq. (2.5).
The mixing angles of S to h and H 1 are strongly constrained by the low-energy highprecision flavor data, depending on the S mass in the LRSM [31]. At the one-loop level, S can decay into two photons, i.e. S → γγ, which is induced by the W R boson and the singly-and doubly-charged scalars [31,35]: where we have neglected the contributions from the SM fermion loops which are all highly suppressed by the mixing angles sin θ 1,2 , and take the limit of light S (compared to the BSM particles in the loop). Note that the partial width does not depend on the gauge coupling g R , as the dependence of W R couplings and W R mass on g R are canceled out. Thus, the partial width of S to diphoton is effectively suppressed only by the v R scale, independent of the mixing angles θ 1,2 .
As detailed in Refs. [31,35], if S is below the GeV-scale, it tends to be long-lived. In the limit of sin θ 1, 2 → 0, the dominant decay mode of S is the diphoton channel, and its lifetime only depends on its mass m S and the v R scale [cf. Eq. (3.6)]. A long-lived S in the LRSM with lifetime bcτ S 3 m can be a good candidate for the KOTO anomaly. Setting sin θ 2 = 0, the preferred parameter space of m S and the S − h mixing angle sin θ 1 for the KOTO anomaly is shown by the shaded green region in Fig. 2. 4 As in Fig. 1, the dashed green line corresponds to the central value of the KOTO result, while the shaded green band is the 95% CL favored region from Eq. (1.1). The region with S mass m S > 180 MeV is not included here, because it does not have any overlap with the KOTO signal region [6]. For concreteness, we set the VEV ratio ξ = κ /κ = m b /m t which is natural  1). For comparison, we show the exclusion regions from a previous KOTO search for K L → π 0 νν (brown shaded) [3], NA62 search for K + → π + νν (blue shaded) [10], E949 search for K + → π + X (magenta shaded) [11], and beam-dump experiment at CHARM (orange shaded) [12]; cf. Sec. 3.2. The grey and pink shaded vertical regions are excluded by the ∆N eff and supernova constraints respectively (cf. Sec. 3.3). Also shown are the future prospects at NA62 (dashed blue) [68], DUNE (dashed purple) [69], the long-lived particle searches at LHC (dashed red) and FCC-hh (dot-dashed red) [31]. For the CHARM limits we choose four benchmark values of v R = 10, 20, 50 and 100 TeV. Similarly, the solid and dashed lines for the supernova limits correspond respectively to the luminosity of 2 × 10 53 erg and 3 × 10 53 erg for a fixed v R = 10 TeV (cf. Fig. 3).
for the known hierarchy of bottom and top quark masses. We have evaluated the effective branching ratio in Eq. (2.7) for different v R values and found that it is almost independent of the v R value, as in the parameter space of interest the typical lifetime of S is much longer than the KOTO detector size of 3 m. As the FCNC couplings of S are at tree-level in the LRSM, the mixing angles sin θ 1 for the KOTO anomaly (and the following constraints) are orders of magnitude smaller than in the generic singlet model (cf. Fig. 1).

Laboratory constraints
As for the generic singlet model in Section 2, the most stringent limits for the parameter space relevant for the KOTO anomaly are from the searches of K + → π + νν at NA62 [10], K + → π + X at E949 [11], K L → π 0 X at KOTO [3], and the e + e − , µ + µ − and γγ decay products from kaon decay at the CHARM beam-dump experiment [12]. Evaluations of these limits are quite similar to those in the generic singlet model, as discussed in Sec. 2.2 and we do not repeat them here. As in Fig. 1, the current E949, NA62 and KOTO limits are shown respectively by the magenta, blue and brown shaded regions in Fig. 2, and the future NA62 improvement is indicated by the dashed blue lines, which corresponds to a limit down to 2.35 × 10 −11 for the branching ratio BR(K + → π + S) [68]. The limits from the CHARM beam-dump experiment are presented by the orange shaded region in Fig. 2. Unlike the generic singlet model, the most stringent CHARM limit in the LRSM comes from the γγ channel, since this is the dominant decay mode of S for small mixing angles. Therefore the event number depends on the v R scale, as illustrated with four benchmark values of v R = 10, 20, 50 and 100 TeV. For a larger v R value, the lifetime of S tends to be longer and as a result the CHARM limits get weaker.
With an improved proton-on-target number (PoT) of 5 × 10 21 , DUNE can collect 8 × 10 21 kaons, with M pp = 11 and χ s = 1/7 [69]. With the energy E S 12 GeV, the decay length parameters L = 500 m and ∆L = 7 m for the DUNE beam dump set up [69], and setting the Kaon absorption length at DUNE the same as that for CHARM, the current CHARM limits on the mixing angle sin θ can be improved by two orders of magnitude, as shown by the dashed purple curve in Fig. 2. The decay K → πS can also be searched for in the SHiP experiment, but the PoT number 2 × 10 20 is almost one order of magnitude lower than DUNE, and the lifetime that can be probed is also shorter [70]. Thus, we estimate that the prospect of S search at SHiP is weaker than at CHARM and DUNE [31] and is not shown in Fig. 2.
For all the calculations above in the LRSM, we have set the VEV ratio ξ = κ /κ = m b /m t . When the ξ parameter is different, the KOTO region, the NA62 and CHARM limits for sin θ 1 are all universally rescaled by ξ, and this does not help to enlarge the parameter space for the KOTO anomaly. As for the generic singlet model, the limits from the flavor-changing decays K → πχχ (with χχ = e + e − , µ + µ − , γγ) are not applicable to the long-lived S, and the limits from B meson decays, K and B meson oscillations are much weaker than those from the K mesons in the parameter space of interest.
As can be seen from Eqs. (3.4) and (3.5), the decay K → π + S for the KOTO anomaly and the KOTO, E949 and NA62 limits are determined by the scalar mass m S and the mixing angle sin θ 1 , whereas the CHARM limit are mostly from the decay S → γγ which is dictated by the scalar mass m S and the v R scale in the limit of small mixing angles [cf. Eq. (3.6)]. Therefore the LRSM could accommodate the KOTO anomaly while evading the stringent limits from CHARM (and the supernova limits below), which is very different from the singlet scalar model in Section 2.
As in the generic singlet case in Section 2, the long-lived scalar S in the LRSM can also be searched for as LLP in the high-energy colliders [31]. Unlike the singlet scalar case, when the mixing angles sin θ 1,2 are very small (cf. Fig. 2), the scalar S in the LRSM can be produced from the gauge interactions mediated by the heavy W R (and Z R ) bosons, i.e. pp → W R (Z R )S. As a result of the Majorana nature of the heavy RHNs in the LRSM, the smoking-gun signal of W R boson at hadron colliders is same-sign dilepton plus jets without significant missing energy [71], and the current most stringent LHC samesign dilepton limits requires that the W R mass m W R 5 TeV, depending on the RHN mass [72]. 5 For a 5 TeV W R boson, the production cross section of S at LHC 14 TeV is only σ(pp → W R S) 0.025 fb, and we cannot have any LLP prospects for m S < 1 GeV at LHC or MATHUSLA if the SU (2) R gauge coupling is the same as that for SU (2) L (see the left panel of Fig. 17 in Ref. [31]). It is easy to understand: in the limit of small mixing angles sin θ 1,2 , the decay width is proportional to m 3 S /v 2 R (cf. Eq. (3.6)); so for a light S, the decay lifetime is so long that almost no S decays inside the LHC detector. At future 100 TeV colliders FCC-hh and SPPC, the production cross section of S can be almost four orders of magnitude larger than at LHC 14 TeV for m W R = 5 TeV, and we can have LLP prospects for below-GeV scale at FCC-hh and the dedicated LLP detectors therein [31]. Setting m W R = 5 TeV, the LLP prospects at FCC-hh and the dedicated LLP detector is shown in Fig. 2 respectively by the dashed and dot-dashed red lines. The regions to the right side of the red lines can be probed by the LLP searches, which however do not have any allowed KOTO signal region.

Astrophysical and cosmological constraints
As in the generic singlet model case, if S is light, it can be produced in the supernova core and get constrained by the collapse luminosity. In the LRSM, S can be produced in two distinct channels: (i) Nuclear bremsstrahlung process (2.10), which originates from the mixing with the SM Higgs. In this case, the effective couplings of S to nucleons are highly suppressed by the mixing angle sin θ 1 required for the KOTO explanation, thus the corresponding supernova limits are too weak and not relevant to the KOTO anomaly.
(ii) Photon fusion process, i.e. γγ → S, which is highly suppressed by the ratio m 2 S /v 2 R [cf. Eq. (3.6)]. Assuming the photon momentum follows the Bose-Einstein distribution in the supernova core, we follow the calculations in Ref. [74] to estimate the production rate of S which turns out to be just at the order of ∼ 10 53 erg for the benchmark values of v R = (10 − 100) TeV, as shown in Fig. 3. For simplicity, we have set both the scalar mixing angles sin θ 1,2 to be zero. The region of m S for which the luminosity exceeds the observed value of (2−3)×10 53 erg [61] can be excluded. For instance, the supernova limits for v R = 10 TeV are shown by the pink shaded region in Fig. 2, with the solid and dashed lines corresponding to the luminosity of 2 × 10 53 erg and 3 × 10 53 erg respectively, which exclude respectively the mass ranges of 15 MeV m S 27 MeV and 19 MeV m S 23 MeV. If v R goes higher than roughly ∼ 50 TeV, the production rate will be too small such that we do not have any supernova limits.
In the LRSM, even if the mixing angles sin θ 1,2 are extremely small, S can still decay into two photons through the W R boson and the heavy charged scalars. Therefore in the parameter space of interest, the lifetime of S is always much shorter than one second, and we do not have any limits from BBN.
As in the U (1) B−L model case, a light S contributes to the relativistic degree of freedom N eff in the early universe. In the limit of small mixing angles sin θ 1, 2 , S can be in equilibrium with photon if its mass m φ 2 MeV [31]. The current limit of ∆N eff < 0.7 has excluded a scalar particle lighter than 5 MeV [75,76], which is shown by the gray shaded region in Fig. 2.
As can be seen from Fig. 2, although the central value of the KOTO anomaly has been excluded mostly by E949 and NA62, there is still a narrow band left within the 2σ uncertainty of KOTO data between 5 MeV m S 46 MeV with sin θ 1 ∼ (5 − 6) × 10 −8 , even after all the laboratory and cosmological/astrophysical limits are taken into consideration. In addition, the full parameter space can be conclusively tested in the future NA62 and DUNE data.

Conclusion
The three tantalizing events found in the signal region of the flavor-violating decay K L → π 0 νν at the KOTO experiment might be a glimpse of BSM physics. Possible explanation of this by a light long-lived scalar particle which has either tree or loop-level flavor-changing couplings to the s and d quarks and has a lifetime approximately larger than the KOTO detector size of 3 m has been suggested [6]. In this paper, we have studied three possible model implications of this suggestion and constraints on the model parameters from various laboratory measurements and astrophysical/cosmological observations to see if there is any parameter space left to explain this anomaly. In the SM+S model and U (1) B−L model, the whole parameter space explaining the KOTO events has been excluded by the supernova constraint, as shown in Fig. 1. In the LRSM, on the other hand, there still remains a narrow range of parameter space between 5 MeV m S 46 MeV with sin θ 1 ∼ (5 − 6) × 10 −8 which can explain the KOTO anomaly within the 2σ range, while being consistent with all existing constraints, as shown in Fig. 2. This remaining allowed parameter space can be tested in future NA62 and DUNE experiments.
Note Added: While finalizing our manuscript, we noticed Ref. [7], which has some overlap with our SM-singlet scalar case (cf. Sec. 2).