Three Exceptions to the Grossman-Nir Bound

We show that the Grossman-Nir (GN) bound, $\text{Br}(K_L\to \pi^0\nu \bar\nu)\leq 4.3 \, \text{Br}(K^+\to \pi^+\nu \bar\nu)$, can be violated in the presence of light new physics with flavor violating couplings. We construct three sample models in which the GN bound can be violated by orders of magnitude, while satisfying all other experimental bounds. In the three models the enhanced branching ratio $\text{Br}(K_L\to \pi^0+{\rm inv})$ is due to $K_L\to \pi^0\phi_1$, $K_L\to \pi^0\phi_1\phi_1$, $K_L\to \pi^0\psi_1\bar \psi_1$ transitions, respectively, where $\phi_1 (\psi_1)$ is a light scalar (fermion) that escapes the detector. In the three models $\text{Br}(K^+\to \pi^++{\rm inv})$ remains very close to the SM value, while $\text{Br}(K_L\to \pi^0+{\rm inv})$ can saturate the present KOTO bound. Besides invisible particles in the final state (which may account for dark matter) the models require additional light mediators around the GeV-scale.


Introduction
In the SM, the K L → π 0 νν and K + → π + νν decays proceed through the same short distance operator, involving the fields of the quark level transition (s → dνν). The matrix elements for the K L → π 0 νν and K + → π + νν transitions are thus trivially related through isospin, leading to the Grossman-Nir (GN) bound [1] Br(K L → π 0 νν) ≤ 4.3 Br(K + → π + νν). (1.1) The bound remains valid in the presence of heavy New Physics (NP), i.e., for NP modification due to new particles with masses well above the kaon mass. The bound is saturated for the case of maximal CP violation, if lepton flavor violation can be neglected (see Ref. [2] for counter-examples).
In this paper we investigate to what extent NP contributions to K → π+inv decays can violate the GN bound. Simple dimensional counting shows that for large violations of the GN bound the NP needs to be light, of order of a few GeV at most (see Section 2 and Refs. [3,4]). Such light NP faces stringent experimental constraints from rare meson decays and collider/beam dump searches as well as from astrophysics and cosmology. Nevertheless, the couplings needed to modify the rare K → π + inv decays are small enough that interesting modifications of the GN bound are indeed possible. We identify three sample models that achieve this through the following decays: • Model 1: K L → π 0 φ 1 , where the mass of the light scalar, φ 1 , can be anywhere from m φ 1 m K − m π to a few MeV or even less, • Model 2: K L → π 0 φ 1 φ 1 , where the mass of the light scalar, φ 1 , is required in a large part of the parameter space to be above m φ 1 m π /2 in order to avoid constraints from invisible pion decays, • Model 3: K L → π 0 ψ 1ψ1 , with ψ 1 a light fermion whose mass is required to be above m ψ 1 m π /2 in most of the phenomenologically viable parameter space.
The φ 1 and ψ 1 particles are feebly interacting and escape the detector, resulting in the K L → π 0 + inv signature, as does the SM transition, K L → π 0 νν. The NP is thus detected through an enhanced Γ(K L → π 0 + inv) rate. Furthermore, the three models can be distinguished from the SM and each other by measuring the energy distribution of the neutral pion, dΓ(K L → π 0 + inv)/dE π , see Fig. 1 for several sample distributions. While the two body decay in Model 1 results in a fixed pion energy, the three body decays in Model 2 and 3 can be close to the SM distribution for light φ 1 and ψ 1 masses and differ from it for non-negligible masses. Let us mention in passing that the lightness of the scalars could be due to them being a pseudo Goldstone boson of a broken global symmetry whereas for fermions light masses are natural due to chiral symmetry. In all three models the branching ratio Br(K + → π + + inv) remains close to the SM value, Br(K + → π + νν) SM = (8.4±1.0)×10 −11 [5][6][7], and thus below the preliminary NA62 bound Br(K + → π + νν) exp < 1.85 × 10 −10 [8], while Br(K L → π 0 + inv) can be enhanced well above its SM value, Br(K L → π 0 νν) SM = (3.4 ± 0.6) × 10 −11 [5][6][7]. The NP induced K L → π 0 + inv transitions can even saturate the present upper experimental bounds from KOTO Br(K L → π 0 νν) exp < 3.0 × 10 −9 [9]. 1 1 For two body decays the bound is somewhat stronger, Br(KL → π 0 νν)exp < 2.4 × 10 −9 , for m φ 1 mπ [9]. If the new preliminary data are interpreted as a signal, they instead correspond to a rate of Br(KL → π 0 + inv)KOTO = 2.1 +2.0 −1.1 × 10 −9 [10,11] (note though, that some of the observed events are likely due to previously unidentified backgrounds [12]). Furthermore, in the numerics we quote the experimental bounds on three body decays, K → πφ1φ1 K → πψ1ψ1, assuming the experimental efficiencies are the same as for the SM K → πνν transition. In reality, we expect the bounds to be weaker, since the experimental Figure 1. Left: The normalized decay width distributions as functions of the pion energy, E π , for the SM (black line), for the decay dominated by the two body NP transition (Model 1), K L → π 0 φ 1 , for two choices of invisible scalar masses, m φ1 = 1 MeV, 200 MeV (red lines) and for three body NP decay (Model 2), K L → π 0 φ 1 φ 1 , with m φ1 = 1 MeV, 120 MeV (blue lines). Right: The branching ratio distributions, where NP contributions saturate the present KOTO bound [9]. At the kinematic endpoint, E π → m π , we have dΓ SM ∝ p 3 π dE π while for NP Model 2 dΓ K L →π 0 φ1φ1 ∝ p π dE π , where p π = (E 2 π − m 2 π ) 1/2 (similarly for Model 3, Section 5, dΓ K L →π 0 ψ1ψ1 ∼ p π dE π unless the Yukawa couplings y ij are purely real). This follows from the partial wave expansion, dΓ/dE π ∼ p 2l+1 π , adapted to EFTs [22]. In the K → πνν rate the V -A SM interaction induces a negligible Swave contribution proportional to the neutrino mass, whereas the scalar interactions in our model induce a non-suppressed S-wave. The maximum recoil, E π → E π max , in contrast, is controlled by a single power of the ν/φ 1 -velocity in the q µ rest frame, β ν,φ1 = (1 − 4m 2 ν/φ1 /q 2 ) 1/2 (where q 2 = m 2 K + m 2 π − 2m K E π ). For small m ν /φ 1 this velocity goes to 1 for most values of q 2 , leading to a sharp cut-off at E π max .
The three models considered in this work differ from the other proposed NP solutions to the KOTO anomaly in that they allow for large violations of the GN bound at the level of the amplitudes already. In contrast, Ref. [13] relies on the fact that the available phase space is larger for neutral kaon decays due to m K L − m π 0 > m K + − m π + and thus K + → π + X inv decays can be forbidden by a finely tuned choice for the mass of the invisible final state X inv . Ref. [11] instead obtains, in one of the models, an apparent violation of the GN bound from the experimental set-up; the produced light NP particles decay on experimental length-scales, and are not observed in NA62 but are observed in KOTO due to the geometry of the experiments. Finally, the NP models of Refs. [11,[14][15][16][17][18][19][20][21] do not violate the GN bound, but can allow for a large signal in KOTO since NA62 is not sensitive to X inv with a mass close to the pion mass.
The paper is organized as follows. In section 2 a general Effective Field Theory analysis efficiencies are highest for larger values of Eπ, while NP decays considered here are less peaked towards maximal Eπ (as compared to the SM).
is presented. The three models are discussed consecutively in Sections 3, 4 and 5 with the  main plots collected in Figs. 5, 6, in Figs. 15, 16 and in Figs. 19, 20 for Model 1, 2 and 3, respectively, with constraints due to K 0 −K 0 mixing, cosmology and invisible pion decays discussed in the respective sections. The paper ends with conclusions in Section 6, while details on decay rates and integral conventions are deferred to two short appendices.

The EFT analysis
We first perform an Effective Field Theory (EFT) based analysis, assuming that the SM is supplemented by a single light scalar, ϕ, while any other NP states are heavy and integrated out. The light scalar has flavor violating couplings and is created in the K 0 → π 0 ϕ decay. The effective Lagrangian inducing this transition is given by 2 where we only keep the parity-even operators of lowest dimension and work in the quark mass basis. There is a single dimension 4 operator, and the sum runs over the dimension 7 operators, where Γ i , Γ i include both Dirac and color structures. At the quark level the dimension 4 operator induces the s → dϕ transition and thus contributes equally to K 0 → π 0 ϕ and K + → π + ϕ decays, see the first diagram in Fig. 2. The resulting matrix elements for the K L → π 0 ϕ and K + → π + ϕ decays are The K L → π 0 ϕ decay is CP violating and vanishes in the limit of zero weak phases, Im c (4) → 0. These contributions therefore obey the Grossman-Nir relation, The dimension 7 operators, on the other hand, contribute to K 0 → π 0 ϕ and K + → π + ϕ decays in a qualitatively different way. The K L → π 0 ϕ decay can proceed through the weak annihilation type contractions of valence quarks, i.e., through the third diagram in Fig. 2. The K + → π + ϕ transition requires thedd internal line to close in a loop (cf. the 2nd diagram in Fig. 2). Such contractions also contribute to K L → π 0 ϕ. Using at first perturbative counting the latter contributions are suppressed, giving parametric estimates where we neglected m π compared to m K and do not write factors that are parametrically of the same size but may differ by O(1), such as different form factors in the two cases. 2 The full Lagrangian is L = c4 sd ϕ + c 4 sγ5d ϕ + · · · , but we display only the operators that are parity even and thus contribute to the K → π decay.  Figure 2. Contributions from dimension 4 (first diagram) and dimension 7 (2nd and 3rd diagrams) EFT operators to the K → πϕ transition. The last diagram contributes to K L → π 0 ϕ only. The quark-loop diagram in the middle corresponds to the S, E classes of diagrams and the weak annihilation diagram on the right to the W, C classes of diagrams in the lattice computation of Ref. [23].
Depending on the Dirac-color structures Γ ( ) i of the operator one or more gluon exchanges may be required leading to additional (α s /4π) n -factors shown in (2.4).
A priori this leaves two classes of NP models with potentially sizeable violations of the GN bound. The first possibility is heavy NP, with a suppressed c (4) Wilson coefficient such that dimension 7 operators dominate. The other possibility is light NP such that the EFT assumption, on which the above analysis is based on, is violated.
Building viable heavy NP models that violate the GN bound faces several obstacles. First of all, c (4) would have to be heavily suppressed, c (4) m q /Λ, well below naive expectations. If this is not the case, the "heavy" NP scale needs to be quite light. For instance, for c (4) ∼ m q /Λ, c (7) i ∼ O(1) the dimension 4 operator contributions dominate over the dimension 7 ones already for Λ O(3 GeV) (see also the discussion in [3]). Furthermore, even if the hierarchy c (4) c (7) was realised, it is not clear whether the GN bound could be violated by more than a factor of a few. The scaling estimates in (2.4) were based on perturbative expansion, while the kaon decays are in the deep non-perturbative regime of QCD. One can get an idea of the size of the M K ∝ πϕ|O (7) |K matrix elements by linking them to the ones for K → π + − decays that were explored in lattice QCD for light quark masses above their physical value (m π = 430 MeV and m K = 625 MeV) [23]. Figure 5 in Ref. [23] indicates that the quark-loop and weak annihilation contractions, corresponding to the middle and the right diagrams in Fig. 2, lead to contributions of comparable size, contrary to the perturbative expectations in (2.4). If these results carry over to K → πϕ decays, it would seem that the ratio of M (7) (K L → π 0 ϕ)/M (7) (K + → π + ϕ) would not easily exceed a factor of ∼ 2 in models of heavy NP. It is unclear, however, whether this qualitative feature, based on the evaluation of the SM V − A four quark operators [23], would carry over to a model with scalar-scalar four quark operators, originating from a scalar mediator. For instance, for V − A operators the weak annihilation topology is chirally suppressed in the factorisation approximation, while this is not the case for scalar operators.
In conclusion, for heavy mediators the GN bound might or might not be violated in the case c (4) c (7) . In this manuscript we therefore focus on the second possibility, the possibility of light NP mediators, where we can use Chiral Perturbation Theory (ChPT) Figure 3. Diagrams for the K → πφ 1 decay in Model 1 with the GN-violating contribution to the very right. These diagrams enter the matrix elements in Eqs. (3.13), (3.14). Note that the η in the loop contributes to the K L decay only. Diagrams which we neglect, such as the diagrams of O(p 4 ) or O(g 3 qq ), are not shown.
with light NP states as a reliable tool to make predictions.

Model 1 -scalar model leading to two-body kaon decays
In the first example we introduce two real scalar fields, φ 1 and φ 2 . The enhancement of the K → π+inv branching ratio over the SM is due to the K → πφ 1 decay, while K → πφ 2 is kinematically forbidden, i.e., we take m φ 2 > m K − m π . The φ 1 interacts feebly with matter and escapes the detector, resulting in a missing momentum signature 3 . The relevant terms in the Lagrangian are where q, q = {u, d, s} and summation over repeated indices is implied. The couplings g (i) qq are complex, and their imaginary parts trigger the K L → π 0 φ 1 decay. Large violations of the GN bound arise when there is a large hierarchy among the following couplings, g while all other couplings are further suppressed. In our benchmarks these remaining couplings as well as g (1) sd will be set to zero. Before proceeding to predictions for branching ratios and the numerical analysis, it is instructive to perform a naive dimensional analysis (NDA). This will give us insight into why large violations of the GN bound are possible as well as to how large these violations can possibly be.
Taking m S ∼ m φ 2 ∼ m K the NDA estimate for the two decay amplitudes are, sd Im g where the first term in each line is due to the 1st diagram in Fig. 3. The second term in (3.3) is due to the 3rd diagram in Fig. 3, which is absent in the K + → π + φ 1 decay. This is the crucial difference between the two decays and leads to large violations of the GN bound, provided g (1) sd is small. However, violations of the GN bound cannot be arbitrarily large. Even if g (1) sd is set to zero, the K + → π + φ 1 transition is generated at the loop level from the 2nd diagram in Fig.  3, giving the 2nd term in (3.4). Without fine-tuning the ratio M(K L → π 0 φ 1 )/M(K + → π + φ 1 ) is thus at best as large as the loop factor, 16π 2 ∼ 10 3 . Taking into account the present experimental results, this is more than enough to saturate the present KOTO bound while only marginally modifying the K + → π + +inv decay.
In order to simplify the discussion we assume below that the vacuum expectation values (vevs) of the scalar fields φ 1,2 vanish, φ 1 = φ 2 = 0. If this is not the case the K → πφ 1 decays receive additional GN-conserving contributions, see Fig. 4 (right). More precisely, it is the renormalised vevs that are set to zero, φ 1 ren = φ 2 ren = 0, since we work to one loop order. That is, we set the sum of the two diagrams in Fig. 4 to be zero. Had we set them instead to their natural value, φ i ren ∼ m K g (i) dd /16π 2 , our results would not change qualitatively. While M(K + → π + φ 1 ) would be modified by an O(1) factor, in M(K L → π 0 φ 1 ) such contributions are always subleading and one would thus still have large violations of the GN bound.

Estimating the transition rates using ChPT
We use ChPT to calculate the transition rates. In constructing the ChPT we count φ 1 ∼ φ 2 ∼ O(p). 4 As far as QCD is concerned φ 1,2 are external sources and can be treated as spurions [24,25] when building the low energy effective Lagrangian. The QCD Lagrangian, including (3.1), can be conveniently rewritten as, where we keep only the light quarks, q = (u, d, s). The diagonal mass matrix is M q = diag(m u , m d , m s ), while χ (i) S,P are 3 × 3 Hermitian matrices describing the quark couplings to φ 1,2 , Since we set the couplings to the up quark to zero they have the following form 5 (3.7) The off-diagonal couplings in (3.7), where s and p stand for with χ (i) S,P given in (3.7). The LO ChPT Lagrangian, with φ 1,2 included as light degrees of freedom, is given by where the ellipses stand for additional terms in the scalar potential. Here U (x) = exp(iλ a π a /f ) is the unitary matrix parametrizing the meson fields [24,25], B 0 is a constant related to the quark condensate, B 0 (µ = 2 GeV) = 2.666(57) GeV, f is related to the pion decay constant f f π / √ 2 = 92.2(1) MeV [26], with normalization 0|ūγ µ γ 5 d(0)|π − (p) = ip µ f π . The kaon decay constant f K = 155.6 ± 0.4 MeV [27] accommodates SU(3) breaking at times.
In this paper we work to partial NLO order: all LO terms in the chiral expansion O(p 2 ) are kept, as well as the one loop corrections which are of order O(p 4 ) and all finite. 5 For light φ2, which is our preferred scenario, assuming g (2) uu = 0 would not introduce new qualitative features. According to chiral counting, g (2) uu = 0 induces a K + π − φ1-term at O(p 4 ), and is thus subleading to KLπ 0 φ1-terms that we consider. Hence we set g (2) uu to zero for simplicity rather than necessity.
The complete O(p 4 )-expressions for decay amplitudes involves additional contact terms (counter-terms or low energy constants), parametrically of the same size as the one loop corrections. However, since φ 1,2 are propagating degrees of freedom in our EFT the values of the low energy constants in O(p 4 )-ChPT are generally different from the ones in pure QCD and therefore unknown. The associated error in K L → π 0 φ 1 is small, since the NLO corrections are always subleading, while in K + → π + φ 1 they could give O(1) corrections but would not invalidate our conclusions. For simplicity they are set to zero throughout and we do not discuss them any further.
Next we calculate the K → πφ 1 decay amplitudes. Expanding in the meson fields the O(p 2 ) Lagrangian reads where we only kept terms relevant for the calculation of the K → πφ 1 transition, and the analysis of experimental bounds on the φ 1 -couplings.
The NP contributions to the decay amplitude for the K L → π 0 φ 1 and K + → π + φ 1 transitions are, see Fig. 3, are structures occurring in all three models. They depend on the loop function Moreover, the replacement f 2 → f π f K /2 accounts for the main SU(3) breaking effects.
Note that the amplitude vanishes in the limit of no CP violation, Imĝ sd → 0. The first term in (3.13), proportional toĝ (2) sd , is the O(p 2 ) contribution due to the tree level exchange of φ 2 , see the 3rd diagram in Fig. 3. It is isospin violating since it gives rise to the K L → π 0 φ 1 transition but not to K + → π + φ 1 . The first term in the second line of Eq. (3.13) is the remaining O(p 2 ) contribution, due to the emission of φ 1 directly from the meson line, see the 1st diagram in Fig. 3. This contribution is isospin conserving -it is present for both K L → π 0 φ 1 and K + → π + φ 1 transitions. It is proportional toḡ (1) sd and is thus small due to the assumed hierarchy among the couplings, Eq. (3.2).
The hierarchy of couplings |g (2) sd,ds | |g (1) sd,ds | thus leads to maximal violation of the GN bound by NP contributions. However, this violation cannot be arbitrarily large. Even in theḡ (1) sd → 0 limit we still have isospin conserving NP contributions generated at one loop, see the 2nd diagram in Fig. 3, giving the last term in (3.13). If φ 2 is heavy and integrated out these radiative corrections match onto the φ 1 − Kπ vertex, which is then radiatively induced. Moreover the K L → π 0 φ 1 and K + → π + φ 1 decays receive contributions from π 0 − φ 1 mixing where flavor violation comes from the SM K → ππ transition. For our choices of parameters these contributions are always negligible.
The NP contributions add coherently to the SM rate, 17) and the partial decay width due to NP is is the pion's momentum in the K L rest frame and λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz the kinematic Källén function. The expressions for the K + → π + + inv decay is completely analogous. Numerically, this gives (the SM predictions are taken from Refs. [5][6][7]28]) where we kept only the leading term for the NP contribution. The typical values of the inputs parameters for the NP contribution were chosen such that they reproduce roughly the KOTO anomaly (in fact slightly larger, but within 1σ). Note the very high scaling in the φ 2 mass, underscoring that φ 2 needs to be relatively light in order to have large violations of the GN bound. For the charged kaon decay the numerical result is where in the NP contribution we only kept the tree level term and set the value ofḡ (1) sd to be similar to the one-loop threshold correctionḡ (1) sd ∼ḡ (2) sd g (2) dd /8π 2 , cf. Eq. (3.14), with the typical values of the later couplings as in (3.19). While the correction to K + → π + + inv is O(1) of the SM branching ratio, the correction to K L → π 0 + inv can be orders of magnitude above the SM, giving large violations of the GN bound. Note that NP in Model 1 contributes to the 2-body decay K + → π + +X 0 only, and for massless X 0 is subject to the bound Br(K + → π + + X 0 ) < 0.73 × 10 −10 from E949 [29], which is slightly stronger than the preliminary NA62 bounds on the 3-body decay Br(K + → π + + inv) exp < 2.44 × 10 −10 and the 2-body decay Br(K + → π + + X 0 ) exp 1.9 × 10 −10 (for massless X 0 ) [8].

