An imprint of a new light particle at KOTO?

Recently, the KOTO experiment reported their new preliminary result of searching for the decay $K_L\to\pi^0\nu\bar{\nu}$. Three candidate events were observed in the signal region, which exceed significantly the expectation based on the standard model. On the other hand, the new NA62 and previous BNL-E949 experiments yielded a consistent result and confirmed the standard model prediction in the charged meson decay $K^+\to\pi^+\nu\bar{\nu}$. Furthermore, the two decays are bound by a well-motivated relation from analysis of isospin symmetry which is hard to break by new physics of heavy particles. In this work we study the issue by a systematic effective field theory approach with three simplest scenarios, in which the $K_L$ may decay into a new light neutral particle $X$, i.e., $K_L\to\pi^0X$, $K_L\to \gamma\gamma X$, or $K_L\to\pi^0XX$. We assess the feasibility of the scenarios by simulations and by incorporating constraints coming from NA62 and other relevant experiments. Our main conclusion is that only the scenario $K\to\pi XX$ for a long-lived light scalar $X$ has the potential to accommodate the three candidate events at KOTO and the NA62 result simultaneously while the region below the KOTO's blind box provides a good detection environment to search for all three scenarios for a relatively heavy $X$.


I. INTRODUCTION
The flavor changing neutral current decays of the neutral and charged kaons K → πνν provide a clean venue to examine precisely the standard model (SM) and to search for new physics beyond it. Recently, the KOTO experiment reported their preliminary result for the decay K L → π 0 νν using data collected during years 2016-2018 [1]. If the three candidate events observed in the signal region are confirmed in the future, it would imply a decay branching ratio [2]: B K L → π 0 νν KOTO = 2.1 +2.0(+4.1) −1.1(−1.7) × 10 −9 , (1) at 68% (95%) CL, corresponding to the KOTO single event sensitivity, SES π 0 νν = 6.9 × 10 −10 . This value is significantly larger than the SM prediction B K L → π 0 νν SM = (0.34 ± 0.06) × 10 −10 [3], which would pose a strong challenge to the SM.
This tension between the neutral and charged sectors is even exacerbated due to a theoretical relation between the two. Based on the well-motivated assumption that the decays are dominated by the interactions with isospin change ∆I = 1/2, they are related by the so-called Grossman-Nir (GN) bound [6], B(K L → π 0 νν) 4.3B(K + → π + νν).
Together with the bound in the charged sector of Eq. (4) this implies B K L → π 0 νν < 8.8 × 10 −10 at 95% CL, which would lead to at most 1.3 instead of 3 candidate events at KOTO. This bound is hard to break as it is immune to low energy effects at leading order of new heavy particles, and may only be violated by ∆I = 3/2 interactions of higher dimensional operators between quarks and neutrinos; for recent discussions and attempts, see for instance, Refs. [7][8][9][10]. The challenge here is to accommodate the KOTO result while respecting the measurements at NA62 and others. In this work we propose to relax the tension by assuming a light neutral particle that can appear in the final states of kaon decays. We will consider three scenarios, one of which has been suggested in the literature, and will assess by simulations whether they are practically feasible. Since π 0 decays into a pair of photons and the neutrino pair appears as missing energy, there are three simplest scenarios that could mimic the SM decay searched for at KOTO, namely, K L → π 0 X, K L → γγX, and K L → π 0 XX, in which X appears as missing energy.
We will first examine in Sec. II the most popular scenario K → πX [2,[11][12][13][14][15][16][17]. We will assess whether it is feasible by simulating its event distribution and comparing it with that measured at KOTO. We will also incorporate constraints coming from other measurements on K L and K + decays. In Sec. III we consider the loopinduced process of K L → γγX which could disguise itself as a signal event for the decay searched for at KOTO but has no counterpart at NA62. If the new particle couples only in pairs to quarks for one reason or another, it will also appear in pairs in the kaon decays. We will therefore make an analysis on the decays K → πXX in Sec. IV. In the course of our analysis we will make some suggestions that may deserve further study in the KOTO and NA62 experiments. We summarize our main results in the last section V. Our analysis in the following sections is based on the effective field theory approach. The effective interactions at leading order and the amplitudes and branching ratios for various decays are detailed in App. A, B, and C, and the simulation setup and validation is discussed in App. D.
The signal region for the K L → π 0 νν search at KOTO is also suitable for the two-body decay K L → π 0 X, where X acts as missing energy. Assuming parity and working to leading order in low energy expansion, X can be a scalar or vector particle. However, as shown in Fig. 1, the signal distributions in the Z vtx − P T plane of the two decays are different. Here Z vtx is the π 0 decay vertex position projected onto the K L beam direction, and P T is its transverse momentum with respect to the beam. We have performed the simulations by applying the kinematical cuts proposed in KOTO's 2015 data analysis [18] and assuming the signal region chosen in its 2016-2018 data analysis [1]. While the decay K L → π 0 νν has a relatively uniform distribution in the signal region, K L → π 0 X concentrates in a narrow interval of P T . Furthermore, the heavier the X is, the lower the maximal P T is. For a heavy X, it is hard to accommodate the candidate events with a high P T at KOTO.
The branching ratio for K L → π 0 X corresponding to a specific number of candidate events, N signal , that can be accommodated is estimated by the relation Here is the detection efficiency in the signal region at KOTO [1] upon imposing various kinematical cuts [18]. We have estimated the ratio of the two efficiencies by simulations and then used the above quoted SES π 0 νν (with an assumed relative error of 10%) to obtain SES π 0 X . The total background is 0.10±0.02 in the signal region (about one half from the decay K L → π 0 νν and another half from its SM background) and 0.08±0.06 in the region below the blind box [1]. The branching ratio for K L → π 0 X with two-sided 68% CL limits is shown in Fig. 2 as a function of m X in various mass intervals where a specific number N signal of candidate events are accommodated. Note that the choice of the interval delimiters is not sharp but only a rough estimate based on simulations. When X is so heavy (roughly m X > 190 MeV) that none of the candidate events can be accommodated or even all of its signals drop below the blind box (m X > 270 MeV), we then use the above background estimates to set a 90% CL upper bound on the decay branching ratio. Now we examine whether the decay K L → π 0 X offers a feasible interpretation to the candidate events observed in the KOTO signal region by a comprehensive analysis of the decay distribution and the limits set by other experiments.
We consider first the X mass intervals m X ∈ (100, 165) ∪ (260, 354) MeV. These intervals were not taken into account in the NA62 experiment since its sensitivity to the decay K + → π + X is considerably degraded by the large backgrounds K + → π + π 0 (γ) and K + → π + π 0 π 0 , π + π + π − , respectively. Thus the restrictive GN bound is practically not in action. Nevertheless, the second interval is obviously not supported by the KOTO signals as seen in Fig. 1 and the first one cannot provide a perfect solution either. For X in the first interval, i.e., m X ∼ m π 0 , the decay K L → π 0 X can only accommodate two candidate events but not the one with a high P T ∼ 238 MeV. If X is stable and invisible and if we leave aside the last event, we may obtain the branching ratio with two-sided limits at 68% (95%) CL at, e.g., m X = 135 MeV, For X of other mass, it must be so long-lived to be invisible at KOTO while short-lived to decay into SM particles to be vetoed at NA62, thus avoiding the constraint from the GN bound. In order to incorporate all three candidate events, X must be light enough, m X 60 MeV, but still it is difficult to offer a feasible solution. First of all, not all observed events are gracefully in the main distribution region as can be seen in Fig. 1. More importantly, as we will detail below, this light mass region is tightly constrained by other experiments.
For m X 60 MeV, the signals of K + → π + X at NA62 would be predominantly located in its signal region 1 for the K + → π + νν search [4] with missing mass squared m 2 miss ∈ (0, 0.01) GeV 2 , where no events were observed, implying at 90% (95%) CL. As we did not simulate the decays K + → π + X, π + νν at NA62, we obtained the upper limit by assuming an equal SES = (0.346 ± 0.017) × 10 −10 for both decays and by utilizing that the expected signals and background for K + → π + νν are respectively 0.82 ± 0.15 and 1.21 ± 0.14 corresponding to their total expectations in signal regions 1 and 2 at NA62 [4]. A similar limit, B (K + → π + X) E787 < 0.52 × 10 −10 (90% CL) for m X 115 MeV, was obtained by E787 [19]. Now we can employ the above results to obtain a strong constraint on the branching ratio of K L → π 0 X. We first recall that the real, physical branching ratio may differ its measured value in a specific detector (det) in an experiment such as KOTO or NA62 if X decays [2]: where τ X is the lifetime of X, and p/m X and L are the effective boost factor and detector size, respectively. We apply the above relation to both K L → π 0 X at KOTO and K + → π + X at NA62, and employ Eq. (5) for real branching ratios to link the two, so that we have where r = (L/p) NA62 /(L/p) KOTO ≈ 1.49, and for brevity B with a subscript indicates its measured or real value. Plugging in the central value B(K L ) KOTO ≈ 1.44 × 10 −9 from Fig. 1 and the upper bound on B(K + ) NA62 in Eq. (9), we arrive at the lower bound at 90% (95%) CL: This lower bound of order 10 −7 has been strictly excluded by other measurements. Being light, X can only decay to e + e − and/or γγ. For X → e + e − , B K L → π 0 X (suppressing the subscript) is tightly constrained by the measurements of B(K L → π 0 e + e − ), which is expected to be less than 2.8 × 10 −10 at 90% CL at the KTeV experiment [20]. For X → γγ, the limits mainly come from the measurements of B(K L → π 0 γγ) and its spectrum at KTeV [21] and NA48 [22]. The results B K L → π 0 γγ = (1.29 ± 0.03 ± 0.05) × 10 −6 at KTeV and B K L → π 0 γγ = (1.36 ± 0.03 ± 0.03) × 10 −6 at NA48 are consistent with the SM prediction ∼ 1 × 10 −6 [23][24][25][26]. Furthermore, both theoretical calculations and experimental observations give a consistent distribution which is dominated by the invariant mass interval m γγ ∈ (160, 360) MeV. If the above lower limit Eq. (12) held true for a light X, one would expect a peak or at least enhancement around m γγ = m X ∈ (0, 60) MeV in the experiments. But no such enhancement was observed at either experiment. Actually, NA48 reported no signal events for m X ∈ [0, 40] MeV, and set the limit B K L → π 0 X < 6.0 × 10 −9 (90% CL) for m X ∈ [30, 110] MeV [22]. These are great obstacles to the interpretation of the KOTO result in terms of the decay K L → π 0 X, especially when X is a light particle.

III. KL → γγX
In the KOTO experiment the directions of the photon pair were not recorded, and the reconstruction of the decay K L → π 0 νν was based on the assumption that the photon pair arises from the π 0 decay which in turn is a product of the K L decay in the beamline. This motivates us to consider the decay K L → γγX in which the photon pair would be misidentified as coming from the π 0 decay and X appears as missing energy. It is worth mentioning that this process is not constrained by the GN bound as its counterpart in the charged sector is absent.
The decay K L → γγX can appear at the one-loop order. The effective field theory calculation of the decay rate is delegated to App. B. Assuming parity symmetry, X could be a pseudoscalar or an axial-vector. In this section we analyze its phenomenological aspects by simulations. The distribution in the Z vtx − P T plane of the reconstructed events is depicted in Fig. 3 at three typical masses of the pseudoscalar X P . A similar distribution was found for the axial-vector case and thus not shown separately. As we can see from the figure, the KOTO's three candidate events could be covered over a range of m X , but the decay events more tend to be located below the signal region. If we naively insist that the KOTO's three events are faked by K L → γγX, more events should be found below the signal region, which however is not the case. Therefore, we will not pursue this idea further, but employ KOTO's results to work out an experimental upper limit on the decay K L → γγX instead.
We choose as our signal region for K L → γγX specified by Z vtx ∈ [2.9, 5.1] m and P T ∈ [0, 120] MeV, which excludes KOTO's three events but includes most of our decay events with a background of 0.08 ± 0.06. According to the observed zero event in the region [1], we obtain an upper limit on B(K L → γγX P ) and relevant Wilson coefficients as a function of m X P in the top and bottom panels of Fig. 4 respectively. Also included is the projected future detection capability of KOTO (dashed curves) when an SES π 0 νν = 3.0×10 −11 is reached. In the bottom panel, we also display limits coming from other experiments: K + → π + π 0 νν at E787 [27], K L → π 0 π 0 νν at E391a [28], and the K 0 −K 0 mixing [29]. For a light X P , the limit is dominated by K L → π 0 π 0 X P , while for a heavy X P the strongest limit comes from the K 0 −K 0 mixing; in between (roughly for m X P ∈ (220, 350) MeV) the KOTO can yield the best upper bound in the future. Similar limits for an axial-vector X A are shown in Fig. 5.
Finally, we discuss the case in which a pair of X particles appear in the K meson decays for which the kinematics will be very different from that of K → πX that we studied in Sec. II. For a fermionic X this has been suggested in Ref. [30]. We consider here a bosonic X which for one reason or another only couples in pairs to the quarks and may be a scalar or vector of either parity. Our simulation results for the signal distribution in the Z vtx − P T plane are shown in Fig. 6 at three different masses. We find that the events for a light X are almost evenly distributed in the signal region as in the case for K L → π 0 νν but differ significantly from the decay K L → π 0 X as shown in Fig. 1. As X becomes heavier, the maximal P T also drops. But the truncation of the distribution around the maximal P T is not as sharp as in the two-body decay, and this leaves some space to explain the KOTO signal with a high P T . For even larger masses, the distribution will shift below the blind box.
The branching ratio for K L → π 0 XX from above simulation is shown in Fig. 7. The KOTO candidate signals prefer a light scalar X S with m X S 30 MeV, which must be unstable and have an appropriate lifetime to avoid the GN bound. The requirement is that with a high probability the two X S s should be both invisible to KOTO while at least one of them decays into SM particles (γγ or e + e − ) which can be vetoed at NA62. Qualitatively speaking, in terms of kinematics, the experimental constraints on the πX S X S mode are looser than those on the πX S mode. But to assess the feasibility of the scenario, we have to determine the appropriate lifetime and mass of X S . The physical branching ratio for K → πX S X S at KOTO or NA62 can be expressed as where f represents the probability that the two X S s propagate respectively at the momentum p 1 , p 2 for a distance L 1 , L 2 to exit the detector without decay. The integration over the signal region ensures that all situations of the K decay have been taken into account. However, it is difficult for us to manage a systematic simulation to cover all information due to the complexity of displacedvertex simulation. Considering that the two X S s have a large boost at both KOTO and NA62, they should have a high probability of flying in the beam direction with a small angle to the veto plates. For further estimation, Bottom: upper limits on Wilson coefficients: |C P sd | from K + → π + π 0 + invisible by E787 experiment (red solid), |Re(C P sd )| from K 0 −K 0 mixing (green solid), K + → π + π 0 + invisible by E391a experiment (blue solid), KL → γγXP by KOTO's current result (purple solid) and future expectation (orange dashed).
we found that p 1 = p 2 ≡ E and L 1 = L 2 ≡ L serve as a good approximation, which leads to the simplification: The two-dimensional probability function F (L, E) is much simpler compared to the four-dimensional one. For K L → π 0 X S X S at KOTO, it can be extracted from our current simulation; for K + → π + X S X S at NA62, we have adopted a rough simulation in which only some cuts on the signal region are considered.
As an example of simulation we consider the scenario that a scalar X S of mass m X S = 10 MeV decays into a photon pair. The branching ratio measured at KOTO can be read off in Fig. 7, with the central value being B(K L → π 0 X S X S ) KOTO = 2.7 × 10 −9 . For NA62, assuming that K + → π + X S X S has the same acceptance as K + → π + νν, we obtain B(K + → π + X S X S ) NA62 < 1.60 × 10 −10 at 95% CL. By incorporating all these into Eq. (14), we obtain B(K L → π 0 X S X S ) real as a function of the lifetime τ X S shown in Fig. 8 as the green curve, and an upper limit on it (red curve) from NA62 with the aid of the GN bound. This yields the allowed region with τ X S 2.2 × 10 −7 s, within which the veto information can also be gained from the figure: at least 74.6% of the K + → π + X S X S signals are vetoed at NA62, while for τ X S 10 −8 s the two X S s at KOTO nearly do not decay before exiting the detector. On the other hand, the untagged K L branching ratio [29] constrains B(K L → π 0 X S X S ) < 1%, which leads to τ X S 1.7 × 10 −10 s. But we emphasize once again that this is only a rough estimate, and an appropriate deter- mination of the lifetime and mass of X S could only be achieved by a systematic detector simulation. We advocate that the KOTO and NA62 collaborations will take this endeavor in their future experimental analysis.

V. CONCLUSION
We have investigated by simulations the feasibility to interpret the recent KOTO result in terms of a new neutral particle that appears in the kaon decays. Since π 0 decays into a pair of photons and neutrinos appear as missing energy, we have considered three scenarios, i.e., K L → π 0 X, γγX, π 0 XX. Our results can be summarized as follows. The simplest scenario K L → π 0 X is difficult to accommodate all three candidate events at KOTO, especially the one with a high P T . The signal events for a relatively heavy X are mainly distributed below the KOTO's signal region, which have been employed to work out a bound on the decay branching ratio. While the KOTO result favors a light X, our comprehensive analysis on the NA62 and other experiments sets a strong constraint that essentially excludes this potential. Since the scenario K L → γγX has no constraint at NA62, it is free of the GN bound. While the three candidate events can be accommodated, the distributions do not fit: more events would be expected below the blind box. We have used the latter to set a constraint on the branching ratio, and compared it with those from other measurements and the expected KOTO's future capability. In terms of the signal distribution, the KOTO candidate events favor the third scenario K L → π 0 X S X S , where X S is a scalar with m X S 30 MeV. To accommodate the measurements at NA62, X S should be unstable but long-lived, whose lifetime is estimated to be 1.7 × 10 −10 s τ X S 2.2 × 10 −7 s at m X S = 10 MeV. But a more precise result would necessitate sophisticated simulations which we hope the KOTO and NA62 collaborations will perform in their future experimental analysis. Appendix A: Effective field theory framework For our purpose of accounting for the KOTO anomaly we assume a new real neutral particle X of mass below a few hundreds MeV. It may be a scalar (X S ), pseudoscalar (X P ), vector (X V ) or axial-vector (X A ) particle as appropriate to the scenario under consideration. We start with the low energy effective field theory that contains the X particle in addition to the light quarks and leptons and has the QCD and QED gauge symmetries. As we are mainly concerned with the transitions between the down and strange quarks, we only consider their couplings to the X field: where p, r refer to the d, s quarks and the Wilson coefficients are 2 × 2 Hermitian matrices in the d, s space. When studying the third scenario K → πXX, we switch off the above single-X couplings but switch on the following double-X couplings: where X and X µ may have any parity without affecting our discussions in the work. Below the chiral symmetry breaking scale the Nambu-Goldstone bosons (NGBs) become the dynamical degrees of freedom, whose low energy physics is determined at leading order by [31,32] where U = exp i √ 2Φ/F 0 exponentiates the octet NGBs and As usual, F 0 is the decay constant in the chiral limit and B parameterizes the quark condensate qq = −3BF 2 0 . The new field X gets involved in the form of additional terms in the external sources to QCD, which are Hermitian matrices in the light flavor space: in the case of single-X couplings in Eq. (A1), and in the case of double-X couplings in Eq. (A2). Note that the QED interaction is contained as usual in the sources l µ , r µ , which we do not write explicitly. The effective interactions contained in Eq. (A3) will be applied to calculate the quantities in the following appendices.

Appendix B: Decay amplitudes and distributions
The decay width for K L → π 0 X is, for a scalar X S , where p π,z is the π 0 momentum component in the K L beam direction, or, for a vector X V , The decay K L (k) → γ(q 1 )γ(q 2 )X(p) involves only neutral particles, and can only take place at the one-loop order whose Feynman diagrams are shown in Fig. 9. The result is finite, and the (spin-summed) squared matrix element is, for a pseudoscalar X P , where r π,K = s/(4m 2 π ± ,K ± ), s = (q 1 + q 2 ) 2 , and α ≈ 1/137 is the fine structure constant, or for an axial-vector X A , The one-loop function is For simulations, we use the distribution where t = (k − q 1 ) 2 . The decay width for K L (k) → π 0 (p)X(q 1 )X(q 2 ) is easily computed to be, for a scalar X S or vector X V respectively, where r X = s/(4m 2 X ), s = (q 1 + q 2 ) 2 , and t = (k − q 1 ) 2 .

Appendix C: Other experimental constraints
In this appendix we list other experimental constraints that have been used in Sec. III for a comprehensive analysis. The experimental upper limits on the four-body kaon decays are [27,28]: The corresponding amplitudes for the decay K(k) → π(p 1 )π(p 2 )X(q) read, for a pseudoscalar X P , and for an axial vector X A , where s = (p 1 + p 2 ) 2 , t = (k − p 1 ) 2 , u = (k − p 2 ) 2 , and is the polarization vector of X A . Moreover, the K 0 −K 0 mixing can give limits on the couplings. The experimental measured quantities are the K L − K S mass difference ∆M K and the CP violation parameter K [33][34][35], whose current experimental values are [29], ∆M K = (3.484 ± 0.006) × 10 −12 MeV, (C7) | K | = (2.228 ± 0.011) × 10 −3 . (C8) Considering the theoretical uncertainties from longdistance contributions in ∆M K , we require the new contribution do not exceed the experimental value. Relatively, the calculation of K is more credible, hence we require the new contribution to be less than 30% of its experimental value. The limits on the Wilson coefficients then read, for a pseudoscalar X P , , and for an axial-vector X A , In this appendix we will briefly describe how we do the simulations for various K L decays. A systematic simulation is complicated and time consuming, so we content ourselves with a simplified version of it in this work. We find that even in this simple framework we can get accurate results as in a systematic simulation. In the following, we will explain our procedure and compare our result on K L → π 0 X with KOTO's to verify our simulation.
We first generate initial K L particles according to the momentum distribution of K L measured experimentally [36], which will have a certain probability of decay in the detector. All of the K L decay modes in this paper produce two photons, and the distributions of the energies and positions of the photons are largely dependent on the probability distribution functions. For the two-body decay K L → π 0 X, we use a uniform distribution function as in Eqs. (B1) and (B2). For the three-body decay K L → π 0 νν, we adopt the same distribution function as in Eq. (S1) in the supplemental material of Ref. [2]. The distribution functions of K L → Xγγ and K L → π 0 XX are determined by Eq. (B6), and Eqs. (B7) and (B8) respectively. It is worth mentioning that for the decay modes K L → π 0 νν, K L → π 0 X, and K L → π 0 XX, the two photons are originated from the π 0 , while for K L → Xγγ, the two photons are directly generated by K L ; only the former decays will lead to an invariant mass of the two photons m γγ m π 0 . Then the photons are captured by the CsI calorimeter in the detector, and we record their energies and positions in our simulation. In order to better simulate the detector's response to the photons, we include the energy and position resolution of the CsI calorimeter, which can be found in Ref. [37].
By assuming that the two photons produced on the beam axis are from π 0 decay (in the interesting case of K L → Xγγ the two photons are also required to have an invariant mass m γγ m π 0 in order to fake the signals), we can reconstruct the decay location Z vtx and π 0 's transverse momentum P T by combining the information of photons' energies and positions. Then we use the same selection criteria as KOTO [18] to analyze the reconstructed events. We use the signal region in KOTO's analysis of 2015 data [18] for validation, whose result is shown in Fig. 10, but employ the new signal region [1] to analyze our decay modes. It should be noted that we do not consider shape-related cuts and veto cuts, whose efficiencies are considered the same for different processes and can be estimated as r = KOTO π 0 νν / π 0 νν , where π 0 νν is from our rough simulation and KOTO π 0 νν is the real acceptance from KOTO's systematic simulation. We calculate the real acceptance of each process by multiplying our result with r. After applying the reconstruction and various kinematical cuts, we finally obtain the distributions of signal events in the Z vtx −P T plane, and the branching ratios and Wilson coefficients of some processes can be evaluated accordingly.