KOTO vs. NA62 Dark Scalar Searches

The two kaon factories, KOTO and NA62, are at the cutting edge of the intensity frontier, with an unprecedented numbers of long lived and charged Kaons, ~ 10^{13}, being measured and analyzed. These experiments have currently a unique opportunity to search for dark sectors. In this paper, we demonstrate that searches done at KOTO and NA62 are complementary, both probing uncharted territories. We consider two qualitatively different physics cases. In the first, we analyze models of axion-like-particles (ALP) which couple to gluons or electroweak gauge bosons. In the second, we introduce a model based on an approximate strange flavor symmetry that leads to a strong violation of the Grossman-Nir bound. For the first scenario, we design a new search strategy for the KOTO experiment, K_L ->pi^0 a ->4 gamma . Its expected sensitivity on the branching ratio is at the level of 10^{-9}. This demonstrates the great potential of KOTO as a discovery machine. In addition, we revisit other bounds on ALPs from Kaon factories, highlighting the main sources of theoretical uncertainty, and collider experiments, and show new projections. For the second scenario, we show that the model may be compatible with the preliminary analysis of the KOTO-data that shows a hint for New Physics.


Introduction
The Standard Model (SM) of particle physics is a successful description of Nature, especially given the discovery of the Higgs boson at the LHC [1,2]. The SM describes forms of matter which interact via the electro-magnetic, weak and strong forces. However, the SM is incomplete as it can not account for e.g. the observed baryon asymmetry of the universe, neutrino masses and mixings, and the origin of Dark Matter (DM). Motivated by the finetuning problem of the electroweak (EW) scale that conventionally requires TeV new physics (NP) which also characterizes the DM sector, tremendous efforts have been made to search for new states at the energy frontier, and yet, so far there is no conclusive sign of the beyond the SM (BSM) physics. On the other hand, a NP sign could appear as a light weakly coupled state, for instance associated with a pseudo Nambu Goldstone boson (pNGB) field, and the representative example is an axion or axion-like-particle (ALP) 1 . The mass scale of the pNGB can be substantially lighter than the GeV scale, and its interaction strength with SM particles can be suppressed by a higher symmetry-breaking scale. This type of particle can be tested at high-intensity experiments, such as in rare meson decay measurements, at B-factories, beam-damp experiments, and neutrino experiments.
In this paper, we demonstrate that both Kaon factories have great opportunities as discovery machines of light new particle. In Sec. 2, we introduce two qualitatively different physics cases that show that NP searches done via charged Kaon decays at NA62 and via neutral Kaon decays at KOTO are complimentary, as opposed to be strongly linked with each other. This is in contrast to what one would naively expect to be the case due to the GN bound.
First, we consider an ALP (a) with coupling to gluons or W bosons as a representative candidate of pNGB. In this context, we propose a novel search for the KOTO experiment (Sec. 4). Specifically, KOTO can search for K L → π 0 a where the ALP decays to di-photon. This search will be complementary to the K L → π 0 a with an invisible ALP search that is already performed by the collaboration. These two channels together will probe experimentally unexplored parameter space of the ALP coupled to SU(2) gauge bosons (Sec. 5.1) or to gluons (Sec. 5.2), in the mass range from 10 MeV to 350 MeV. NA62 will also probe parameter space through the corresponding K + → π + a decay that we analyze.
The second scenario we analyze in this paper is a theory with an approximate strange flavor symmetry, with an additional light, flavon-like, complex scalar field, φ. We discuss its phenomenology in Sec. 6. The flavor preserving coupling allow for a SM singlet final state, consisting of the real (σ) and imaginary (χ) parts of φ, to be accessible only to K L and not to its charged isospin-partner. Therefore, breaking the GN relation. KOTO is particularly sensitive to such a scenario, once we allow the χ to decay to two photons. Furthermore the expected signal can be made compatible with the preliminary analysis of the KOTO-data that shows a hint for NP [17] (though more investigation of the collaboration is needed). Other explanations of this anomaly can be found in [6][7][8][9][10][11][18][19][20][21][22]].

Light Scalars at Kaon Factories
Here, we describe the two new physics scenarios where the Kaon factories can play a major role probing the parameter space.

Massive axions, ALPs and pNGBs
The Goldstone theorem provides one of the most compelling motivation for the presence of light scalars as their masses are protected by a shift symmetry. The simplest manifestation of the Goldstone theorem is the case of a spontaneously broken U(1) symmetry that leads to the presence of a light ALP. Such a state can be motivated by a solution of the hierarchy problem [23], the strong CP problem [24][25][26], the flavor puzzle [27], and combinations of these with DM physics [28][29][30][31][32][33][34]. For concreteness, to motivate our scenario we focus on the QCD axion case, however, the essence of our reasonings below holds for a broader class of ALP models. The typical breaking scale of the Peccei-Quinn (PQ) symmetry [24], F a , considered in literature is rather high. The standard axion window is 10 9 F a 10 12 GeV.The upper bound is due to the over-production of axion as dark matter, and the lower bound comes from astrophysical observations [35].
However, there is a theoretical concern about the quality of the PQ symmetry with a high decay constant [36][37][38][39][40]. Any global symmetry is believed to be broken by the UV physics: the quantum gravity does not respect global symmetries; or any global symmetry can be an accidental symmetry of the UV physics. In the effective field theory, this conjecture implies that higher dimensional operators suppressed by the UV physics scale, Φ|Φ| D−1 /Λ D−4 U V , can explicitly break the PQ symmetry where Φ is a field which carries a non-zero PQ charge and has a VEV of F a . These operators ruin the PQ mechanism because this operator shifts the minimum of the axion potential away fromθ = 0, where ∆ is a non-aligned CP phase that is generically expected to be of order one. Even though the deviation is suppressed by a high scale Λ U V ≤ M pl , the effect in theθ can be significant because of two factors: (1) the original axion potential is not very steep, m 2 a F 2 a ≈ m 2 π F 2 π ; (2) the precision of the neutron EDM measurement is accurate, δθ 10 −10 . Therefore, operators up to D 10 need to be absent to maintain the PQ mechanism. This situation is unsatisfactory from the low energy point of view. Some mechanism should maintain the quality of the global PQ symmetry to be extremely good to solve the strong CP problem. This problem is not unique to the QCD axion but is also common to other solution to the QCD CP problem [41] and other mechanisms that strongly rely on precise global symmetries [42][43][44].