Constraints onĝ
The K 0 −K 0 mixing is an important constraint on the model. The contributions to the meson mixing matrix element are where the tree-level exchanges of φ 2 is with the ellipses denoting higher order terms (we also neglect the NP contributions to the absorptive mixing amplitude since it only enters at one loop). The replacement f → f K / √ 2 accounts for the SU(3) breaking.
We consider two constraints, ∆m K and K which are CP conserving and CP violating respectively. Using the relation ∆m K = 2ReM 12 and conservatively assuming, due to the relatively uncertain SM predictions of ∆m K , that the NP saturates the measured ∆m K , we obtain in the limit m φ 2 m K , and with the experimental value ∆m expt. K = 3.484(6) × 10 −12 MeV [27], this translates to To obtain the bounds on non-SM CP violating contributions to K 0 −K 0 mixing we use the normalized quantity For the theoretical prediction of K we use the expression [30] K = e iφ sin φ We take the values for ∆m K = m L − m S , ∆Γ K = Γ S − Γ L , and φ = arctan(2∆m K /∆Γ K ) from experiment [27]. With the SM prediction for | K | from [31], and the NP contribution to M 12 , Γ 12 from Eq. (3.22) we get The global CKM fit by the UTFit collaboration results in 0.87 < C K < 1.39 at 95% CL [32,33], which translates to the following 1σ bounds (3.29) These bounds will improve in the future, once the improved prediction for K [31] is implemented in the global CKM fits.