Heavy Axion as a Consequence of the Quality Problem
To construct theories that are protected against Planck suppressed operators of D ≥ 5, the favored decay constant is necessarily low. Assuming the standard relation of axion mass and decay constant, m a ≈ m π F π /F a , and requiring a small deviation, δθ < 10 −10 , one can obtain the bound on the effective decay constant and the mass, This parameter space is similar to the original Weinberg-Wilczek axion [25,26]. This motives us to search for axions with a mass at and above the MeV scale. The low decay constant along the standard QCD axion relation has been excluded by astrophysical observations and beam-dump experiments. However, the bounds do not apply if there is an additional contribution to the axion mass. Many phenomenological studies for ALPs show that parameter space with very low decay constant are poorly constrained if the mass is heavier than ∼ 50 MeV (see recent works, for example, Refs. [45][46][47][48][49][50][51][52]). Indeed, there are models of heavy axions without the standard axion relation, where the strong CP problem can be addressed [42,51,[53][54][55][56][57]. As an example of non-minimal heavy axion, we consider a scenario with a mirror strong sector. The mirror sector shares a strong CP phase and quark phases with the SM ensured by a Z 2 symmetry, and a single axion relaxes the two CP phases of the SM and mirror sectors. With soft breaking of the Z 2 symmetry, a higher confinement scale in the mirror sector is achieved, which is an extra source of the axion mass.
We now revisit the quality problem and the bound on the axion mass based on the scenario with the mirror strong sector. First, the higher dimensional operator should be sufficiently suppressed as in Eq. (2.2), The axion mass is dominated by the contribution from the mirror sector because its confinement scale Λ is much higher than the SM one, Λ QCD , where m q are SM quark masses. Generally, we expect a hierarchy Λ F a because Λ is generated by dimensional transmutation, but the confinement scale can be up to Λ = F a , and consequently m 2 a < m q F a . Combining this with Eq. (2.4), we get .
(2.6) Specifically, the bounds for D = 5 case are 1 MeV m a 5 GeV m q 10 MeV (2.8) The above parameter space is only weakly covered by the current experimental probes. Since part of this mass range is within the range of Kaon experiments (particularly the lower mass range), it is very important to develop a search program to discover heavy ALPs at Kaon factories. For the phenomenological study of heavy axions, we consider two simplified models: the first involves a ALP coupled to the electroweak sector of the SM, and the second a ALP coupled to gluons (for more information see Sec. 5).

The generalized GN bound and how to avoid it
Under fairly general assumptions, the K L → π 0 νν rate can be strongly constrained by the K + → π 0 νν rate via the Grossman-Nir (GN) bound [5]: The numerical factor comes from the difference in the total decay widths of K L and K + , isospin breaking effects, and QED radiative corrections [4,58]. The GN bound only relies on the following assumptions [5]: First, the isospin symmetry, which relates the decay amplitudes of K ± to the ones of K 0 andK 0 . Second, the ratio of the K andK 0 decay amplitudes to the corresponding sum of final states is close to unity, where if the final state is CP eigenstate it means no CPV in the decay. For the πνν final state, within the SM, it is expected to be an excellent approximation. The above assumptions are not easy to be violated even when going beyond the SM. Inspired by [59], we shall construct a model based on an approximate global flavor symmetry, that avoids the GN bound via exploiting strong isospin breaking (see [12][13][14][15][16] for relevant discussions). To realize the idea, we add a light complex scalar, φ, which carries a half strange (or second generation doublet) flavor charge. This implies that we expect the following operator to be allowed by the symmetry and present in the effective theory, in the down quark mass basis, y 1 HQ 1 sφ 2 /Λ 2 and/or y 2 HQ 2 dφ 2 /Λ 2 + h.c. , (2.10) where the first (second) operator corresponds to φ 2 carries a units (Q 2 ) flavor charge, and we assume φ = 0. In the broken electroweak phase, this effective Lagrangian leads to an effective operator y 1,2s dφ 2 + h.c. that induces the K L → σχ decay, with σ = Re(φ)/ √ 2 and χ = Im(φ)/ √ 2 (here, for simplicity, we assume an approximate CP conservation in the decay). Using NDA, from Eq. (2.10) we expect However, due to conservation of charge there is no analogous 2-body decay of the charged Kaon unless additional charge pions are added to the final state. This implies that the charged Kaon decay is suppressed, by two-vs-three-body (and possibly kinematical) phase space factors which implies a strong violation of the effective new physics GN bound. As discussed in Sec. 6, we find that the NP charged Kaon decays are suppressed by at least two orders of magnitude relative to the K L one. Thus, in such a scenario, it is possible that while, at present, the KOTO detector is sensitive to a NP signal, the NA62 one is not. The model, as presented above, has an exact φ-parity symmetry which renders the φ state stable. To achieve a visible signal at Kaon experiments, we add a CP conserving coupling, that is responsible to the decay of χ into two photons. Up to small symmetry breaking effects, to be discussed below in Sec. 6, σ would be stable and hence the final state of the K L → σχ(γγ) is similar to the K L → π 0 νν, which KOTO is searching for.

Overview
KOTO is an experiment searching for the rare neutral Kaon decay, K L → π 0 νν, whose branching ratio is expected to be (3.0 ± 0.3) × 10 −11 [3,4]. In the past, the E391a experiment, at KEK, set the most stringent limit on the branching ratio at 2.6 × 10 −8 [60]. The first KOTO analysis based on data collected in 2015 was able to set a bound at BR(K L → π 0 νν) KOTO < 3 × 10 −9 [61]. This is relatively close to the bound obtained from the charged decay, K + → π + νν, using the Grossman-Nir bound: KOTO is a fixed target experiment that utilizes a 30 GeV proton high intensity beam extracted from the J-PARC main ring accelerator. The produced Kaons are purified by a 20m-long beam line and enter in the detector of Fig. 1, as indicated by the arrow, where the beam axis is denoted as the Z direction. The flux of Kaons was measured by an engineering run in 2015 at Z ∼ −1.5m [62]. The actual detector consists of a CsI calorimeter (Ecal) at the front target and various veto detectors for charged particles and photons. The measured momentum distribution of the incoming K L flux is shown in black in Fig. 2 and it peaks at around 1.5 GeV. Then the Kaons decay in the decay volume at 2 m < Z < 6.148 m to produce pions or neutrinos, and the momentum distribution of the decayed K L is shifted towards lower values as shown by the orange histogram in Fig. 2. The neutral pions are reconstructed through the identification of photons that hit the CsI calorimeter with E γ > 2 MeV.
The output from the detector is the position of the photon energy deposition on the ECAL (on the plane perpendicular to the beam direction), and the timing of the hits. What is known is the energy of photons rather than their four-momenta because the decay vertex of the Kaon (effectively same as the pion), Z vtx , is unknown and the ECAL can measure only the photon energy. Furthermore, the final states of interest include no charged particles, which could provide directional information. In order to reconstruct the decay vertex, the standard technique is to impose at least one additional assumption regarding the invariant mass of the parent particle or intermediate particles [62,64] (see appendix A for a brief review). This procedure still has multi-fold ambiguities, but the correct vertex can be picked, at least based on statistical merit, by requiring that the reconstructed event describes the physical process. This challenge holds for SM processes such as K L → π 0 (γγ)νν, K L → π 0 π 0 → 4γ, as well as possible new physics processes, to be discussed below.

Flux, signals, and future plans
In Tab. 1, we summarize the most important numbers that characterize the intensity of the KOTO experiment. We compare the amount of data collected in 2015, to the one collected in 2016-2018. We then report the amount of data that is aimed to be collected in the coming years to reach the measurement of several SM events for K L → π 0 νν. The K L flux at the beam exist is usually reported. We refer to it as N 0 K L , which is calculated by protons on target (POT) [ K L → πνν within 3 m < Z < 4.7 m 1.48 × 10 11 2.3 × 10 11 ∼ 6.6 × 10 12 S.E.S of K L → πνν 1.3 × 10 −9 6.9 × 10 −10 ∼ 2.9 × 10 −11  Since most of the entering Kaons do not decay, we need to translate N 0 K L to the fraction of decays inside the relevant detector volume. The probability for the K L to decay inside the entire detector region, 2 m < Z < 6.148 m, calculated at truth level is 7.9%. The fraction of K L → π 0 νν relevant to the actual analysis, within the region of 3 m < Z < 4.7 m, using our reconstruction-level simulation was found to be 3.2% (consistent with [65]), while the one associated with multiple pion final state analysis, which corresponds to the region of 2 m < Z < 5.4 m was found to be 6.5%, at the truth level (consistent with [62]).
For the following discussion, it is useful to examine the search for the K L → π 0 νν decay as done by the KOTO experiment in more detail. The relation between the flux, acceptance and the S.E.S is given by The event selections are given in [61], and most of them are included in our analysis, except the veto and shower-shape cuts. To keep these cuts into account, we choose uniform efficiencies: for the veto cut we use 0.17, and for the shower-shape cut we use 0.52 [61]. The signal region after these cuts is defined by the transverse momentum of the reconstructed π 0 (γγ) and its decay vertex, Z vtx , that is, p min T (Z vtx ) < p π 0 T < 250 MeV and 3 m < Z vtx ]. Note that the vertex reconstruction in the K L → π 0 νν search occasionally has two-fold physical solutions for Z vtx , so we choose the one further from the ECAL. We have checked that this is the physical one in most cases. Even if we discards the events with the two-fold solutions, the acceptance changes by only O(1%). In Sec. 6, we will discuss how the requirement to have photons in the signal region affects the acceptance of a NP model that leads to a K L → σχ, χ → γγ decay.

New ALP searches at KOTO
The KOTO experiment can look for heavy axions or ALPs produced from K L decays. Particularly, as we will argue, searches could be designed to identify the decay topology K L → π 0 a, a → γγ, that we will analyze in detail in this section.

Reconstruction of the four photon signature
The signal of our interest is four photon final state from K L → π 0 (γγ)a(γγ). As already mentioned, in order to reconstruct the decay vertex, the standard technique is to impose at least one assumption on the invariant mass of the parent particle or intermediate particles.
This procedure still has multi-fold ambiguities, but the correct vertex can be picked within an error of a few percent by consistency checks. However, note that this technique works only for the anticipated decay topologies.

Axion reconstruction
Extending the standard technique, we propose a new algorithm to reconstruct the signal process K L → π 0 a → 4γ without knowing the axion mass, as follows.
1. Divide the four photons into two pairs which can be associated to mother particles, say the pair A consists of γ 1 and γ 2 and the other pair, B, consists of γ 3 and γ 4 . There are six such combinations of photon-pairs.
2. Obtain a candidate vertex, Z vtx , assuming that the mother of the photon-pair B is a neutral pion. This assumption holds for the signal as well as physical background processes K L → π 0 π 0 , π 0 γγ. Consequently, the invariant mass of the pair B can be written as, There are up to two solutions for Z vtx . One solution is often unphysical, since the reconstructed vertex is outside from the decay volume, or it is an imaginary number, and thus discarded.
3. Repeat the above steps for all the possible parings, which lead to twelve-fold ambiguities in a four-photon event.
4. If the pairing and the reconstructed vertex are found to be physically consistent, the four photon invariant mass, m 4γ , is also required to be peaked around the K L mass, and the combination that minimizes the resulting |m K L − m 4γ | is selected.
This reconstruction algorithm works quite well, allowing us to select the correct pair in the signal process. The fraction of the (preselected candidate) events that are correctly selected is 90% (70% for m a m π 0 ). The di-photon invariant mass of one pair is expected to have a peak at around the ALP mass, while the background events K L → π 0 π 0 have a peak in the same variable around the neutral pion mass.

Simulation
We develop a MC simulation based on our reconstruction algorithm, and estimate the acceptance of the signal and the K L → π 0 π 0 , K L → π 0 γγ backgrounds. We start from the known K L flux, and then let K L decay to π 0 π 0 , π 0 γγ, π 0 a, and subsequently decay the π 0 and a to γγ. For the photon energy measurements, we include the dominant smearing effects due to the ECAL.

K L momentum and vertex reconstruction
Our simulation aims at obtaining the two-dimensional distribution of the Kaons in terms of the reconstructed momentum and the reconstructed decay point, Z vtx . We generate the Kaon momentum according to its measured distribution at the beam exit Z −1.5 m, shown in black in Fig. 2. The Kaons then decay according to their lifetime of approximately 0.51 ns. We assume that the Kaons are fully aligned with the beam axis. We collect the decayed K L within the decay volume of the detector, 2 m < Z K L < 6.1 m where Z K L is the actual point of K L decay, and the decay probability is about 7.9% 3 . As a cross-check, we calculate the decay probability in the fiducial volume of K L → πνν analysis, that is, 3.2%, and it is consistent with the reported probability in [65].
Starting from K L momentum and decay vertex, {p K L , Z K L }, the flow of our simulation for the two body decays is is as follows, In the last step, the information of the photon momenta and the K L decay position is mapped onto the photon positions (x, y) on the plane of the ECAL. The three body decay K L → π 0 γγ is treated in a similar fashion.

Detector's finite resolution
To take into account the detector effects, that is the photon's finite energy and positionresolution, we smear each photon's energy and position following the detector resolution [68], 4 Taking these as a standard deviation of a gaussian distribution, we smear each photon hit as (4.7) Thus, outputs of our MC samples are (x, y, E) detect γ 1,2,3,4 , where the energy smearing dominates the total smearing. There are other detector effects, such as photon inefficiency, shower-shape, and timing-resolution, but these are beyond our simulation setup and the effects are expected to be minor because we reproduced shapes and normalizations of several measurements (see Appendix B).

Preselection of four photon events
Our preselection of four-photon events is similar to the one used for the four-photon analysis for the K L → π 0 π 0 decay [62]. We employ a series of basic cuts on photon energies and positions.
1. The four photons should hit the front ECAL, which is a circle of 1 m radius. No photons hit the main barrel of the detector (see Fig. 1).

As the CsI calorimeter has
Moliere radius of about 3.5 cm for the electromagnetic shower, the four photons are required to be inside a 90 cm radius, 3. The position of the innermost photon should be outside of the beam hole, max[|x i |, |y i |] ≥ 7.4 cm.
4. The photons should be well separated such that d min = min | r i − r j | ≥ 15 cm. This ensures that there are at least four clusters of hits, and, thus, events with four or more photons.

The minimal energy of each single photon should be min[E
6. The total photon energy should be i E γ i ≥ 350 MeV.
The efficiency of the above preselection cuts is about 7% for both signal and background, except the efficiency for K L → π 0 a with m a ≤ 20 MeV that is about 1% or less due to CUT4 and CUT5 described above.

Cuts after reconstruction
After the preselection defined in Sec. 4.3.1, the search for axion decay to a pair of photons, within the multi photon events proceeds via the following set of cuts 5 : 7. The reconstructed vertex should be within 2 m < Z vtx < 5.4 m , which defines the fiducial volume of this analysis.
8. The four-photon invariant mass should match the K L one, namely, |m 4γ − m K L | < 20 MeV.
9. The invariant mass of the photon pair which corresponds to the non-pion candidate, m γ 1 γ 2 , is required to be away from the neutral pion mass, |m γ 1 γ 2 − m π 0 | > 10 MeV. This cut particularly removes most of the K L → π 0 π 0 background. 10. To further remove the π 0 π 0 background, we examine all the possible di-photon pairings to check if any of them reproduces the K L → π 0 π 0 decay topology. The event is discarded if, for any pair assignment, the pair A satisfies |m γ 1 γ 2 − m π 0 | < 20 MeV and |m 4γ − m K L | < 50 MeV . Only 1.6% of π 0 π 0 background remains after this cut while the other decay topologies are almost unchanged.
When we compute the sensitivity to axion masses near the pion mass, then CUT10 is excluded. This cut would, in fact, substantially reduce the signal for m a = (130 − 140) MeV, The overall efficiencies of CUT1 -CUT10 are 9 × 10 −5 and 5% for K L → π 0 π 0 and π 0 γγ, respectively. The signal efficiency is 3-6% except for m a 20 MeV or m a ∼ m π 0 . Our background π 0 π 0 MC statistics is poor for m γγ < 50 MeV, we thus treat this region as a single bin for m a = 1, 10, 20, 30, 40 MeV. are also shown. We use a larger bin size for the π 0 π 0 background except near the pion mass due to the poor MC statistics of the remaining events. For this figure, we also include the efficiencies of veto and shower cut ( veto = 17% and shower = 52%).
So far, our simulation setup does not incorporate the veto cuts and the shower shape cut that are adopted by the KOTO analysis [61]. To take into account the efficiencies of the veto and shower shape cuts, we multiply the kinematic acceptance obtained above by veto = 17% and shower = 52% [61] for both signal and background 6 .
After all cuts and efficiencies, the remaining background events are mainly from combinatorics of K L → π 0 π 0 and an irreducible K L → π 0 γγ. The corresponding distributions are shown Fig. 3 in gray and red, respectively. In the same plot we also show the expected signal for axion masses of 50 MeV and 200 MeV respectively and BR(K L → π 0 a) = 10 −9 . The number of events is computed assuming the future KOTO luminosity of N 0 K L = 2 × 10 14 . Next we estimate the typical size of the ALP di-photon peak. Given the above cuts, we have found the RMS of the peak as a function of the mass by fitting the result of our MC simulation to an approximate functional dependence. By fitting the typical resulting width of the peak around m a to its RMS value, the peak region is defined as Fig. 4).
���� ���������� Before ending this section, let us briefly discuss another potential source of background: the three-body decay K L → 3π 0 . Although this Kaon decay mode has a large branching ratio, BR(K L → 3π 0 ) 0.2, the photon multiplicity is six. The impact of this background in the 4γ analysis will crucially depend on the photon inefficiencies of the Main Barrel detector ( γ ), which is at the level of 10 −3 [69,70]. We estimate the total efficiency of the 3π 0 background as where the factor of 30 is from combinatorics, and other are efficiencies other than photon one. Eq. (4.8) can be compared to the other physics backgrounds, K L → π 0 γγ and K L → 3π 0 γγ: BR(K L → π 0 γγ) × π 0 γγ 8 × 10 −6 and BR(K L → π 0 π 0 ) × π 0 π 0 8 × 10 −8 . Thus, K L → 3π 0 background is subdominant but can be larger than K L → π 0 π 0 , which may affect the sensitivity at low mass m a 100 MeV (see Fig. 3). The simulation of K L → 3π 0 background requires the full detector simulation, which is beyond the scope of the paper.

Displacement and energy of axion decay
Light ALPs tend to be long-lived because hadronic final states are kinematically forbidden. Decay with up to 5 cm displacement is effectively prompt decay in our analysis for the KOTO experiment. In fact, the resolution of the reconstructed vertex due to the finite energy resolution of the photons is typically 5 cm. This is shown in Fig. 5, where a comparison between the location of the true and reconstructed vertex location is shown for different axion masses.
In order to remove displaced decay with a displacement larger than 5 cm 7 , each signal where τ a is the ALP mean life-time and (E a /m a ) is the relevant boost factor, βγ. The ALP energy distribution used to compute the boost factor is shown in Fig. 6 for different values of the ALP mass. Figure 6. Distribution of the energy of the ALP produced in K L → π 0 a decays for different ALP masses after cuts. The original MC sample size is 5 × 10 4 , and the different normalization of the several curves indicate the difference in the cut efficiency depending on the ALP mass.

Expected sensitivity to the four photon search
The solid blue line in Fig. 7 shows the future reach of the KOTO experiment to the K L → π 0 a → 4γ signature, as a function of the ALP mass for both m a < m π and m a > m π . We assume that the sensitivity is determined by the statistical uncertainty, and the decay of the ALP to di-photon is treated as a prompt decay. To obtain this curve, we require S > 2 √ B where S is the number of signal events from K L → π 0 a and B is the number of background events from K L → π 0 π 0 , π 0 γγ. The analysis is done with the future KOTO luminosity of N 0 K L = 2 × 10 14 . From the figure we observe that KOTO can be sensitive to branching ratios as small as few×10 −9 .
This proposed search can have systematic uncertainties from the determination of the SM background. In Fig. 7, we show the cases of 1% and 10% systematic uncertainties as dashed purple and green lines, respectively. We expect these two curves to be very conservative. In fact, the expected signal has a reasonably narrow peak shape, which will allow data-driven background subtraction such as side-band technique.

Axion simplified models
In this section, we study the KOTO sensitivity to the K L → π 0 (γγ)a(γγ), four-photon final state, in terms of several ALP simplified models. We also compare the reach to other past and present high intensity experiments.

Introduction to the model
We consider a simplified model where the ALP couples only to the field strengths of the SU(2) W gauge bosons: where W a µν is the SU (2) field strength tensor,W a αβ = 1 2 αβµν W a,µν , and the g aW coupling is the leading term in the EFT expansion. This coupling is responsible of Kaon decays into ALPs, through W-loop penguin diagrams. In particular, the charged and neutral Kaons will have a decay width [71]: where The effective coupling g asd is given by . In our numerical analysis, we use the CKM elements as taken from the CKMfitter Group [72]. In Fig. 8 we show the branching ratio of K L → π 0 a as a function of the ALP mass, as well as of the g aW coupling (gray dashed curves). In this scenario, the branching ratio of K + → π + a is correlated with the K L one through the isospin relation, BR(K + → π + a)/BR(K L → π 0 a) ∼ 1.8 .
Once produced, the ALP will decay back to SM particles. In particular, below the pion mass, the axion will decay to photons with a width: In Fig. 8 we show the proper lifetime of the ALP in meters (red curves). Similarly, after electroweak symmetry breaking, the ALP will also couple to Zγ and ZZ. In particular, g aZγ 4 aZ µνF µν , g aZγ = g aW sin θ cos θ; g aZZ 4 aZ µνZ µν , g aZZ = g aW cos 2 θ . (5.6) As we will discuss in the next section, the former coupling can induce a signal at the LEP experiment, since it induces an exotic decay of the Z boson, Z → γa, with width given by

KOTO sensitivity and comparison with other experiments
The KOTO model independent bound presented in Fig. 7 can be interpreted in terms of this ALP simplified model. We compute an "effective branching ratio" for K L → π 0 a(γγ) from Eq. (5.3), taking into account the probability for the ALP to decay within 5 cm from the Kaon decay vertex: where τ a is the proper lifetime of the ALP, as shown by the red curves in Fig. 8. γ a is the boost factor of the ALP that can be easily extracted from Fig. 6. The corresponding reach is shown in the right panel of Fig. 9 by the region delimited by the red dashed line. This bound corresponds to the "Future Sensitivity" shown in Fig. 7 for the model independent bound. The bound is relatively flat above the pion mass and at around g aW ∼ (5 − 8) × 10 −5 /GeV. It becomes quite weaker at ALP masses m a 50 MeV, because of the weaker bound on the BR(K L → π 0 a → 4γ) (see Fig. 7) and because the life time of the ALP becomes quickly macroscopic.
We now compare this bound to the bounds that we can obtain from other present, past, and future high intensity experiments, and, in particular, with the NA62 experiment. For additional phenomenological analyses of similar benchmark scenarios, see e.g. [73].

Past Kaon experiments
Other past experiments looked for an ALP produced from either charged or neutral Kaon decays. The charged Kaon experiments E949 and NA48/2 set an upper bound on the Present and future bounds on the parameter space. In gray, we present the present bound (as shown in the left panel); in red and magenta, and in purple and blue, we present the future bounds at KOTO (4γ and 2γ + invisible signatures), and at NA62 (π + + 2γ and π + + invisible signatures).
branching ratio of K + → π + γγ [74,75] that can be used to set a constraint on a prompt ALP. Similarly, the E949 and NA62 bounds on the SM K + → π + νν decay [65,76] can be reinterpreted in terms of a constraint on a long lived ALP. Finally, the KTeV analysis for K L → π 0 γγ [64] can be utilized to set constraints on a prompt ALP, and the KOTO analysis for K L → π 0 νν [61] to set constraints on an invisible ALP.

• NA48/2, π + γγ analysis
We utilize the NA62/48 measurement of K ± → π + γγ in the kinematic range z = (m γγ /m K ) 2 > 0.2 [75] to set a bound on the ALP parameter space. Our analysis is similar to the one done in Ref. [71], even if we use a different statistical method. In particular, as a conservative bound, we require that the expected signal is less than the observed data plus two sigma uncertainty. We use Fig. 4 of [75] to set the bound on the branching ratio as a function of the ALP mass for m a ∈ (220 − 350) MeV. We require that the ALP decays in the detector volume, and, more specifically, that the decay length in the lab frame is less than 10 m. We include the corresponding weight factor (1 − exp[− 10 m τa(Ea/ma) ]) where E a is taken to be 37 GeV (i.e. half of the Kaon energy). Our bound is shown in violet in the left panel of Fig. 9.
• E949, π + γγ analysis The E949 experiment searched for K + decays at rest with a pion momentum p π + > 213 MeV. This analysis was re-interpreted in terms of K + → π + a, a → γγ with the ALP decaying within 80 cm of the stopped Kaon [71]. The corresponding bound is shown in purple at m a < 110 MeV in the left panel of Fig. 9.
• NA62, π + + invisible analysis In the K + → π + νν analysis [65], there are two distinct signal regions at low and at high missing mass: 0 MeV < m miss < 100 MeV (R1) and 161 MeV < m miss < 261 MeV (R2), respectively. We calculate the acceptance of K + → π + a(invisible), A π + a = 5.2% and 7.0% in R1 and R2 signal regions, respectively 8 . In addition, to compute the total yield, we adopt the same trigger efficiency, trigg = 0.90, and veto efficiency, veto = 0.76 as for the SM decay. The decay K + → π + a is effectively a K + → π + + invisible decay at NA62, as long as the ALP has a decay length of at least 150 m. Therefore, we compute the number of signal events in the two signal regions as N K + × BR(K + → π + a) × A π + a × veto × trigg × e −150 m/τ γ , where N K + is the number of Kaons decaying in the fiducial region (N K + ∼ 1.2 × 10 11 with the present dataset). The mean boost, γ , is calculated for each ALP mass using our signal Montecarlo events that pass the cuts on the geometrical acceptance for the charged pion, assuming that all Kaons are produced with an energy of exactly 75 GeV. Based on the NA62 observed number of events, we require the number of signal events to be less than 3.0 in the R1 signal region and less than 4.74 in the R2 signal region. The corresponding bound is shown in blue in the left panel of Fig. 9.
• E949, π + + invisible analysis The E949 collaboration has interpreted their analysis for the SM K + → π + νν decay in terms of a bound on a new stable massive particle produced from K + → π + X [76]. We utilize this result to set a bound on our ALP parameter space. We require that the effective branching ratio for K + → π + a is smaller than the one presented in Fig. 18 of [76]. To compute this effective branching ratio, we compute the probability for the ALP to escape the detector, i.e. to have a life-time longer than 1.5m, starting from a Kaon decaying at rest 9 . Our bound is presented in cyan in the left panel of Fig. 9.
• KTeV, π 0 γγ analysis The KTeV analysis for K L → π 0 γγ [64] has been utilized to set a bound on a prompt (i.e. decaying within 1 m from the K L decay) ALP decaying into two photons [71]. The corresponding bound is shown in red in the left panel of Fig. 9.
invisible NP particle with mass below ∼260 MeV. The analysis utilizes 3.68 × 10 11 K L decaying inside the detector (see Table 1). Branching ratios as small as ∼ 2.2 × 10 −9 have been tested under the assumption of a 100% invisible decay. We reinterpret this search in terms of our ALP model, requiring that the ALP has a lifetime long enough to decay after the detector. We obtain the distribution of the energy and decay point of K L based on our MC simulation with the analysis defined in [61]. Our bound is shown in pink in the left panel of Fig. 9.

Past Colliders and beam dumps
In addition to Kaon experiments, the LEP and Tevatron colliders also set constraints on the parameter space of this model. This is shown by the two green regions in the left panel of Fig. 9.
Because of the ALP g aZγ coupling (see Eq. (5.6)), the Z boson can decay to aγ, inducing a multi-photon signature in the LEP detectors [77,78]. The total and differential cross sections for the process e + e − → γγ was measured by the L3 collaboration at around √ s = 91 GeV [79]. In particular, the L3 experiment set the bound BR(Z → γγ) < 5.2 × 10 −5 . This bound is directly applicable to our model at light ALP masses, since the photons from the ALP decay would be collimated in the L3 detector. The dark green region in left panel of Fig. 9 is the bound we obtain from this branching ratio, asking the ALP to decay before the L3 ECAL (and therefore with a decay length smaller than ∼ 0.5m [80]).
Similarly, the CDF collaboration searched for the decay of a Z boson into two photons [81]. In particular, the collaboration set a bound BR(Z → γγ) < 1.46 × 10 −5 and BR(Z → γπ 0 ) < 2.01×10 −5 , with the pion detected as a single photon. We apply this more conservative bound on the decay into a photon and pion for ALP masses below the pion mass. To obtain the corresponding bound, we require the ALP to decay before the CDF central electromagnetic calorimeter located at 6.8 in from the collision point [82] 10 . The bound is represented by the dark green region at m a < m π , in left panel of Fig. 9.
Finally, past electron and proton beam dump experiments set a bound on the coupling of the ALP with photons g aγγ (see Eq. (5.5)). We take these bounds from [47]. They are represented in gray in the left panel of Fig. 9.

Future measurements at Kaon experiments
Next, we compare the future sensitivity of KOTO to the four-photon final state (see red region in the right panel of Fig. 9) to other projection of searches of NA62 and KOTO.
In particular, the purple region in the figure represents our projection of the NA48/62 K + → π + γγ analysis. To produce this region, we scale the NA48/62 K + → π + γγ uncertainty by the √ L where L is the ratio of NA62 and NA48/62 number of Kaon decaying in the fiducial volume. The NA62/48 have used 1.59 × 10 9 K ± decays in the fiducial volume [75], while for the future projection we assume that the NA62 will collect 10 13 K + decays with a downscaling trigger factor of 400 for K + → π + γγ [65].
The blue region in the figure represents the projection of the NA62 π + +invisible bound utilizing the full future luminosity. The bound corresponds to 12 events obtained with 10 13 K + decaying in the fiducial region 11 .
Finally, the pink region in the figure represents the projected bound for the KOTO K L → π 0 + invisible analysis. To obtain this bound, we scale the bound on the branching ratio in [61] with the √ L (L = 1.6 × 10 13 /3.68 × 10 11 , see Table 1).

Gluon coupled axions 5.2.1 Introduction to the model
The axion solution to the strong CP problem makes benchmark scenarios with an ALP coupled to gluons particularly interesting. The effective Lagrangian at the low energy scale µ ∼ m a 12 , can be written as where F a is the ALP decay constant andG a µν = 1 2 αβµν G a αβ . Since the effective theory can be valid up to a scale 4π/g ag = 4π(2πF a /α s ) ∼ 4TeV(F a /10 GeV) (while the cutoff can be 4πF a in many models), we focus on the ALP phenomenology and ignore the bounds from heavy states.
To obtain the form of the effective theory below the Λ QCD scale we resort to chiral perturbation theory. For convenience, we perform a chiral rotation of light quarks to remove the aGG coupling [83] (see also [35,45]) and generate the derivative couplings with the three light quarks (up, down, and strange) at leading order in the chiral Lagrangian: In our analysis, we also keep the strange quark, since, as we further discuss below, the mixing with the η meson also plays an important role when computing Kaon to ALP decay processes [84].
These derivative couplings induce a kinetic mixing between the SM mesons and the ALP. The π 0 and η states receive a small admixture of the physical ALP state, such that π 0 π 0 phys + θ πa a phys , η η phys + θ ηa a phys , (5.11) where, at the leading order, the mixing angles are given by , (5.13) 11 We have obtained 12 events via rescaling the number of SM single event sensitivity (0.267) and background (0.152) events observed now by NA62 with 1.2 × 10 11 K + [65]. 12 In Appendix C.5, we will briefly discuss additional UV contributions that can affect the K → πa rate if this effective Lagrangian is, instead, valid at a higher energy scale.  where we have defined and F π is the pion decay constant given by F π ≈ 93 MeV. θ ηη is the η-η mixing, whose value has a large uncertainty and lies in the range −(10 • -20 • ) (see e.g. [85][86][87]). Note the different m a dependence in the ALP-η mixing of the cos θ ηη and sin θ ηη terms. This is due to the fact that the sin θ ηη term arises from mass mixing, the cos θ ηη from kinetic mixing. At the same order in the chiral Lagrangian, the physical masses of the ALP, pion, and eta mesons are unaffected.
From the ALP mixing with neutral light mesons and the known operators for hadronic decays of the Kaons in the chiral Lagrangian (see Appendix C), we can calculate the Kaon decay widths at the leading order (similar calculations can be found in [88]). For simplicity, in the following we will fix sin θ ηη = −1/3 [49]. We will comment in the text, how the results will change if we had fixed a different value of θ ηη in the −(10 • -20 • ) range.
where the CP violating parameter in the Kaon mixing is given by K = 2.23 × 10 −3 , and | p a | is the absolute value of the momentum of the ALP. The corresponding effective couplings are The G 8 and G 27 couplings are the coefficients in front of the operators responsible for thes →d transition, which transform like (8 L , 1 R ) and (27 L , 1 R ) (see Appendix C). From lattice calculations, we know that the G 8 coefficient is significantly larger than G 27 . In our numerical analysis we will use the leading order values [67], where V ud , V us are CKM elements. In Fig. 10, we show the BR(K L → π 0 a) and BR(K + → π + a) as a function of m a and of the decay constant F a . As we expect from the K suppression in Eq. (5.15), the branching ratio of the neutral mode is generically suppressed, if compared to the one of the charged mode. There are also some accidental cancellations of the charged Kaon branching ratio. The position of the cancellation at low mass m a ∼ 80 MeV largely depend on the particular value chosen for θ ηη , whose uncertainty is sizable. The position of the cancellation at higher mass m a ∼ 210 MeV, instead, depend importantly on both the exact values of the quark masses, and the mixing angle θ ηη .
The decay of the ALP is controlled by the di-photon coupling and it is generated by the chiral rotation and the mixing with the mesons, g aγγ 4 aF µνF µν , The trace runs on the three-flavor space, Q q is the diagonal matrix with the electric charges of the quarks on the diagonal, N c = 3 is the number of colors, λ 3 , λ 8 are Gell-mann matrices (normalization tr[λ a λ b ] = δ ab ), and λ 0 = 1 √ 3 diag{1, 1, 1}.

KOTO sensitivity and comparison with other experiments
In Fig. 11, we show the current bounds (left panel) and future reach (right panel) on the parameter space of this simplified model. We compare the bound from the KOTO experiment to the bounds from other Kaon experiments, as well as other present and future accelerator experiments. The discussion for the bounds and projections is almost parallel to Sec. 5.1.2 for the SU(2)-coupled ALP simplified model. The most relevant differences arise for LEP, beamdump experiments, the GlueX experiment, and PIBETA experiment, which we comment in the following.

• LEP
The GG coupled ALP does not have a coupling to Zγ unlike the SU(2) coupled ALP. Still LEP set a constraint on this benchmark model through the process e + e − → γ * → γa where the di-photon from the ALP decay is collimated and seen as a single photon.
In [89], the bound on the aFF operator was derived from the OPAL inclusive 2γ search [90]. We show this bound in dark green in the left panel of Fig. 11.

• Proton and electron beam dump experiments
In the proton beam dump experiments, the GG coupled ALP can be produced through the meson mixings and decay by the effective photon coupling. The bound was studied in Ref. [91] using the CHARM result [92]. In our figure, we also include the bound from the electron beam dump experiments, E141 and E137, where the induced photon coupling is responsible to both the production and decay (see [47] and references therein). Both bounds are shown in gray in the left panel of Fig. 11.

• GlueX experiment
The GlueX experiment can be used to set a bound on the ALP parameter space [50]. The experiment utilizes a 9 GeV photon beam colliding against a fixed target. The ALP can be produced from the decay of vector mesons such as ρ and ω and observed via its decays to photons. The bound was derived using 1/pb data of [93], and it is shown in yellow in the left panel of Fig. 11.
• PIBETA experiment The precision measurement of π + → π 0 (→ γγ)eν at the PIBETA experiment [94] gives a constraint on θ πa for 100 MeV m a m π 0 [95]. The corresponding bound is shown in light green in the left panel of Fig. 11 ("πβ"). We have checked that the measurement of π + → eν by the PIENU collaboration [96] does not give an additional constraint in the region of parameter space shown in Fig. 11 [95].
The main updates in Fig. 11 are the bounds from Kaon decays. In the figure, we only present the bounds obtained from visible searches (K ± → π + γγ and K L → π 0 γγ), since the invisible ones do not extend the reach of the beam dump experiments. We represent each Kaon bound in the figure with a band. This quantifies the uncertainty coming from varying the quark mass ratios in the range m u /m d = 0.47 (+0.06, −0.07) MeV, m s /m = 27.3 (+0.7, −1.3) MeV [97]. This uncertainty particularly affects the NA62 bound at m a ∼ 230 MeV, close to the accidental cancellation for BR(K + → π + a) shown in Fig. 10. We do not show the uncertainty on the bound coming from the uncertainty in the determination of θ ηη . This will particularly affect the E949 [K + ] bound at low mass. This bound on F a can change by a factor of ∼ 2 varying θ ηη in −[10 • , 19 • ]. Note that the bounds on the K L decays do not suffer of large uncertainties, as we discuss in more details in Appendix C.4.
As shown in the right panel of Fig. 11, in the future, both the NA62 K + → π + a ("[2γ]"), and the proposed KOTO K L → π 0 a ("[4γ]") searches will significantly extend the probed parameter space, especially at m a > m π 0 . It is also interesting to note that, in the future, the parameter space of the gluon coupled ALP will be fully probed up to decay constants F a ∼ few TeV for m a m π . The regions of parameter space not yet probed in the left panel of Fig. 11 will be, in fact, probed by the GlueX experiment with 1/fb data [50] and by the SeaQuest experiment at Fermilab [98].

Breaking the Grossman-Nir bound
In this section, we describe in details the GN-breaking simplified model introduced in Sec. 2.2, highlighting the unique sensitivity of the KOTO experiment in probing its parameter space.

Chiral Lagrangian analysis
The effective scalar (φ)-quark couplings given in Eq. (2.10) can be embedded in the chiral Lagrangian in the form of mass terms as where B 0 = m 2 π /(m u + m d ), F π is the pion decay constant, F π ≈ 93 MeV, y 1,2 are real couplings in the mass basis, v is the vacuum expectation value of the Higgs, v = 246 GeV, and Σ is the common exponential pion field matrix λ sd is the three by three matrix leading to s → d flavor violation Expanding the Lagrangian in (6.1) in powers of Π, we find the matrix elements for the several meson to φ transitions: These terms will lead to exotic K + , K L and K S decays, where we are working on the phase convention K L On top of this effective Lagrangian, we add the effective operator in Eq. (2.12), χ Λχ F µνF µν , that is responsible of the χ decay into two photons. As long as Λ χ 50 TeV mχ 120 MeV 2 , the decay length of χ is smaller than ∼ 10 cm in the analyzed mass range, and, hence, its decay is effectively prompt [6].

New Kaon decays
From the effective Lagrangian in (6.4), we can compute the matrix elements for the several transitions. We find (6.5) where | p σ | is the absolute value of the momentum of σ in the center of mass frame. Also the charged Kaons will inherit new exotic decay modes. However, due to charge conservation, only decay modes with three (or more) final states will be generated (see (6.4)). In particular, Analogously, K ± can also decay to π ± χχ with the amplitude given by (6.6) with the replacement dm 2 σσ dm 2 σπ → dm 2 χχ dm 2 χπ . K L also acquires new three-body decays: and correspondingly for K L → π 0 χχ. Finally, new (two or three-body) decays of K S will be also induced: and, similarly, one can obtain the width for K S → χχ with the replacement | p σ | → | p χ |. These decay modes will be obviously more suppressed due to the larger width of K S .
The loop is quadratically sensitive to the internal momentum, Λ cutoff . The loop momenta that characterize the pion-photons coupling decrease significantly above the QCD scale. Therefore, the above estimate of the χ − π 0 mixing shows that this effect can be neglected. As for σ, it can decay to four photons (e.g. via its coupling to χ and a neutral Kaon which couples to two photons) however this coupling is suppressed by CKM factors, extra loop and 1/Λ χ . Therefore it is safe to consider σ effectively stable.
The efficiency for K L → σχ to end up in the KOTO signal region depends crucially on the mass of the σ and χ particles. In the left panel of Fig. 12, we show in blue the efficiency has a function of m χ , that, for convenience, we fix to be = m σ . A sizable efficiency is reached as long as the χ mass is not too far away from the mass of the pion. In Fig. 13, we also show the distribution of our montecarlo events for K L → σχ for different values of the m χ = m σ mass. As we can observe, the events fall nicely in the signal region (the region delimited in red) as long as 100 MeV m χ = m σ 160 MeV.
Using the efficiency of the left panel of Fig. 12 and the widths discussed in the previous section, we can compute the sensitivity of KOTO to our model, as well as the corresponding predictions for the other exotic K + and K S decay modes. In the right panel of Fig. 12, the blue lines represent the BR(K L → σχ) needed to produce 3 events in the KOTO signal region using the data collected in 2016-2018 (solid line), or future KOTO data (dashed blue). Note that 2016-2018 data is already able to probe a branching ratio as small as BR(K L → σχ) ∼ 1.3 × 10 −9 . This corresponds to a GN breaking scale as high as Λ GNV / (y 1 + y 2 ) ∼ 10 7 GeV.  The other lines in the right panel of Fig. 12 are the corresponding predictions for BR(K L → π 0 σσ) and BR(K L → π 0 χχ) (light blue), BR(K ± → π ± σσ) and BR(K ± → π ± χχ)(red), BR(K ± → π ± σχ) (yellow), BR(K S → σσ) and BR(K S → χχ) (green), and BR(K S → π 0 σχ) (purple), once we demand the model to produce 3 events in the KOTO signal region using the data collected in 2016-2018. For the latter three curves, we have fixed y 2 = 2y 1 (see the parametric dependence of the several widths discussed in Sec. 6.2).

Discussion and overview
Rare Kaon decay modes have been always considered among the few holly grails of flavor physics because of their rareness, and because of our ability to control them well theoretically. This made rare kaon decays singular in their ability to probe new physics (NP) models.
The current time is rather unique as both the KOTO and the NA62 experiments are collecting high quality data aiming to reach unprecedented precision in the measurement of neutral and charged Kaons, respectively, providing a direct test to the SM predictions. What makes all this possible is the huge fluxes of Kaons achieved at J-PARC and at CERN. The fact that the KOTO and NA62 detectors have access to these astronomical fluxes makes them sensitive to other types of dynamics that we denote as "dark sector physics". By dark sectors we refer to a class of models with light particles that couple only weakly to the Standard Model (SM) fields. In this work, we have shown that such dark sectors can be probed in regions that could not have been searched for so far, which is rather exciting.
In this paper, we particularly highlight the complementarity of the two experiments. Naively, the fact that NA62 already probes charge Kaon decays with branching ratios as small as 10 −10 , while KOTO is an order of magnitude behind in the corresponding neutral decay mode, makes one conclude that KOTO is only providing us with a secondary validation of the searches done at the NA62. This statement may be enforced by the Grossman-Nir (GN) relation that bounds the size of the NP contributions in the neutral mode via the charged one. We, however, demonstrate that the physics of dark sectors do not necessarily follow this pattern. The reason is two fold: i) on the experimental side, the two experiments are different in several essential aspects, in terms of kinematics, acceptance and sensitivity to different final states; ii) on the theoretical side, when examining dark sectors one find that the NP-GN relation can be effectively violated by as much as several orders of magnitude.
We show this by considering two qualitatively different physics cases. In the first, we consider models of axion-like-particles (ALPs) which couple to electroweak gauge bosons or to gluons. We find that ALP-diphoton decay mode, K L → π 0 (γγ)a(γγ), can be very efficiently searched for at KOTO. We layout a new search strategy that, if adopted by the collaboration, would allow KOTO to probe uncharted territories of ALP-physics. At the same time we also find that the corresponding final state at NA62, K + → π + a, while being equally interesting, can suppressed. This is the case of the ALP-coupled to gluon model where cancellations between different contributions when including the η and η contributions can happen, albeit with large theoretical uncertainties. This probably calls for a more detailed theoretical analysis of the K + → π + a decay, going beyond leading order in the chiral Lagrangian, and also carefully including the uncertainties related to quark masses and to the η − η mixing. This might be an interesting study to be performed on the lattice which would then be freed from the uncertainties related to the chiral expansion. In the second, we introduce a model based on approximate strange flavor symmetry, that effectively leads to a strong violation of the Grossman-Nir bound. We find that this benchmark model can be discovered by the KOTO experiment looking for K L decaying into two photons plus invisible. It is also worthwhile to mention that this benchmark could also account for the potential-candidate events seen at KOTO and the absence of signals at NA62 at the same time. In fact, the corresponding charged Kaon signals at NA62 would be several orders of magnitude suppressed, and effectively hidden to this latter experiment.
Our main messages of this paper are quite general and motivate model independent searches at both Kaon factories.
Note added: while this work was at its final stage of completion, Refs. [99][100][101][102] that have some overlap with with the topics discussed above, appeared.
3. Pick the combination where the two reconstructed vertices, Z vtx,1 and Z vtx,2 , are the most consistent. A "pairing variance" is introduced to evaluate the consistency of the two vertices, where n pair = 2 for four photon events, and dz is the distance from the ECAL (Z vtx = 6.148m−dz). The combination that minimizes χ 2 dz 2 is picked, and then the decay vertex of K L is identified. Since dz 2 = dz 2 ( r 1 , r 2 , E γ 1 , E γ 2 ), the variance is obtained by the combination of resolutions of photon position and energy: where σ E and σ r are the energy and position resolution, respectively.