Constraints on representative benchmarks
To highlight the typical values of couplings that can lead to sizable correction in K → π+ inv decays, while passing all other constraints, we form a benchmark 1 (BM1) and a benchmark 2 (BM2), These depend on two real parameters, g dd and g sd , parametrizing couplings of φ 2 to quarks. All the remaining couplings of φ 2 to quarks as well as all the direct couplings of φ 1 to quarks are set to zero in accordance with previous discussions. The triple scalar coupling is fixed to λ S m S = 1 GeV (and other potentially relevant scalar couplings assumed to be small, see Section 3.4.1). The mass of φ 1 is taken to be small, m φ 1 = 1 MeV, while m φ 2 is kept as a free parameter that is varied in the range m φ 2 ∈ [0.4, 1.5] GeV, cf. footnote 4. The form of couplings in BM1, Eq. (3.30), is such that the NP contributions to K are maximized. This benchmark is thus representative of the parameter space that is most constrained. Fixing g dd = 10 −3 the allowed regions are shown in Fig. 5. The red regions are excluded by the NA62 bound on Br(K + → π + φ 1 ) exp 1.9 × 10 −10 [8], the E949 bound Br(K + → π + φ 1 ) exp < 0.73 × 10 −10 [29] and by the KOTO bound Br(K L → π 0 φ 1 ) < 2.4 × 10 −9 [9]. The E949 and NA62 bounds shown are for massless φ 1 , which is a good approximation for our benchmarks, where m φ 1 = 1 MeV. For heavier masses, above m π , the bound is expected to become significantly weaker and completely disappear for m φ 1 m π 0 , as in [34]. The green bands denote the 1σ bands of the branching ratio [10,11] that corresponds to the anomalous KOTO events. The blue line denotes the GN bound, showing that large violations of the GN bound are possible in this model. This violation is most apparent in Fig. 5 (right) which gives the allowed values of g sd as a function of m φ 2 , with the dashed lines denoting contours of the ratio Br(K L → π 0 + inv)/Br(K + → π + + inv). The present KOTO bound is saturated by values for this Figure 5. The parameter space for Model 1, BM1, for g dd = 10 −3 in (3.30). The GN bound is denoted with blue lines, while the green regions give the 1σ bands corresponding to KOTO anomalous events [10,11]. Left: the predictions for Br(  [29] and KOTO [9]. Right: Contours of Br(K L → π 0 + inv)/Br(K + → π + + inv) (dashed lines) as functions of g sd , m φ2 , with the hatched regions excluded by K 0 −K 0 mixing and π 0 → inv bounds. The region around the kaon mass is masked out (gray region).
ratio of around 20, while still satisfying the K constraint, Eq. (3.29), and the π 0 → inv constrain discussed below, see Eq. The solid black lines in Fig. 5 (left) show the values of Br(K L → π 0 +inv) and Br(K + → π + + inv) for g sd = 5 · 10 −10 and g sd = 2 · 10 −9 , varying m φ 2 ∈ [0.4, 1.5] GeV, while fixing g dd = 10 −3 (the grey dotted parts of the lines are excluded by a combination of K 0 −K 0 and π 0 → inv constraints). The SM predictions for the two branching ratios, Br(K + → π + νν) SM = (8.4 ± 1.0) × 10 −11 and Br(K L → π 0 νν) SM = (3.4 ± 0.6) × 10 −11 [5][6][7], are denoted with blue bands (1σ ranges). For the larger value, g sd = 2 · 10 −9 , the prediction is still quite far away from the SM for m φ 2 = 1.5 GeV, but would of course tend to the SM for m φ 2 → ∞. For larger values of g sd deviations from the SM prediction for Br(K + → π + + inv) at the level of a few are predicted for this benchmark and subject to the indicated constraints from E949, while for smaller values of g sd the deviations in Br(K + → π + + inv) become negligibly small. That is, it is possible to explain the KOTO anomalous events without having any appreciable NP effects in the charged kaon decay nor in K 0 −K 0 mixing. We next move to BM2. The form of couplings in Eq. (3.31) was deliberately chosen such that there is no NP CP violation in K 0 −K 0 mixing, in order to avoid the K bound. The bound from ∆m K , Eq. (3.24), on CP conserving contributions to K 0 −K 0 mixing is much weaker, giving the hatched excluded region in Fig. 6 (right). This means that for the same mass of φ 2 the flavor violating couplings to quarks can be much larger than in BM1. In Fig. (right) 6 we show the g dd = 3 × 10 −5 slice of the parameter space, in which case g sd can be as large as 10 −7 . Furthermore, the form of couplings in BM2, Eq. (3.31), is such that there is no NP effect at all in Br(K + → π + +inv), to the order we are working, and the E949 bound is completely avoided. In contrast, the effect on Br(K L → π 0 +inv) can be very large and easily saturate KOTO's present upper bound, as shown for two representative couplings g sd = 3 × 10 −9 , 8 × 10 −8 (black lines, with dashed parts of the lines excluded by ∆m K ). BM2 comes with enhanced symmetry; φ 2 is a pure pseudoscalar and φ 1 a pure scalar. This has implications for flavor conserving couplings of φ 1 , to which we turn next.

Constraints on the φ 1 -couplings
So far the scalar mass φ 1 has been fixed to 1 MeV. Next, we show that for the two benchmarks the radiatively generated couplings of φ 1 to pions, nucleons, and photons are all well below the bounds for a large range of φ 1 masses (including m φ 1 = 1 MeV). Figs. 5 and 6 are thus valid for a larger set of φ 1 masses, as long as m φ 1 m K .
Finally, the π 0 → φ 1 φ 1 decay could also proceed at tree level via an additional interaction term in (3.11) of the form δL = λ m S φ 2 φ 2 1 . Whereas, contrary to Model 2, λ plays no role in the K → πφ 1 decays per se, it is potentially dangerous for the invisible pion decay. In the absence of a UV completion we may choose its initial value to be sufficiently small (zero in practice) to pass the constraint.

φ 1 − π 0 mixing
The φ i mix with light pseudoscalars through the g (i) qq couplings, Eq. (3.1). The φ 1 − π 0 part of the mass matrix to one loop receives contributions in Fig. 7, and is parametrized by the Lagrangian, m φ 2 m φ 1 ,π,η , with the effective φ 1 − π 0 coupling given by dd Re g dd L(m π ) ss Re g dd L(m η ) + g where we have exceptionally kept the g ds -terms since they are leading in BM2. The first term is due to tree level mixing, see Fig. 7 (left), the second term are the one loop corrections due to diagram in Fig. 7 (right). The loop function L(m X ) ≡ −m 2 In the following we will take for simplicity this limit, which provides a reasonable approximation for the parameter region of Figure 7. The leading order and one loop induced φ 1 − π 0 mixing. Figure 8. CP violating contributions to φ 1 → γγ, matching onto the coupling g 1γγ .
interest, since L π 0.8, L K,η 0.4 for m φ 2 = 400 MeV. In the two benchmarks (3.30), (3.31), the effective φ 1 − π 0 couplings are In BM2 there is no φ 1 − π 0 mixing φ 1 is a pure scalar and parity is conserved. Working in the mass insertion approximation for the off-diagonal mass term, Eq. (3.36), the φ 1 − π 0 mixing angle, s θ ≡ sin θ ≈ θ, between the interaction states φ 1 and the mass eigenstate Note that this expression for the mixing angle is only valid for m φ 1 sufficiently far away from m π . For the two benchmarks, we have The φ 1 − π 0 mixing is thus very small in most of the viable parameter space, justifying the use of the mass insertion approximation.

Couplings of φ 1 to photons
The dominant decay channel of φ 1 is to two photons. In the limit m φ 1 m π the interactions with two photons are described by the effective Lagrangian The dominant contribution to the CP violating coupling g 1γγ is from the π 0 anomaly term Figure 9. CP conserving contributions to φ 1 → γγ, matching onto the coupling h 1γγ . Figure 10. Tree level and one loop contributions matching onto the effective couplings g 1ππ /g 1KK .
via the φ 1 − π 0 mixing, see Fig. 8. Working in the mass insertion approximation for the off-diagonal mass term, Eq. (3.36), gives with g 1π given in (3.37). The CP conserving h 1γγ coupling receives the first relevant contributions from radiative corrections with K + and π + running in the loop cf. Fig. 9. For our benchmarks the first nonzero contributions arises at two loops, while for BM2 the numerically most important contribution arises at three loops In the m φ 2 m K (m π m φ 1 by assumption) limit the one and two loop contributions, in Fig. 9, assume the form whereas the effective couplings of φ 1 , L eff ⊃ φ 1 (g 1ππ π + π − + g 1KK K + K − ), to two light charged mesons evaluate to and g 1KK = g 1ππ | dd→ss . The first term in (3.47) is the tree level term from (3.12), see Fig. 10 (left). In both benchmarks, BM1 and BM2, this contribution was set to zero. The second term in (3.47) is the one loop correction, see Fig. 10 (right). We kept the flavor violating contribution proportional toḡ sd even though it is numerically negligible. For the three loop contribution to h 1γγ we resort to a NDA estimate, still in the m φ 2 m K limit, dd )(Im g (2) ss ) , where the O(1) factors are not displayed. Finally we are in a position to assemble the results for the benchmarks. Using g sd g dd , the φ 1 -photon couplings evaluate to in BM1, while for BM2 they turn out to be and we remind the reader that λm S = 1 GeV for reference. The g 1γγ coupling vanishes in BM2 since φ 1 is a parity even scalar in that benchmark. The value quoted for h BM2 1γγ is the NDA estimate of the flavor conserving 3 loop contribution. For representative values of g dd in the two benchmarks we used the values in Figs. 5 and Fig. 6 for BM1 and BM2, respectively.
The above couplings of φ 1 to photons are sufficiently small that for both benchmarks the φ 1 is stable on collider scales. More concretely, the φ 1 → γγ partial decay width is given by such that φ 1 is stable on solar to cosmological timescales. For such small couplings the laboratory constraints from, e.g., π + → φ 1 e + ν decays [37] are irrelevant, whereas astrophysical and cosmological constraints are important (cf. figure 1 in Ref. [38]) and further discussed in Section 3.4.5.