B Validation of our analysis
We validate our simulation and reconstruction algorithm by cross-checking the measured quantities at KOTO [61,62].

B.1 Detector effects
First, we check the the detector smearing we include in our analysis reproduces the KOTO results. At KOTO, π 0 → γγ was measured in a special run and the shape of the di-photon invariant mass was reported in Fig. 7 of [62]. Assuming the decay vertex is known, we reproduce the shape with the energy and position resolution. We use the detector parameters given in [62]. We find that the position resolution is only a minor effect. This implies that, when the vertex is reconstructed, the main source of uncertainties is from ECAL smearing.

B.2 Reconstruction of four photons
K L → 4γ was measured by the KOTO collaboration, and the reconstructed four photon invariant mass is shown in Fig. 11 of [62], . The peak region is dominated by K L → π 0 π 0 . We simulate K L → π 0 π 0 events and perform the SM reconstruction discussed in Appendix A. The shape of the peak region is well reproduced.
Ref. [62] provides the acceptance for the performed analysis, A = 1.48 × 10 −3 . We reproduce this acceptance at the 10% level, as it is shown in Table 2 Table 2. Decay probability, cut flows, and acceptance of our K L → π 0 π 0 analysis.

B.3 Analysis of the decay K L → π 0 νν
We have used the KOTO K L → π 0 νν result to set constraints on the two ALP benchmarks of Sec. 5, as well as for the K L → σχ analysis. First, we cross checked the efficiency of K L → π 0 X(→ invisible) with [6]. We calculated the KOTO correction factor (K L → π 0 a). Also, we checked the overall acceptance. With the efficiencies of veto and shower cut ( veto = 17% and shower = 52%), we get an acceptance ∼ 30% higher than the reported acceptance reported in [61]. This information is useful to estimate the uncertainty of our proposed K L → π 0 a → 4γ search, where the result is fully based on our simulations.