Couplings of φ 1 to nucleons
The couplings of φ 1 to protons and neutrons are tree-level and loop-level induced by g (1) dd and g (2) dd respectively, cf. Fig. 11. One can use Heavy Baryon Chiral Perturbation Theory (HBChPT) [39] to organize different contributions. We only keep only the leading terms which are (in relativistic notation) with t a = σ a /2, a = 1, 2, 3 and σ a are Pauli matrices, N = (p, n) the isospin doublet of nucleons, and In the heavy φ 2 limit the following effective Lagrangian provides a good description of the φ 1 -nucleon system. Assuming m φ 2 ,N m φ 1 , the diagrams in Fig. 11, evaluate to where Y N 1 stands for the nucleon-nucleon entries in (3.54). In theF (r) term in (3.56) we in addition assumed the m π m φ 1 limit. The real-valued loop functions F (r),F (r), with r = m 2 φ 2 /m 2 N , are given by 6 The first term in (3.56) is due to the 1st diagram, while the one loop corrections are due to the 2nd and 5rd diagram in Fig. 11. For the pseudoscalar coupling to nucleons,g 1N N , we keep the pion exchange term (dropping the η-exchange) in the 4th and the 5th diagram in Fig. 11 resulting in the tree level and one loop terms in (3.57). To simplify the expressions we show the one loop contribution in (3.57) only in the heavy m φ 2 limit.
Numerically, we have for BM1, setting m φ 2 = 1 GeV, 60) 6 It is noted that the 3rd diagram, in Fig. 11, does not introduce any infrared (IR) divergences in the limit mπ → 0. This is a consequence of the derivative couplings of pions, cf. Eq. (3.53). We note in passing that for a double insertion of this interaction term one cannot use the naive EOM and replace gA(N γ µ γ5t a N )∂µπ a → −2mN gAN γ5t a N π a . For a concise technical discussion we refer the reader to Ref. [41]. Use of the naive EOM leads to the IR divergence that is linked to the absence of the derivative coupling in that case. The same applies to the single insertion of the gA-term in 4th and 5th diagram. Figure 11. The leading order and one loop induced φ 1 -couplings to nucleons grouped into parity conserving coupling g 1N N (even in g A ) and parity violating couplingg 1N N in (3.56) (odd in g A ). The 3rd diagram is the only non-vanishing contribution to g 1N N in BM2. Below we analyse the combined constraints from the previous two subsections.

Combined analysis of φ 1 -constraints
The most important constraint on the φ 1 -couplings comes from the neutrino burst duration observed in the supernova SN1987A. The interactions of φ 1 with matter inside an exploding supernova are dominated by its couplings to nucleons. For m φ 1 = 1 MeV, used in our benchmarks, the SN1987A observations exclude g eff 1N N ≡ (g 2 1N N + (3/2)g 2 1N N ) 1/2 in the range 7 · 10 −10 GeV −1 g eff 1N N 4 · 10 −6 GeV −1 [42]. For larger values of g eff 1N N the φ 1 gets trapped inside the proto-neutron star (PNS) and does not contribute to the cooling. This is the case for BM1, see Eqs. The photon couplings of φ 1 are less relevant for SN1987A since the Primakoff emission of φ 1 is always subdominant relative to the emission of φ 1 in nucleon-nucleon scattering. Figure 13. The constraints on the g dd coupling in BM1 (left) and BM2 (right) due to couplings of φ 1 to photons and nucleons as a function of the φ 1 mass. The purple regions are excluded by beam dump searches, E949 (K + → π + X) and NA62(π 0 → inv), the red region by SN1987, while the dashed line shows the upper bound from cosmology in the absence of any other light states or φ 1 -couplings. The star denotes the values of g dd and m φ1 in Fig. 5 (Fig. 6) for BM1 (BM2). The region around m φ1 m π 0 is masked out (gray region). This is best illustrated by the fact that SN1987A would exclude the range 10 −8 GeV −1 g 1γγ , h 1γγ 10 −5 GeV −1 , if φ 1 were to coupled to photons only. The induced couplings of φ 1 to photons are at the lower edge of this range for BM1 and well below for BM2, cf. Eqs. (3.48) and (3.49) respectively. This should be contrasted with nucleon couplings which for BM1 traps φ 1 inside the PNS as it is above and not below the exclusion window.
The constraints from the SN1987A neutrino burst duration are shown for a range of φ 1 masses for benchmarks BM1 and BM2 in Fig. 13 (left) and (right) as red regions, respectively. According to the analysis of Ref. [42], the bounds are relevant all the way up to m φ 1 300 MeV, though we truncate the plots at 200 MeV. These bounds may however depend on the details of the SN1987A explosion, and may even be absent if this was due to a collapse-induced thermonuclear explosion [43].
In addition, Fig. 13 shows with purple shading the constraints from beam dump experiments (we use the combined limit as quoted in [42]), and from the invisible pion decay by NA62 [8]. The φ 1 − π 0 mixing angle s θ needs to be smaller than about 2 × 10 −5 in order to satisfy the K + → π + X constraints from E949 [29] and NA62 [8]. This imposes a constraint on g dd that is comparable but slightly less stringent than the beam dump limit. The upper bound from cosmology, i.e., the impact of φ 1 decays on big bang nucleosynthesis and distortions of cosmic microwave background, are shown with a dashed line [44]. This bound is very sensitive to the details of the model. For instance, if the φ 1 decays predominantly to neutrinos these bounds would be drastically modified and thus potentially irrelevant. Figure 14. The diagrams inducing the K → πφ 1 φ 1 decays in Model 2, with the matrix elements shown in Eqs. (4.3) and (4.4). The 3rd diagram violates the GN bound.