C.1 ALP-meson mixing
We want to compute the ALP-meson mixing arising from the the effective Lagrangian in (5.10). We work in the framework of the chiral Lagrangian and perform a chiral rotation such to remove the mass-mixing between the ALP and the light mesons [83]. The remaining ALP interactions with SM mesons is through the kinetic mixing that is given by whereκ q is the diagonal matrix with κ q = 1 mq / q ( 1 m q ) on the diagonal, and the Σ is the non-linear meson field. However, due to the non-negligible mixing between the η and η mesons, the η meson has a mass mixing with the ALP through the axial anomaly.
To compute this effect, we follow the prescription given in [49]. We keep only the light state η in the mass basis and decouple the η , assuming sin θ ηη = −1/3 (see, however, Eq. (5.13) for the generic expression of the ALP-η mixing). Then the non-linear meson field is given by where we are adopting the pion decay constant F π ∼ 93 MeV. The chiral Lagrangian has now both kinetic mixing terms of the ALP with the pion and the η, and mass mixing terms of the ALP with the η: where m is the matrix m = exp(iκ q a 2Fa γ 5 )m q exp(iκ q a 2Fa γ 5 ) and B 0 = m 2 π /(m u + m d ). After diagonalizing this system, the physical ALP and meson eigenstates are given by where