Model 2 -scalar model leading to the three-body kaon decays
Model 2 has the same field content as Model 1, except that we impose a Z 2 symmetry under which the scalar φ 1 is odd, φ 1 → −φ 1 . The relevant terms in the Lagrangian are Note that the coupling (q L q R )φ 1 is forbidden by the Z 2 -parity. Because of the Z 2 parity the φ 1 always appears in pairs in the final state and we thus focus on the K → πφ 1 φ 1 transitions with leading diagrams shown in Fig. 14. The 1st diagram in Fig. 14, proportional to the trilinear coupling λ , gives the same contribution to both K + → π + φ 1 φ 1 and K 0 → π 0 φ 1 φ 1 transitions in accordance with isospin. Since we are interested in violations of the GN bound, we impose the hierarchy λ , λ λ 4 , (4.2) and assume m S = O(m K ). For simplicity we further assume that φ 1,2 do not have vevs, or that they are negligibly small (cf. related discussion for Model 1 in Section 3).
Keeping the leading diagrams in the λ and λ 4 g dd expansion, i.e., the diagrams in Fig.  14, the K L → π 0 φ 1 φ 1 decay amplitude reads with F (2) L given in (3.15), while the K + → π + φ 1 φ 1 decay amplitude is , and q 2 = (p 1 + p 2 ) 2 is the invariant mass squared of the φ 1 φ 1 final state system. As for Model 1, f 2 → f π f K /2 in order to account for the main SU(3) breaking effect.
The structure of the two decay amplitudes is reminiscent of the results in Model 1 in Eqs. (3.13), (3.14). The main difference is that there is no direct coupling of φ 1 to quarks due to the Z 2 symmetry. The φ 1 φ 1 pair couples to d → s current instead through the off-shell tree level exchange of φ 2 , see the 1st diagram in Fig. 14. This leads to isospin symmetric contributions to K + → π + φ 1 φ 1 and K L → π 0 φ 1 φ 1 , proportional to the trilinear λ coupling. Hence, in the λ → 0 limit, the K + → π + φ 1 φ 1 transition only receives loop contributions, and the GN bound is maximally violated. Note that λ cannot be arbitrarily small, since it is generated at one loop through φ 2 loop, λ ∼ λ 4 λ /(16π 2 ), and at two loops with φ 2 and π 0 , η running in the loop: λ ∼ λ 4 (g (2) dd ) 3 /(16π 2 ) 2 . For our benchmarks this gives a vanishingly small λ and thus this contribution can be safely ignored in our analysis provided the bare value of λ , λ are set to zero. In this limit the first isospin conserving contribution is at one loop due to the 2nd diagram in Fig. 14. The GN-violating contribution instead arises at tree level, see the 3rd diagram in Fig. 14 and the first term  in (4.3).
The total rate of K L → π 0 φ 1 φ 1 adds coherently to the SM K L → π 0 νν rate. The differential rate for K L → π 0 φ 1 φ 1 is given by where p π = E 2 π − m 2 π and E π = (m 2 K + m 2 π − q 2 )/(2m K ) are the pion's momentum and energy in the K L rest frame, while

Benchmarks for Model 2
The bounds from K 0 −K 0 mixing on flavor violating φ 2 -couplingĝ (2) sd are exactly the same as for Model 1, Section 3.2. To illustrate the available parameter space we therefore use the same two benchmarks for the φ 2 -couplings, Eqs. while m φ 2 is kept as a free parameter. The results in Figs. 15, 16 are fairly independent of the φ 1 mass as long as it is taken to be small, m φ 1 m K , and thus does not modify the final phase space. The choice of benchmark value m φ 1 = 100 MeV is driven by the constraints of the invisible pion decays, see Section 4.2. BM1 and BM2 thus have three free parameters: g dd , g sd and m φ 2 . BM1, shown in Fig. 15, has a well restricted {g sd , m φ 2 } parameter space, since the tree level exchanges of φ 2 contributes a new CP violating source in K 0 −K 0 mixing. This then restricts g sd to be below the hatched region in Fig. 15 (right), see also Eq. (3.29). However, large enhancements of Br(K L → π 0 + inv) over Br(K + → π + + inv) are still possible in significant parts of the parameter space. For instance, setting g dd = 5 · 10 −2 , the KOTO upper bound Br(K L → π 0 νν) exp < 3.0 × 10 −9 [9] (red region in Fig. 15) are obtained for g sd O(10 −9 ) and m φ 2 O(1 GeV). Fig. 15 (left) shows that in the relevant region of parameter space the deviations in Br(K + → π + + inv) from the SM prediction are Figure 15. The parameter space for Model 2, BM1, Eq. (4.6). The color coding is the same as in Fig. 5. In the predictions for Br(K + → π + + inv), Br(K L → π 0 + inv) on the left plot (black lines) we vary m φ2 ∈ [0.55, 1.5] GeV for two values of g sd = 2 × 10 −10 , 2 × 10 −9 .
negligible, while Br(K L → π 0 + inv) can be enhanced well above the GN bound (blue line). For illustration we vary m φ 2 ∈ [0.55, 1.5] GeV, the same range as is shown in Fig. 15 (right), fix g dd = 5 · 10 −2 and show predictions for two choices of g sd = 2 · 10 −10 , 2 · 10 −9 (black lines). The resulting range in Br(K L → π 0 + inv) is denoted with arrows. For g sd = 2 · 10 −9 the K bound is reached, and the exclusion range is shown with gray dotted lines. For both choices of g sd the enhancements can easily be in the range of the KOTO anomaly (green band) without violating any other bounds.
For λ 4 ∼ O(1), which is required for large violations of the GN bound, this excludes φ 1 masses m φ 1 m π /2 for BM1. In BM2 the π 0 → φ 1 φ 1 is forbidden due to parity, so that φ 1 can be light, as long as the parity breaking term λ m S is sufficiently small.
The beam dump and SN constraints in BM1 and BM2 are very similar to the ones shown in Fig. 13 for Model 1, but with rough identification g dd Model 1 → 1/(4π) × g dd Model 2 and m φ 1 Model 1 → 2m φ 1 Model 2 , since the transitions now involve two φ 1 particles in the final state. In particular, for the choices of parameters in Figs. 15 and 16 the collider and SN constraints are presumably satisfied. Figure 17. Graph dominating the annihilation cross section of φ 1 φ 1 → π 0 π 0 and ψ 1ψ1 → π 0 π 0 in Model 2 and Model 3 respectively.

φ 1 as a dark matter candidate
Since φ 1 is odd under the Z 2 -parity, it is absolutely stable and could be a dark matter (DM) candidate. If m φ 1 > m π 0 the φ 1 φ 1 → π 0 π 0 annihilation channel is open. Below we shall see that the annihilation cross section is large enough, in part of the parameter space, such that the correct DM relic abundance is obtained. Note however, that this restricts φ 1 to a rather narrow mass range, The annihilation cross section for φ 1 φ 1 → π 0 π 0 process is dominated by the λ 4 vertex and φ 2 − π 0 -conversion, while the φ 2 s-channel resonance contribution is subleading, since λ λ 4 , cf. (4.2) and Fig. 17. Assuming a non-relativistic φ 1 , as is the case at the time of freeze-out, the leading thermally averaged cross section is given by 10) where in this approximation p π = (m 2 ψ 1 − m 2 π ) 1/2 . Taking m φ 1 = 160 MeV as a representative value gives which is of the right size to get the correct DM relic abundance (3 · 10 −26 cm 3 /s ≈ 1 pb).
For m φ 1 < m π 0 the φ 1 φ 1 → π 0 π 0 annihilation cross section is kinematically forbidden. In that case the dominant annihilation channel becomes φ 1 φ 1 → γγ. The resulting annihilation cross section is so small, that if this were the only annihilation channel, the φ 1 would overclose the universe [45]. This means that φ 1 should also couple to other light states. For instance, φ 1 could annihilate into light SM particles, e.g. φ 1 φ 1 → e + e − or φ 1 φ 1 → νν. Alternatively it could annihilate away to other light dark sector particles or dark photons φ 1 φ 1 → γ D γ D (if φ 1 was gauged under a dark U (1)). Since none of these couplings are related to K → πφ 1 φ 1 decays we do not explore the related phenomenology any further, beyond stating the obvious -that φ 1 could well be a thermal relic for appropriate values of these additional couplings. Figure 18. The K → πψ 1ψ1 in the fermion Model 3 with contribution evaluated in Eqs. (4.3) and (5.8) respectively. The third diagram only contributes to K L → π 0 ψ 1ψ1 , as the notation suggests, and is therefore responsible for potential violation of the GN bound in Model 3.

Model 3 -light dark sector fermions
In this model we introduce a real scalar, φ, of mass m φ , and two Dirac fermions, ψ 1 , ψ 2 , with masses m ψ 1 ,ψ 2 , where the couplings relevant for the K → π+inv decay are The fermion ψ 2 is massive enough such that the decays of K → πψ 2ψ2 and K → πψ 1ψ2 are kinematically forbidden. In contrast and crucially, the decay K → πψ 1ψ1 is assumed to be kinematically allowed. The couplings of φ to the quarks are assumed to have a hierarchical flavor structure g dd,ss , sd,ds can be arbitrarily suppressed in accordance with (5.2). Keeping the leading diagrams in y ij and g (φ) dd , shown in Fig 18, gives the following K L → π 0 ψ 1ψ1 decay amplitude where 2P R,L ≡ 1 ± γ 5 , we have shortenedū ≡ū(p), v ≡ v(p), while q 2 = (p +p) 2 is the invariant four momentum of the fermion pair. The∆ (µ) stands for combinations of fermion propagators∆ , with the latter defined in (3.15). Its arguments are given in terms of loop integrals, where is the same integral with an additional Lorentz-vector k µ in the integrand.
The decay amplitude K + → π + ψ 1ψ1 is analogous, but without the 3rd diagram in Fig. 18. This gives where F (φ) with the later defined in (3.16). A formula for the rate, in differential form, is given in Appendix A.

Benchmarks for Model 3
The new elements of Model 3 are the Yukawa couplings between φ and ψ 1,2 , as well as the absence of the light-scalar φ 1 . In order to ease comparisons with Model 1 and Model 2, we use g while all the other couplings are set to zero. In particular, the only nonzero Yukawa couplings of ψ i fermions for φ are the flavor violating ones, y 12,21 , while the diagonal ones ����� �� ���� ���=�·�� -� � �ψ � =������ Figure 19. The parameter space for Model 3, BM1, Eq. (5.9). The color coding is the same as in are assumed to be vanishingly small, and set to y 11,22 = 0. The mass of the lightest fermion is set to m ψ 1 = 100 MeV. The benchmarks are thus described by four continuous variables: the masses m φ , m ψ 2 and the real parameters g dd , g sd .
The flavor violating coupling,ĝ (φ) sd , is constrained by K 0 −K 0 mixing. The bounds are the same as for φ 2 in Model 1, Section 3.2, and are thus obtained from Eqs. (3.29), (3.24) through theĝ (2) sd →ĝ (φ) sd , m φ 2 → m φ replacements. BM1 is severely constrained by K since tree level exchange of φ induces a new CP violating contribution to K 0 −K 0 mixing. Fig.  19 (middle) and (right) show that large enhancements of Br(K L → π 0 + inv)/Br(K + → π + + inv) are possible only for small values of m φ and m ψ 2 , comparable to the kaon mass. Still, such light NP states are not excluded experimentally and can saturate the present KOTO bound. Fig. 19 (left) shows that in this regime it is possible to have values for this ratio well above the GN bound, in the range of the anomalous KOTO events (green band).
BM2, on the other hand, does not lead to tree level contributions to K . The constraints from K 0 −K 0 mixing are therefore relaxed compared to BM1 as they are only due to ∆m K . As shown in Fig. 20, it is thus possible to saturate the present KOTO upper bounds over a much larger set of parameter space, with masses of m φ and m ψ 2 up to ∼ 1 GeV for g dd = 5 · 10 −2 . Next, we discuss the constraints on the ψ 1 -couplings.
whereM is a shorthand for 2 quoting H(0) = 1 and H(1) = 0.5 as representatitve values. The total rate is easily obtained from the matrix element squared given in (A.4) and the 1 → 2 decay rate (3.33) (without the symmetry factor 1/2) replacing f → f π / √ 2 and adapting β ψ 1 = (1 − 4m 2 ψ 1 /m 2 π ) 1/2 . Assuming Im g which are close to the upper experimental bound Br(π 0 → φ 1 φ 1 ) < 4.4 × 10 −9 [8]. Clearly the tree graph is leading and imposes a constrain |y 11 | O(5 · 10 −6 ) on the Yukawa couplings to the light fermion for the two benchmarks in Figs. 19 and 20. We observe that, for both BM1 and BM2 the lightest fermion is required to be heavy enough that invisible pion decay is kinematically forbidden, m ψ 1 m π 0 /2. Reducing somewhat the value of |y 12 y 21 | ∼ O(10 −2 ), very light ψ 1 are possible. Even in this case the KOTO bounds could be saturated (at least for the BM2 flavor structure of the couplings).