C.2 ∆S = 1 transitions
Based on Cirigliano et al [67] (see also references therein), at the low energy the two operators responsible for ∆S = 1 transitions are 13 The coefficients G 8,27 can be determined by the measurement of Kaon decays to pions (see Eq. (5.18) for their value). In order to study the width of K → πa arising from ALP-meson mixing, we need to obtain the trilinear interactions of K-π-π 0 /η. Below, we expand the two relevant terms in the chiral Lagrangian. 13 As mentioned in Sec. 6.1, we are working on the phase convention KL K 0 +K 0 √ 2 and KS the G 8 contribution coming from the η-ALP mixing also leads to an important contribution [84,88].
In the left panel of Fig. 14, we numerically compare the naive estimate from Eq. (C.15) (pink line) to our full result (black line). The result obtained keeping only the mixing θ πa (θ ηa ) is also shown in blue (red). A large difference between the pink and the black lines is particularly observable at low values of m a . This is due to the fact that the η contribution comes in part from ALP-η mass mixing that, contrary to kinetic mixing, does not go to zero for m a → 0. Furthermore, generically, the two new G 8 contributions from θ πa and θ ηa have a similar size and can lead to large cancellations depending on the value of the ALP mass. However, the value of m a at which this cancellation happens strongly depend on the various parameters. In Fig. 11 of the main text, we estimated the uncertainty on the K + experimental bounds coming from the uncertainty on the quark masses. In addition, there are additional sizable uncertainties coming from the uncertainty on the η-η mixing angle.
On the other hand, the prediction on BR(K L → π 0 a) is rather stable against these uncertainties because the dominant contribution comes from the G 8 contribution from the ALP-pion mixing, and even the naive formula analogous to Eq. (C.15) can typically capture it (see the right panel of Fig. 14). Therefore, the prediction on BR(K + → π + a) has a large uncertainty while the prediction of BR(K L → π 0 a) is more theoretically stable.