ψ 1 as a dark matter candidate
In the minimal version of Model 3, presented in this work, ψ 1 and ψ 2 are odd under the Z 2 -parity. The lightest fermion, ψ 1 can thus be a DM candidate. The situation is similar to Model 2. For ψ 1 in the mass range 135 MeV m ψ 1 181 MeV, ψ 1ψ1 → π 0 π 0 is kinematically allowed, and can lead to the correct relic DM abundance. For lighter ψ 1 only the ψ 1ψ1 → γγ annihilation is allowed. However, if φ were to couple to electrons or neutrinos, the resulting annihilation cross sections can be large enough such that ψ 1 can be the DM.
For now, let us assume that ψ 1ψ1 → π 0 π 0 is kinematically allowed. Then at leading order there are two relevant diagrams as shown in Fig. 17. The corresponding matrix element reads wherev ≡v(p), u ≡ u(p), s ≡ q 2 = (p +p) 2 , and by crossing symmetry from the right diagram in Fig. 18:∆ (µ) =∆ (µ) | p,p→−p,−p in (5.6). The cross section is obtained from the spin-averaged squared matrix element (including a symmetry factor for identical final states) where β π,ψ 1 = (1 − 4m 2 π,ψ 1 /s) 1/2 are the respective velocities in the centre of mass frame. The thermally averaged cross section at leading order in the non-relativistic expansion is given by 17) where in this approximation p π = (m 2 ψ 1 − m 2 π ) 1/2 , and which is of the right order of magnitude to produce the required relic abundance (1 pb ≈ We have presented three models that can lead to large deviations in Br(K L → π 0 + inv), while leaving Br(K + → π 0 + inv) virtually unchanged from the SM expectation. The three models are: Model 1 where the invisible decay is the two body transition K L → π 0 φ 1 , Model 2 with K L → π 0 φ 1 φ 1 and Model 3 with K L → π 0 ψ 1ψ1 three body transitions, can be viewed as representatives of a larger class of models. The scalar φ 1 or fermion ψ 1 that escape the detector could be replaced by a dark gauge boson, or more complicated dark sector final states, without affecting our main conclusions.
Common to all these possibilities is that in addition to the invisible final state particles (in our case φ 1 and ψ 1 ), there has to be at least one additional light mediator with a O(1 GeV)-mass in order to have large violations of the Grossman-Nir bound. In Models 1 and 2 the mediator is another scalar, φ 2 , while in Model 3 there are two mediators, the fermion ψ 2 and scalar φ. The scalar mediators mix with K L and π 0 , which then leads to enhanced Br(K L → π 0 + inv) rates. The required mixings are small, and thus for large parts of parameter space the most stringent constraints are due to the present KOTO upper bound on Br(K L → π 0 + inv). If the anomalous events seen by KOTO turn out to be a true signal of new physics, then these models are natural candidates for their explanation.
In Models 2 and 3 the lightest states, φ 1 and ψ 1 , can be dark matter candidates. For the restricted mass range m π ≤ m φ 1 ,ψ 1 ≤ (m K − m π )/2 and suitable parameter ranges these particles can be the thermal relic. For lighter φ 1 or ψ 1 new annihilation channels are required. For example the mediators could couple to either electrons or neutrinos, in addition to the couplings to quarks.
In the numerics we followed the principle of minimality and switched on the minimal set of couplings required for large violations of the GN bound. We took great care to ensure that the radiative corrections do not modify the assumed flavor structure and potentially invalidate our conclusions. In the future, it would be interesting to revisit our simplified models in more complete flavor models which fix all the couplings to quarks. An even more ambitious possible research direction could be to explore whether the light mediators could be tied to the SM flavor puzzle itself, e.g. along the lines of Ref. [46]. We leave this and related open questions for future investigations.
A The K → πψ 1ψ1 decay rate For completeness we give here the explicit expression for the K → πψ 1ψ1 differential rate in Model 3. The expressions become rather involved because of the presence of the fermions in the final states. For instance, the analytic expression for the rate can only be given as a double differential rate since both variables enter the loop diagram, cf. Fig. 18, in a non-trivial way. At the end of the appendix we also comment on how this decay defies the helicity formalism used for semileptonic and flavor changing neutral currents.
The conversion from a form M = LūP L v + RūP R v + L µū γ µ P L v + R µū γ µ P R v, used in  with X L ≡ Imĝ sd∆ f K f π − Imḡ sd 16π 2 F (φ) L (Ī 4 ) and X µ L = X L |Ī 4 ,∆→Ī µ 4 ,∆ µ , whereas for K + → π + ψ 1ψ1 (5.8) the decomposition reads V [A] µ + = B 2 0 /2 |y 21 | 2 ± |y 12 | 2 X µ + , (A. 8) with X + ≡ −ḡ sd 16π 2 F (φ) + (Ī 4 ) and X µ + = X + |Ī 4 →Ī µ 4 . It seems worthwhile to point out that this decay cannot be cast in the Jacob-Wick helicity formalism, since it does not corresponds to a chain of 1 → 2 decays. This also applies, e.g., to the generalisation of the formalism to effective theories used for B → K + − [48]. The issue is that in Model 3 the φ particle breaks factorization of the fermion part and the rest in the same way as the photon does between the lepton and the quarks, cf. Section 5.3. in Ref. [48]. On a technical level, this can easily be seen from the decomposition of the vector matrix element which necessitates all independent momenta of the decay. In the B → K + − case, induced by the standard dimension six effective Hamiltonian and no QED or electroweak corrections, the amplitude only depends on p π , p +p but not the difference p −p. However, using such a decomposition in the expressions given above does allow in practice for a fast numerical evaluation of the differential rate thereby retaining one of the main advantages of the helicity amplitude formalism.