C.5 Possible UV contributions to K → πa
Our analysis in Sec. 5.2 for the ALP coupled to gluons assumed that the effective Lagrangian in (5.10) is given at the low energy scale µ ∼ m a . Starting with this Lagrangian, we have shown that the BR(K L → π 0 a) will be quite suppressed, if compared to BR(K + → π + a) because the former is CP violating (see the K suppression in Eq. (5.15)).
However, UV completions of this effective Lagrangian generically lead to an enhancement of the branching ratio of the K L mode at the two-loop order with direct CP violation, which we schematically described in the following.
The aGG coupling will also induce a coupling of the ALP with quarks. This is given by g eff aqq aqγ 5 q, g eff aqq = − α s π m q g ag log µ 2 m 2 q − 11 3 + g(τ q ) , (C. 16) where τ q = 4m 2 q /m 2 a and g(τ ) is a loop function that can be found e.g. in [45]. Thanks to these induced coupling, the ALP will be produced in K L → π 0 a through penguin diagrams with the ALP radiated from the quark loop. The corresponding partial width can be estimated as [103] Γ(K L → π 0 a)≈ G 2 where f L is the form factor in the K L decay, and V αd , V αs are CKM factors. This decay is induced, for example, by the dimension five operator ∂ µ a(s L γ 5 γ µ d L ) as in the SU(2) coupled ALP case of Sec. 5.1. In this case, we will not have a K suppression in the decay. Thus, this two-loop contribution can potentially significantly enhance the BR(K L → π 0 a) (see also [104]). Similarly, also the rate of the charged mode K + → π + a can be enhanced by these UV contributions: The enhancement is, however, not as sizable as the one in the neutral mode, because of the absence of the K suppression in the IR contribution.