§ REPRESENTATIVE MODELS

We start by reviewing the different BSM models featuring light new states that can be produced in rare hyperon and kaon decays. For each of the model we recast the current experimental constraints and comment briefly on the future reach, with a more detailed discussion of experimental issue relegated to Section <ref>. Whenever possible, we set well-motivated target branching ratios, explore their phenomenological implications, and comment on the implication of kaon and hyperon measurements on new physics scenarios addressing the shortcomings of the SM.

§.§ Higgs portal scalar

Authors: Gori, Egana-Ugrinovic, Homiller, Knapen, Meade
[For each section we list the contributing authors, with the corresponding author underlined.]

The Higgs portal is the minimal extension of the SM by one additional singlet scalar particle, $S$ <cit.>.
In addition to its simplicity, the Higgs portal has received considerable attention both because it can serve as a simple, renormalizable dark sector portal <cit.>
and because it is similar to UV-motivated models with light dilatons or radions <cit.>. The Lagrangian for the minimal Higgs portal model is given by

\begin{equation}
\mathcal{L} = \mathcal{L}_{\textrm{SM}} + \frac{1}{2}\partial_{\mu} S \partial^{\mu} S - V(S, H) .
\end{equation}

The full potential $V(S,H)$ can be found e.g. in Ref. <cit.>. For rare kaon and hyperon decays the only relevant feature is that the potential gives rise to mixing in the mass matrix between the singlet $S$ and the CP-even, neutral scalar component of the Higgs doublet, $H$. The resulting mass eigenstates are identified as the 125 GeV Higgs boson $h$, and a new scalar boson $\varphi$ with mass $m_{\varphi}$. The couplings of $\varphi$ to the SM fields are equal to the SM Higgs couplings up to a multiplication by a universal mixing angle $\theta$ (or rather by $\sin\theta \simeq \theta$). Therefore, both the production and decay rates of the scalar are each proportional to $\theta^2$. Despite its simplicity, the Higgs portal scalar has interesting phenomenology — in particular at beam dumps, meson factories and other low energy experiments — since its decay is governed by very small inherited couplings proportional to the SM electron and muon Yukawa couplings, while its production is governed by inherited couplings to heavier quarks.[The potential $V(H,S)$ also leads to couplings between the new singlet and the SM-like Higgs, leading for instance to exotic Higgs decays to a pair of scalars. Searches for these decays may be relevant for heavier singlets, but on general ground it can be shown that a large amount of fine-tuning would be required to have sizeable ${\cal B}( h \to \varphi\varphi)$ for a singlet in the mass range of interest for rare kaon decays <cit.>.]

At kaon factories, the new scalar would be produced in the two-body $K^+\to\pi^+\varphi$ and $K_L\to \pi^0 \varphi$ decays via its coupling to up-type quarks and gauge bosons in SM-like penguin diagrams <cit.>.
The corresponding branching ratios for $m_X<300$ MeV are approximately equal to

\begin{equation}
{\cal B}(K^+ \to\pi^+\varphi) = 1.7\times 10^{-3}\sin^2\theta, \quad
{\cal B}(K_L \to\pi^0\varphi) = 5.7\times 10^{-3}\sin^2\theta,
\label{eq:BrKpiX}
\end{equation}

and decrease for larger $m_X$ values (for the explicit formulae, see Eqs. 4.4 and 4.5 of Ref. <cit.>). Contours of constant ${\cal B}(K^+\to\pi^+\varphi)$ are shown in Fig. <ref>.

In the $m_\varphi<m_K$ regime relevant for kaon experiments,  decays of the new scalar particle to charged lepton and pion pairs ($e^+e^-$, $\mu^+\mu^-$, $\pi^+\pi^-$, $\pi^0\pi^0$) should be considered, above the corresponding thresholds. Their partial decay widths are given by <cit.>

\begin{equation}
\label{eq:higgs_portal_widths}
\Gamma(\varphi \to \ell^+ \ell^-) = \sin^2\theta \frac{m_{\ell}^2 m_{\varphi}}{8\pi v^2} \bigg( 1 - \frac{4 m_{\ell}^2}{m_{\varphi}^2}\bigg)^{3/2},
\quad
\Gamma(\varphi \to \pi \pi) = \sin^2\theta \frac{3 | G_{\pi}(m_{\varphi}^2)|^2 }{32 \pi v^2 m_{\varphi}^2} \bigg( 1 - \frac{4 m_{\pi}^2}{m_{\varphi}^2}\bigg)^{1/2},
\end{equation}

where $v = 246$ GeV is the electroweak scale, the $\pi\pi$ decay width is computed using chiral perturbation theory, and the relevant form factor is $G_{\pi}(x) = \frac{2}{9}x + \frac{11}{9}m_{\pi}^2$. The $\varphi\to\gamma\gamma$ decay is always subdominant and can be neglected in practice.

For $2m_e < m_\varphi < 2m_\mu$, the scalar decays primarily to electron-positron pairs and is long-lived with a lifetime $c\tau = \mathcal{O}(10^{-2}~\textrm{m} / \sin^2\theta)$, therefore being largely invisible and experimentally indistinguishable from the $K\to\pi\nu\bar\nu$ decays. Above the di-muon threshold, $\varphi\to\mu^+\mu^-$ becomes the dominant decay mode and lifetime drops to $c\tau= \mathcal{O}(10^{-6}~\textrm{m} / \sin^2\theta)$. Above the di-pion threshold, $\varphi\to\pi\pi$ take over as the dominant decay modes, and the lifetime drops by another order of magnitude. As a result, the rest-frame decay length reaches $10.4~\textrm{m}$ for $m_{\varphi} = 300~\textrm{MeV}$ and $\sin\theta = 10^{-4}$. Further discussion of the lifetime and decay branching fractions can be found in Ref. <cit.>.

Constraints on the Higgs portal scalar in the mixing angle vs mass plane. The NA62 and E949 limits extend down to $m_\varphi=0$ with almost no deterioration.

Experimental constraints on the Higgs portal are aggregated in the $\sin\theta$ vs $m_{\varphi}$ plane in Fig. <ref>. Constraints from $K\to\pi X_{\rm inv}$ searches  obtained by the NA62, KOTO and E949 experiments are discussed in Section <ref> where $X_{\rm inv}$ represents a particle invisible in each detector, while the prospects for searches for visible scalar decays ($\varphi\to\mu^+\mu^-$) are discussed in Section <ref>. Constraints from the CHARM beam dump experiment at the SPS <cit.> are recast following Ref. <cit.>, including the effects of absorption in the target <cit.>,  geometric efficiency and $\eta$ decay contributions.
A similar but slightly weaker constraint from the NuCal/U70 experiment <cit.> is omitted for clarity.
In a similar fashion to beam dumps, complementary limits can be placed using neutrino experiments by searching for the visible decays of $\varphi$. Taking advantage of the large number of kaon decays at rest at the NuMI absorber, the MicroBooNE collaboration set strong limits on $\theta$ with a dedicated search for $\varphi\to e^+e^-$ <cit.>. Similar limits were obtained with a rescast of LSND <cit.> as well as of PS191 data <cit.>. Although not shown explicitly in the figure, T2K can marginally improve on the latter by using the existing search for HNLs <cit.>.

Constraints arising from searches for displaced vertices in $B \to \varphi(\mu^+\mu^-)$ decays at LHCb <cit.>, adapted from Ref. <cit.>, are the strongest in the $m_{\varphi} > 2m_{\mu}$ range. Searches in $B$ decays from Belle <cit.> and BaBar <cit.>, as well as constraints from the $K_L \to\pi^0\mu^+\mu^-$ decay at KTeV <cit.> are sub-leading in the parameter space considered.
Finally, the proposed SHiP experiment <cit.> is expected to be able to cover much of the parameter space for $m_\varphi > 2m_\mu$ down to mixing angles $\sin\theta\sim 10^{-5}$, and reach a sensitivity of $\sin\theta\sim 10^{-4}$ for smaller masses <cit.>.

There are also cosmological and astrophysical constraints on light scalars. If the singlet lifetime is too large, the decay of scalars after nucleosynthesis spoils the light element abundances. The resulting Big Bang nucleosynthesis (BBN) constraint adapted from Ref. <cit.> is shown in Fig. <ref>. The BBN constraint is sensitive to the metastable abundance of singlets, set via thermal interactions with the SM bath at temperatures above an MeV. This introduces a very mild dependence with the singlet-higgs quartic $S^2 H^\dagger H$, which is in one-to-one correspondence with the $h\to \varphi\varphi$ branching ratio. The shaded region indicates the weakest constraint from Ref. <cit.> at a given mass, while the dot-dashed line indicates the constraint for ${\cal B}(h\to\varphi\varphi)=10^{-3}$, which allows the singlet mass to not be fine tuned.
For even smaller values of the mixing angle, there are also constraints from the observation of SN1987a neutrinos, which sets a bound on the amount of energy that could have been carried away by scalars produced in the core. An older calculation of this limit in Refs. <cit.> was recently revisited in Ref. <cit.>, finding substantial disagreement.
The SN1987a constraints are sub-leading to the BBN bounds in the entire parameter space shown, and are not included in Fig. <ref>.

Thus far we have focused on the minimal Higgs portal scenario, when the scalar inherits all its couplings via mixing with the Higgs.
In many well-motivated UV scenarios, however, dimension-5 couplings between the scalar and gauge bosons (in particular, the photon) can arise.
This is the situation in models of light radions or dilatons, or in theories where there are additional vector-like fermions charged under the SM gauge group.
These dimension-5 couplings can compete with the small, inherited couplings of the scalar to light fermions, and significantly change the phenomenology, especially the decay patterns, of the scalar.
While the minimal Higgs portal for $m_{\varphi} < 2 m_{\mu}$ is long lived, with only a percent-level branching ratio to $\gamma\gamma$ compared to $e^+e^-$, more general scalars could decay predominantly to $\gamma\gamma$ and have much shorter lifetimes that significantly change the constraints from beam dump and fixed-target experiments.
While the underlying models in the scalar case are very different from the general pseudoscalar (or axion-like particle, ALP) case, the phenomenological signatures at kaon and hyperon experiments are virtually identical.
We thus leave a detailed discussion of light particles with varying couplings to photons to the next section.

In case $\varphi$ is a pure scalar, the $K\to\pi\pi\varphi$ is generically suppressed compared to $K\to\pi\varphi$ and the search in this channel can be ignored. If $\varphi$ is not purely CP even, as is the case in some well motivated theories <cit.>, then $K\to\pi\pi\varphi$ probes a different chiral structure of the light scalar couplings to the SM. This will be discussed in detail in Sec. <ref>.

§.§ QCD axion and axion-like particles (ALPs)

Authors: Alonso-Alvarez, Alves, Bauer,  Datta,  Gori, Kvedaraitė, Martin Camalich, Neubert, Redigolo, Renner, Schnubel, Soreq, Thamm, Tobioka, Ziegler

Axion-like particles (ALPs) are ubiquitous in BSM models. Whenever an approximate global symmetry is spontaneously broken, there is an associated pseudo Nambu-Goldstone boson (pNGB) which is often called an ALP in the recent literature. Being a pNGB, the ALP mass can naturally be much lighter than the UV scale of symmetry breaking. As such they are the phenomenologically interesting remnants of the high scale new physics dynamics, potentially accessible at experiments. Arguably the most motivated example of an ALP is the QCD axion, which can simultaneously solve the strong CP problem <cit.> and account for the observed dark matter abundance <cit.>. Already early on, in the early 1980s, it was realized that the QCD axion paradigm strongly motivates searches for massless ALPs at kaon factories <cit.>, and even more so if it is connected to a model of flavor <cit.> (for earlier discussions of flavor-violating axions see Refs. <cit.>).  When experimental searches are performed, however, one should keep an open mind and not needlessly limit the potential of such searches. For instance, non-minimal implementations allow for the strong CP problem to be solved by heavier axions than one might naively expect <cit.>. Moreover, other theoretically motivated scenarios such as low scale supersymmetry, composite Higgs models, models of dark matter freeze-out, and models of electroweak baryogenesis may contain light ALPs (see for instance <cit.>). The plethora of diverse theory motivations calls for a wide experimental exploration of the accessible ALPs masses and their couplings to the SM particles.

Let us assume that a global $U(1)$ symmetry is spontaneously broken by a condensate of size $f_a$ (such as the vacuum expectation value of a scalar field) and results in an ALP at low energies. The effective Lagrangian describing ALP interactions with the SM at lower energies
is given by

\begin{equation}
\mathcal{L}_{\text{ALP}}=\frac{(\partial a)^2}{2}-\frac{m_a^2a^2}{2}+\mathcal{L}_{\text{ALP-gauge}}+\mathcal{L}_{\text{ALP-F}} \label{L_ALP}\, .
\end{equation}

The interactions with the SM fields, $\mathcal{L}_{\text{ALP-gauge}}$, $\mathcal{L}_{\text{ALP-F}}$, are of dimension 5, and are obtained after the heavy fields with masses of the order $\Lambda_{\text{UV}}\sim 4\pi f_a$ are integrated out.
The ALP interactions include  in general couplings to all the SM gauge bosons: the gluons, $G^a$, the $SU(2)_L$ gauge bosons, $W^i$, and the hypercharge boson, $B_\mu$,

\begin{equation}\label{eq:alpgauge}
\mathcal{L}_{\text{ALP-gauge}}=\frac{N_3\alpha_s}{8\pi f_a} a G_{\mu\nu}^a\tilde{G}^{a\mu\nu}+\frac{N_2\alpha_2}{8\pi f_a} a W_{\mu\nu}^i\tilde{W}^{i\mu\nu}+\frac{N_1\alpha_1}{8\pi f_a} a B_{\mu\nu}\tilde{B}^{\mu\nu}\ .
\end{equation}

Here $\tilde G^{a\mu\nu}=1/2\epsilon^{\mu\nu\rho\sigma} G_{\rho\sigma}^a$, with $\epsilon^{0123} = 1$, and similarly for $\tilde W^{i\mu\nu}$, $\tilde B^{\mu\nu}$. The $N_i$ coefficients are  scale independent since they are related to the mixed 't Hooft anomalies involving the spontaneously broken global $U(1)$ symmetry and the SM gauge group.[In the QCD axion literature the $N_3$ coefficient is often absorbed in the definition of the decay constant $f_a$, $f_a/N_3\to f_a$, so that $\mathcal{L}_{\text{QCD axion}}=\frac{\alpha_s}{8\pi f_a} a G_{\mu\nu}^a\tilde{G}^{a\mu\nu}$.  We keep the $N_3$ factor explicit.]
In general, these receive contributions from both the SM fermions and from the heavy fermions that were integrated out at $\Lambda_{\text{UV}}$. After electroweak symmetry breaking the couplings to $W^i$ and $B_\mu$ gauge bosons result in ALP couplings to photons
ℒ_ALP-photon= α/8πf_a C_γγ a F_μνF̃^μν,
with $C_{\gamma\gamma}=N_1+N_2$, see also Eq. (<ref>) below.

The ALP couples derivatively to the SM fermions,

\begin{equation}\label{eq:ALP-f}
\mathcal{L}_{\text{ALP-f}} =
\frac{\partial_\mu a}{2 f_a} \, \bar f_i \gamma^\mu \big( C^V_{f_i f_j} + C^A_{f_i f_j} \gamma_5 \big) f_j  \, ,
\end{equation}

where $f=u,d,e,$ and the sum over repeated generational indices, $i,j=1,2,3$, is implied. $C^{V,A}_{f_i f_j}$ are hermitian matrices in flavor space which we take to be real for simplicity. Because of gauge invariance $C^V_{d_i d_j} - C^A_{d_i d_j} = V^*_{u_k d_i } (C^V_{u_k u_l} - C^A_{u_k u_l}) V_{u_l d_j}$, so that in total the ALP interactions with the charged SM fermions and gauge bosons are specified by 33  parameters. Below we also use the notation $F_{f_i f_j}^{V/A}\equiv 2 f_a/C_{f_i f_j}^{V/A}$.

For kaon decays the $s\to d a$ transition is the dominant ALP production mode even if the flavor violating couplings, $C_{sd}^{V/A}$, are highly suppressed compared to the flavor diagonal ones. The dominant kaon decay channels resulting in ALP production are thus

\begin{equation}
K\to \pi a, \quad
K\to \pi \pi a.
\end{equation}

The two-body (three-body) decays probe the vectorial (axial) ALP coupling $C_{sd}^V$ ($C_{sd}^A$), as can be easily shown using  the parity conservation of QCD interactions.

If we take the couplings in Eq. (<ref>) to be defined at the EW scale, $\mu_W = {\mathcal O}(M_Z)$, the renormalization group (RG) evolution from  $\Lambda_{\text{UV}}$ to $\mu_W$ should be taken into account and can be extracted from Refs.  <cit.>. The one loop corrections generate the CKM-suppressed flavor violating ALP couplings at the EW scale, even if the ALP couplings in the UV are flavor diagonal. The log-enhanced contribution, which dominates the one loop correction for $\Lambda_{\text{UV}}\gg M_Z$, follows from the anomalous dimension controlling the RG running and is given by

\begin{equation}
C_{ds}^{A,V}(M_Z)= \frac{y_t^2}{16\pi^2}V_{td}^*V_{ts}C_{tt}^{A} \log\frac{\Lambda_{\text{UV}}}{M_Z} \simeq 2\times 10^{-6} C_{tt}^{A} \log\frac{\Lambda_{\text{UV}}}{M_Z}. \label{eq:FVRGE1}
\end{equation}

The above leading-log expression for $C_{ds}^{A,V}(M_Z)$ is of only limited phenomenological importance. Due to the numerically large CKM suppression the current experiments probe the loop induced $C_{ds}^{A,V}(M_Z)$ only for values of $f_a$ relatively close to the weak scale. The finite loop corrections are therefore of similar size as the  leading log term in Eq. (<ref>); note that $\Lambda_{\text{UV}}$ can be at most $4\pi f_a$. The complete expression for $C_{ds}^{V}(M_Z)$, including finite corrections, is given in Eq. (<ref>).

§.§.§ ALP decays

To be produced in rare kaon decays, the ALP needs to be lighter than $m_K-m_\pi \approx 350$ MeV. This has phenomenological implication for the ALP total width, possible decay channels and parameter space in general. The most important decay channels are $a\to \gamma\gamma$, and, if kinematically allowed, the decay modes $a\to \ell^+\ell^-$, $\ell=e,\mu$. The corresponding decay widths are given by

= α^2  m_a^3/256π^3 f_a^2 |C_γγ^eff|^2  ,
=   m_a  m_ℓ^2/8πf_a^2 (C_ℓℓ^A )^2 (1-4m_ℓ^2m_a^2)^1/2  ,

where the effective coupling of the ALP to photons is <cit.>

= N_1+N_2 +∑_q  6  Q_q^2 C_qq^A(μ_0) B_1(τ_q)
+ 2∑_ℓ C_ℓℓ^A B_1(τ_ℓ)

- (1.92±0.04) N_3 - m_a^2/(m_π^2-m_a^2)
[ N_3 m_d-m_u/m_d+m_u + C_uu^A-C_dd^A ]  .

The loop function $B_1(\tau_f)$ depends on $\tau_f\equiv 4m_f^2/m_a^2$, and is $B_1\simeq 1$ for light fermions ($m_f\ll m_a$) while it decouples as $B_1\simeq -m_a^2/(12 m_f^2)$ for heavy fermions ($m_f\gg m_a$). The explicit form of $B_1(\tau)$ can be found in Ref. <cit.>.
We neglected the threshold corrections coming from $W$-loops since these are negligible for $m_a\ll m_W$. The second line in Eq. (<ref>) encodes the leading order chiral perturbation theory (ChPT) contribution of $N_3$ to the ALP di-photon coupling. Similarly, the ALP coupling to electrons is generated at one-loop from the ALP couplings to gauge bosons <cit.>.  A common notation in axion literature is also $g_{\gamma\gamma}=C_{\gamma\gamma}^\text{eff} \alpha /(2\pi f_a)$.

For $m_a<2 m_e$, ALP decays are dominated by $a\to \gamma\gamma$, and the decay length is given by

\begin{align}
\approx 2.9\times10^7\,{\rm meters} \left( \frac{f_a}{\rm TeV} \right)^2
\left(\frac{1\,{\rm MeV}}{m_a} \right)^3
\left(\frac{1}{|C_{\gamma\gamma}^\text{eff}|}\right)^2  \, .
\label{eq:ALPlifetime:photon}
\end{align}

For heavier ALPs, the decay channel to lepton pairs opens up and typically dominates over the diphoton decay mode, as long as $C^A_{\ell\ell}\gtrsim  C_{\gamma\gamma}^\text{eff} m_a/{\rm GeV}$. We can then estimate the ALP decay length as

\begin{equation}
\label{eq:ctau_a:leptons}
\approx
\left\{
\begin{aligned}
& 1.9\text{ meters}\left( \frac{f_a}{\rm TeV} \right)^2
\left(\frac{10\,{\rm MeV}}{m_a} \right)
\left(\frac{1}{C_{ee}^A}\right)^2\ , & 2 m_e<m_a< 2 m_\mu,
\\
& 1.5\times 10^{-6}\text{ meters}\left( \frac{f_a}{\rm TeV} \right)^2
\left(\frac{300\,{\rm MeV}}{m_a} \right)
\left(\frac{1}{C_{\mu\mu}^A}\right)^2\ ,& 2 m_\mu<m_a.
\end{aligned}
\right.
\end{equation}

The above scalings show that an ALP lighter than the di-electron threshold is unlikely to decay promptly on detector scales, while a heavier ALP can be  decaying promptly for $f_a\lesssim \text{TeV}$.

The hadronic ALP decay modes are typically negligible. The $a\to 3\pi$ decay is kinematically forbidden, if rare kaon decays with an ALP in the final state are kinematically allowed. The $a\to \pi\pi\gamma$ decay could be important in scenarios where the ALP mixes with $\eta'$, due to ${\cal B}(\eta'\to\pi\pi\gamma)\approx 30\%$. In practice, however, the $a-\eta'$ mixing is suppressed by $(m_a/m_{\eta'})^2$, and for $m_a<m_K-m_\pi$ the $a\to \pi\pi\gamma$ decay channel is always subdominant compared to $a\to\gamma\gamma$ <cit.>. The remaining kinematically allowed decay modes are $a\to\mu^\pm e^\mp$, which is lepton flavor violating (Section <ref>), and $a\to\pi^0 e^+e^-$, which is CP violating and hence expected to be highly suppressed.

§.§.§ ALP production in kaon decays

The decay amplitudes for charged and neutral kaon decays into a pion and an ALP receive an IR contribution at energies $\mu<2$ GeV from the matching of the ALP to the ChPT  Lagrangian <cit.>, and UV contributions which are encoded in $C_{ds}^V$. The latter can be induced at tree level <cit.> or at one-loop via ALP couplings to the up-quarks <cit.> or to the $W$-boson <cit.>. The value of $C_{ds}^V$ at the charm threshold, $\mu_c=2$ GeV, is therefore given by

\begin{equation}
\label{eq:CdsV}
C_{ds}^V(\mu_c)=C_{ds}^V(\Lambda_{\text{UV}}) + \sum_{F=t,c}\frac{y_F^2V_{Fd}^*V_{Fs}C_{FF}^{A}}{16\pi^2}\left(\log\frac{\Lambda_{\text{UV}}}{\mu_F}-f_F(x_F)\right)-\frac{g_2^4N_2}{256\pi^4}V_{td}^*V_{ts}f_W(x_t),
\end{equation}

where $\mu_t=M_Z$, while $x_F=m_F^2/m_W^2$ for $F=t,c$. The loop functions are
f_F(x) = 3(1-x+logx)/2(1-x)^2+1/4,
f_W(x) = 3 x[1-x+xlogx]/2(1-x)^2.
All the Wilson coefficients on the r.h.s. of Eq. (<ref>) are defined at the scale $\Lambda_{\rm UV}$. The log enhanced term in Eq. (<ref>) is the RGE contribution already discussed in Eq. (<ref>), where we used the fact that no resummation of electroweak logs is needed. The other two contributions are the finite threshold corrections from the top, charm and the $W$ boson running in the loop, with the ALP either coupling to the quark or the $W$, respectively. In concrete UV completions, the UV threshold corrections might also be important <cit.>. For simplicity we assume these to be subdominant. Note that only the diagonal axial couplings to quarks contribute to the one loop correction of $C_{ds}^V$.

Numerically, the sizes of the different contributions in Eq. (<ref>) are given by (dropping the trivial UV contribution)

\begin{equation}
-10^{-6}\! \left[(2+i)C_{tt}^A\!\left(\log\frac{\Lambda_{\text{UV}}}{\mu_t} +0.02
\right)\!+\! 0.08 \, C_{cc}^A\!\left(\log\frac{\Lambda_{\text{UV}}}{\mu_c}+11\right)\!-(4+2i)10^{-3} N_2\right]. \label{eq:EWestimate}
\end{equation}

If the UV ALP couplings are flavor diagonal, the FV coupling $C_{ds}^V$, relevant for the kaon decays, is highly suppressed by the product of loop and CKM factors. This means that in the case of UV flavor universality, the $K\to\pi a$ decays probe significantly lower values of $f_a$ than for the case of flavor violating UV ALP couplings. Interestingly, the Yukawa suppression of the charm contribution relative to the top gets partially compensated by the less severe CKM suppression, so that the two terms can be of similar numerical size.

Following Refs. <cit.>, we can write down the $K\to\pi a$ decay amplitudes. The decay widths are given by
$\Gamma(K\to\pi a) = \left|\mathcal{A}(K\to\pi a) \right|^2/ 16\pi m_K $, ignoring the phase space factor due to nonzero pion and ALP masses. Furthermore, for $K_L$ decays the $K^0$ and $\bar K^0$ decay amplitudes interfere, since $K_L = \big[(1+\epsilon) K^0 + (1-\epsilon)\bar K^0\big]/\big[2(1+|\epsilon|^2)\big]^{1/2}$.
At leading order in the ChPT expansion, the decay amplitudes are given by

\begin{align}
\label{eqn:chargedKtopia}
\begin{split}
\!\! i{\cal A}(K^-\to\pi^- a)
&\ \ = \frac{N_8}{4f_a}\,\Big[
8\, N_{3}\,m_{K-\pi}^2 \xi_a  + \big(4\, C_{ss}^A+6 \xi_a C_{uu+dd-2ss}^A\big) \,m_a^2
\\
& \quad
+ C_{2uu+dd+ss}^A\,m_{K-\pi-a}^2  +C_{dd-ss}^V\,m_{K+\pi-a}^2 \Big]  - \frac{m_{K-\pi}^2}{2 f_a} C_{ds}^V \,,
\end{split}
\\
\begin{split}
\label{eqn:neutralKtopia}
\!\! - i\sqrt2{\cal A}(\bar K^0\to\pi^0 a)
& = \frac{N_8}{4f_a}\,\Big[
\big( 8\, N_{3} \xi_a + C_{3dd+ss}^A\big)\,m_{K-\pi}^2
+ \big( C_{2uu-dd-ss}^A - 2 \xi_a C_{uu+dd-2ss}^A\big)m_a^2
\\
&   - 2 C_{uu-dd}^Am_a^2 \,\frac{\,m_{K-a}^2}{m_{\pi-a}^2}
+ C_{dd-ss}^V m_{K+\pi-a}^2 \Big]  - \frac{m_{K-\pi}^2}{2 f_a} C_{ds}^V ,
\end{split}
\end{align}

while the $K^0\to\pi^0a$ decay amplitude is obtained via the replacement  ${\cal A}(K^0\to\pi^0 a) = - {\cal A}(\bar K^0\to\pi^0 a) \big|_{C_{ds}^V \to C_{ds}^{V*} = C_{sd}^V}$. Above, we use the shorthand notation, $C_{c_1 uu+c_2 dd+c_3 ss}^{V}=c_1 C_{uu}^V +c_2 C_{dd}^{V}+c_3 C_{ss}^{V} $,  and similarly for $C_i^A$, for instance, $C_{2uu-dd-ss}^{A}=2 C_{uu}^A- C_{dd}^{A}-C_{ss}^{A} $, with all the Wilson coefficients defined at $\mu_c=2$ GeV. We also defined $m_{M_1-M_2}^2=m_{M_1}^2-m_{M_2}^2$, $m_{K\pm \pi-a}^2=m_K^2\pm m_\pi^2-m_a^2$, and $\xi_a= {(m_K^2-m_a^2)}/{(4m_K^2-m_\pi^2-3m_a^2)}$.
The overall  prefactor multiplying the first terms in Eqs. (<ref>), (<ref>), is given by

\begin{equation}
N_8 = -{G_F} V_{ud}^* V_{us} g_8 f_\pi^2/\sqrt{2}\simeq - 1.5 \times 10^{-7}\ .\end{equation}

If the $K\to \pi a$ transitions are induced by the flavor diagonal couplings of ALP to quarks, with no large hierarchies between the couplings, then the ALP couplings to light quarks and gluons (the terms proportional to $N_8$) in Eqs. (<ref>), (<ref>) and the contributions generated via EW loops from ALP coupling to top, charm and $W$, Eq. (<ref>), are of roughly similar sizes. A more quantitative comparison is shown in Fig. <ref>, where the branching ratios for $K\to \pi a$ decays are plotted as functions of the ALP mass for a number of benchmark scenarios, in each of which only one ALP coupling at a time is taken to be nonzero, as indicated.

For a massless ALP we can write the $K^+\to \pi^+ a$ branching ratio as

\begin{equation}
\label{eq:K+:pi+:ALP:numerical}
{\cal B}(K^+\to \pi^+ a)=1.1\times 10^7 \left(\frac{\text{TeV}}{f_a}\right)^2\big[(C_{ds}^V)^2+2 N_8C_{ds}^V \tilde N_3+ N_8^2 \tilde N_3^2\big],
\end{equation}

where we neglected $m_\pi$ and defined $\tilde N_3 =N_3+\big(C_{2uu+dd+ss}^{A}+C_{dd-ss}^V\big)/2$. The form of Eq. (<ref>) motivates us to distinguish two limiting cases.

* If $C_{ds}^V(\Lambda_{\rm UV})$ is nonzero this implies strong limits on the symmetry breaking scale $f_a$. For $C_{ds}^V\sim {\mathcal O}(1)$ the value of $f_a\sim 10^7~\text{GeV}$ would imply very large modifications of kaon decay branching ratios, ${\cal B} (K\to \pi a)\sim \mathcal{O}(1)$. As discussed in Section <ref>, much larger decay constants $f_a$ are already probed by the experiments, thanks to the fact that the irreducible SM background of $K^+\to\pi^+ a$ is $K^+\to\pi^+\nu\bar\nu$ which is highly suppressed in the SM. In this scenario, the ALP is naturally long-lived and, if light enough, can be the Dark Matter.

* If $C_{ds}^V(\Lambda_{\rm UV})=0$ then ${\cal B}(K^+\to \pi^+ a)\sim 10^{-6}-10^{-7}$ for $f_a\sim \text{TeV}$. In this scenario the ALP can decay promptly into the SM particles, a possibility discussed further in Section <ref>. Intriguingly, for $m_a<m_K-m_\pi$ the region where the ALP decays promptly can be fully probed in the near future.

Branching ratios of the $K^+\to\pi^+a$ (left) and $K_L\to\pi^0a$ (right) decays as  functions of the ALP mass, taking only a single coupling in the ALP Lagrangian (<ref>), (<ref>) to be nonzero as shown by the labels, and setting $f_a=1$ TeV and $\Lambda_{\rm UV}= 4\pi f_a$. For a different value of the ALP decay constant, the branching ratios should be rescaled by $1/f_a^2$.

§.§.§ Long-lived ALPs from kaon decays

Flavor-violating ALP decays obey the generalized Grossman-Nir bound <cit.>

\begin{equation}
\label{eq:2body:GN}
{\cal B}(K_L\to\pi^0 a)\leq4.3\cdot{\cal B}(K^+\to\pi^+a).
\end{equation}

The main effect in the numerical factor is due to the kaon lifetime ratio $\tau_L/\tau_+=4.1$. This makes the $K^+\to\pi^+ a$ decay more sensitive to invisible ALPs than the $K_L\to\pi^0 a$ decay (which is also more challenging experimentally). In the remainder of this section we therefore focus on ALP production in the $K^+\to\pi^+a$ decays. Note however that the above rule does not apply generally. Certain new physics models feature an enhanced signal in rare $K_L$ and not $K^+$ decays, as discussed in Section <ref>.

Experimental status and projections

The experimental status and prospects of $K\to\pi X$ searches, with $X$ escaping the detector or decaying invisibly, are discussed in Section <ref>. For the $K^+$ case, the upper limits of ${\cal B}(K^+\to\pi^+X)$ at 90% CL as a function of $m_X$ up to the kinematic endpoint $m_{K}-m_{\pi} = 354$ MeV are shown in Fig. <ref>. Three distinct kinematic regions are determined by the experimental conditions.

* The regions 1 and 2 of the $K^+\to\pi^+\nu\bar\nu$ search (results shown as red shaded areas). The largest background comes from the ultra-rare $K^+\to\pi^+\nu\bar\nu$ decay, leading to very stringent bounds on ${\cal B}(K^+\to\pi^+a)$. The deterioration of the bound at very low $m_a$ masses owes to the definition of the $K^+\to\pi^+\nu\bar\nu$ signal region.

* The $K^+\to\pi^+\pi^0$ region, where the results are obtained from a dedicated search for the $\pi^0$ decays into a fully invisible final state (results shown as a blue shaded area).

* The region $m_{a}>260$ MeV, where the background is dominated by the $K^+\to\pi^+\pi^+\pi^-$ decays. This region serves as a control region for the $K^+\to\pi^+\nu\bar\nu$ measurement (and thus is not masked during analysis), and has not yet been  explored by the experiments.

Constraints on ${\cal B}(K^+\to\pi^+X)$ with the $X$ particle decaying invisibly or escaping the detector as a function of $X$ mass up to the kinematic endpoint (vertical dotted line). Red regions are excluded by the NA62 search for $K^+\to\pi^+X$ <cit.>, the blue region by $\pi^0\to{\rm inv}$ <cit.>, and grey regions by the BNL E949 experiment <cit.>.

The constraints in Fig. <ref> can be re-interpreted as bounds on the ALP couplings entering the $K^+\to \pi^+ a$ decay amplitude, see Eq. (<ref>). This is usually done either in specific UV completions or by assuming that only a single coupling is nonzero (a conventional approach for obtaining bounds on couplings in the effective field theory).
As illustrative benchmarks we consider two limiting cases:

* Flavor Anarchy: The FV couplings are of the same order as the diagonal ones. In numerical estimates we set at the UV scale $C_{sd}(\Lambda_{\text{UV}})=1$, and similarly for the diagonal couplings, $C^V_{qq} = C^A_{qq} = N_3 = 1$.

* Minimal Flavor Violation: The FV couplings vanish at the UV scale, $C_{sd}(\Lambda_{\text{UV}})=0$. For flavor diagonal couplings we take $C^V_{qq} = C^A_{qq} = N_3 = 1$ at $\mu=\Lambda_{\rm UV}$.

The resulting constraints on the axion decay constant $f_a$ from the NA62 searches are shown for the flavor anarchic and minimal flavor violating benchmarks in Fig. <ref>, overlaid with the SN1987A constraints. In the case of nonzero flavor-violating couplings the most stringent constraints arise mainly from Supernova (SN) cooling through hyperon decays <cit.>, while for the flavor-conserving benchmark the main cooling mechanism is via nucleon bremsstrahlung <cit.>. For the flavor anarchic benchmark the SN cooling constraints are roughly comparable with the constraints from $K-\bar K$ mixing due to tree level ALP exchanges <cit.>. These constraints are UV sensitive, since UV dependent dimension-6 operators suppressed by the UV scale $\Lambda_{
\text{UV}}\sim 4\pi f_a$ are expected to give contributions of the same parametric size as the IR contribution from ALP exchanges. The constraints from $K-\bar K$ mixing in Fig. <ref> should thus be viewed with this caveat in mind.

Fig. <ref> (left) illustrates the fact that for generic ALP couplings, with no peculiar flavor structure, the strongest constraints on the axion decay constant are provided by NA62 and are $f_a\gtrsim 10^{12}~\GeV$ for ALP couplings of ${\mathcal O}(1)$. Remarkably, the NA62 constraints are the leading ones in the phenomenologically very interesting region where $m_a\lesssim {\mathcal O}(\text{keV})$, and the ALP can be cold DM produced through the misalignment mechanism. Most importantly, the ALP lifetime in this region, $\tau_a\gtrsim~10^{26}~\textrm{s}$, is sufficiently long to satisfy the stringent bounds on decaying DM. In this context, see also the more detailed discussion of the QCD axion in Section <ref>.

For the minimal flavor violation benchmark (Fig. <ref>, right), the bounds on $f_a$ are a factor of ${\cal O}(10^{-5})$ less stringent. This is a direct consequence of the CKM suppression required for the flavor violating transition, cf. Eq. (<ref>). In this case the constraints from $K^+\to\pi^+a$ decays are relevant exactly in the SN trapping regime $(f_a \lesssim 10^6~\GeV)$ where SN cooling through ALP emissions is no longer possible, since ALPs get trapped inside the proto-neutron star. At low masses ($m_a\lesssim 1$ MeV) this region is excluded by cosmological constraints coming from CMB and BBN <cit.>. At higher masses, a detailed study of the BBN constraints for the ALP couplings considered here is still lacking. However, as an indication of roughly where the BBN constraint are expected to be applicable we delineate in Fig. <ref> with the solid line the region below which the ALP lifetime is $\tau<1$ s. This as a plausible guess of the parameter space that would survive the BBN constraints without invoking extremely small reheating temperatures or other highly non-standard cosmologies. See Ref. <cit.> for a careful assessment of the dependence of BBN constraints on the cosmological history.

Constraints on the axion decay constant $f_a$ for the flavor anarchic (left) and minimal flavor violation (right) benchmarks. The constraints from NA62 are showed as colored regions, cf. Fig. <ref>. The SN1987A constraints (grey region) are taken from Ref. <cit.> for SN cooling through nucleon bremsstrahlung $NN \to NN +{\rm inv}$, and from Ref. <cit.> for cooling via hyperon decays $\Lambda\to n +{\rm inv}$. The $K-\bar K$ mixing constraints assume that the dimension-6 UV contributions are small. Also shown in the right panel are the contour lines below which the ALP is no longer long-lived ($c \tau_a < 1$ m) and above which its lifetime is stringently constrained by cosmology ($\tau_a > 1$ s).

Consequences for the QCD axion

The prototypical example of a light ALP is the QCD axion addressing the strong CP problem in the SM. As a consequence, the mass of the QCD axion predominantly arises from its couplings to SM gluons even though exceptions are possible in more elaborate models. In minimal scenarios the QCD axion mass can be directly related to the axion decay constant <cit.> (redefining $f_a\to N_3 f_a$ one recovers the standard QCD axion notation):

\begin{align}
m_a = 5.7~\mu{\rm eV} \left( N_3\frac{10^{12}~\GeV}{f_a} \right).
\end{align}

Phenomenological constraints in explicit models usually imply $f_a \gg 10^6~\GeV$, so that the QCD axion is essentially massless and stable for collider purposes. By definition the QCD axion couples to gluons, $N_3 \ne 0$, while couplings to fermions are model dependent. Flavor-violating couplings generically arise when Peccei-Quinn (PQ) charges are flavor non-universal. This possibility is especially well motivated if the PQ symmetry is (partially) responsible for explaining the  hierarchies in the SM Yukawa couplings <cit.>.

In Fig. <ref> we compare the most stringent constraints on the QCD axion couplings at $\mu = \mu_c$. As illustrative examples we consider three distinct cases: i) only the FV vector coupling $C^V_{sd}$ is nonzero, ii) only the FV axial-vector coupling  $C^A_{sd}$ is nonzero, and iii) the combination of flavor-diagonal couplings $C_{qq} \equiv N_3 = C^V_{qq} = C^A_{qq}$ are all nonzero. For all three cases we show  the constraints both from SN cooling (“SN1987A”) as well as the constraints from the laboratory experiments (“Lab”). The SN cooling constraints are due to hyperon decays (for flavor-violating cases) and from nucleon bremsstrahlung (for the flavor-diagonal scenario). For nucleon bremsstrahlung the constraints come from observations of the SN neutrino flux, implying that proto-neutron star (PNS) did not have effective non-standard cooling mechanisms. The bounds from SN cooling apply for values above $f_a \gtrsim 15~{\rm TeV}$, since otherwise axions become trapped inside the PNS <cit.>.
Note that the constraints from axion production through hyperon decays do not have a trapping regime <cit.> and therefore are valid also in the large coupling (small $f_a$) regime. The laboratory constraints are due to the bounds on $K^+\to\pi^+X_{\rm inv}$ for the $C^V_{sd}$ and $C_{qq}$ cases, and from $K^+\to\pi^+\pi^0 X_{\rm inv}$ for the $C^A_{sd}$ case <cit.>, using the results of the ISTRA+ experiment <cit.>. The stronger constraints on $K_L\to\pi^0\pi^0 X_{\rm inv}$ decays from E391a <cit.> do not apply to the $m_a=0$ case.

Fig. <ref> shows that NA62 searches are sensitive to the QCD axion (or more generally, a very light invisible ALP) if it has sizable flavor-violating vector couplings. For flavor-violating axial couplings or the flavor-diagonal case the constraints from SN1987A are instead stronger than the NA62 reach by several orders of magnitude. This model-independent observation allows to test the way Peccei-Quinn $U(1)$ symmetry is embedded in the SM flavor group. In particular, if this abelian symmetry is identified with the symmetry explaining the Yukawa textures à la Froggatt-Nielsen then the resulting QCD axion could be probed simultaneously as DM by the cavity experiments such as ADMX and  via the $K^+\to \pi^+ a$ decays at NA62 <cit.>. If  $U(1)_{\rm PQ}$ is part of a non-abelian flavor group such as $U(2)$ <cit.> the flavor-violating axion couplings can instead be strongly suppressed <cit.>.

The NA62 sensitivity in the ($m_a$, $g_{a\gamma\gamma}$) plane obtained setting $C^V_{sd} = 1$ is shown in Fig. <ref>. In this scenario NA62 can provide constraints on the axion decay constant $f_a$ of the order of $10^{12}$ GeV, which is about two orders of magnitude more stringent than the astrophysical bounds from white dwarf (WD) cooling (for $C^A_{ee} = 1$) and SN1987 (for $C_N = 1$), and about four orders of magnitude stronger than the constraints from CAST/HB stars (for $N_1 = N_2 = 0$). Therefore NA62 can potentially constrain the region in the parameter space for which the QCD axion yields the observed DM abundance via misalignment in the pre-inflationary scenarios, complementing the standard axion searches with haloscopes.

Laboratory and astrophysiccal constraints on the inverse axion couplings for the case of a massless ALP, $m_a = 0$. These  apply unchanged also to the standard QCD axion. The SN1987A bounds are taken from Ref. <cit.> for the case of $C^V_{sd}$ and $C^A_{sd}$ (cooling via hyperon decays) and from Ref. <cit.> for the case of $C_{qq}$ (cooling via nucleon bremsstrahlung).
In contrast to the bounds from hyperon decays the nucleon bremsstrahlung has a trapping regime for smaller values of $f_a$. These values are excluded by Kamiokande constraints <cit.>,
leaving an open parameter space only for very small values of $f_a$, outside the plotted range. The laboratory bounds for $C^V_{sd}$ and $C_{qq}$ ($K \to \pi$ decays) are discussed in the main text,
while the bounds on $C^A_{sd}$ ($K \to\pi\pi$ decays) <cit.> are from the ISTRA+ experiment <cit.>.

Constraints on the QCD axion in the $m_a(f_a)$–$g_{a \gamma \gamma}$ plane. Shown in gray are the usual constraints from existing and planned halo- and helioscopes and from HB star cooling, overlaid is the standard QCD axion band in yellow. In addition, we show the constraints on $m_a/f_a$ for various couplings of order unity; displayed in blue are the HB/CAST bounds <cit.> for $N_1 = N_2 = 0$ (i.e. the “hadronic axion” which couples to photons through its mixing with $\pi^0$), the bounds from star cooling via nucleon bremsstrahlung in SN1987A for $C_N =1$ <cit.>, and via electron bremsstrahlung in white dwarfs (WD) for $C^A_{ee} =1$ <cit.>. Red solid (dashed) lines show the present (future) constraints from NA62 for $C^V_{sd} =1$. In orange we indicate the values of $m_a$ (and $f_a$) for which the QCD axion can fully account for the observed dark matter abundance, depending on whether the PQ  symmetry is broken before or after inflation (taken from Ref. <cit.>).

Probing axial flavor violating ALP couplings

Three-body decays $K^+\to\pi^+\pi^0 a$ and $K_L\to\pi^0\pi^0 a$ are crucial laboratory probes of the axial FV coupling of the ALP. Assuming for simplicity $m_a\ll m_\pi$, the decay widths are given by <cit.>

\begin{align}
\label{eq:Gamma:3body:K+}
& \frac{\text{d}\Gamma(K^+\to \pi^+\pi^0 a)}{\text{d}s}=\frac{(C_{sd}^A)^2}{f_a^2}\frac{(m_{K^+}^2-s)^3}{6144\pi^3m_{K^+}^5}\beta(F_p^2+\beta^2G_p^2+2\beta F_p G_p),\\
\label{eq:Gamma:3body:KL}
& \frac{\text{d}\Gamma(K_L\to \pi^0\pi^0 a)}{\text{d}s}=\frac{(C_{sd}^A)^2}{f_a^2}\frac{(m_{K_L}^2-s)^3}{4096\pi^3m_{K_L}^5}\beta F_s^2,
\end{align}

where $s=(p_{\pi_1}+p_{\pi_2})^2$, $\beta=1-4m_\pi^2/s$, and the functions $F_s$, $F_p$ and $G_p$ can be extracted from the $K^+\to\pi^+\pi^-e^+\nu$ decay form factors. The decay widths satisfy a generalized Grossman-Nir bound

\begin{equation}
{\cal B}(K_L\to\pi^0\pi^0 a)< 31\cdot {\cal B}(K^+\to\pi^+\pi^0 a),
\end{equation}

where the numerical factor on the r.h.s. is enhanced compared to the two-body generalized Grossman-Nir bound, Eq. (<ref>), by phase space effects and form factors.

The experimental status is summarized in Section <ref>. The strongest limit on the massless ALP comes from the ISTRA+ experiment <cit.>:

\begin{equation}
\label{eq:K+:3body}
{\cal B}(K^-\to\pi^-\pi^0 a)<0.9\times 10^{-5}\quad \Rightarrow\quad 2f_a/C^A_{sd}>3.6\times 10^7\text{ GeV}.
\end{equation}

This is expected to be improved with the NA62 Run 1 dataset. The projected NA62 reach would still be below the nominal bounds from the SN cooling due to hyperon decays <cit.> (Fig. <ref>).
However the SN cooling bounds carry an important caveat that they
are likely affected by large systematic uncertainties related to the SN modeling (it is even possible to argue that the bound may be completely absent <cit.>).
In addition to the above bound on the $K^-$ decay, there is a phenomenologically interesting experimental bound on ${\cal B}(K_L\to \pi^0\pi^0 a)$ from the E391a experiment <cit.>, which is set to be improved by KOTO. This bound does not cover the massless ALP (see Section <ref>), while it leads to $2 f_a/C_{sd}^A\gtrsim 1.7\times 10^7$ GeV for a light, yet still massive, ALP <cit.>.

For massive ALP, the expressions for the decay widths in Eqs. (<ref>), (<ref>) need to be modified, and include a new form factor that cannot be extracted from data (see Appendix C in Ref. <cit.>). While a prediction for this form factor can be made using chiral perturbation theory <cit.>, understanding of theoretical errors for the massive ALP case could benefit from further studies.
In summary, experimental exploration of the axial flavor violating coupling of the ALP calls for updates on the experimental and on the theoretical fronts, both of which can be achieved  in the near future.

§.§.§ Prompt ALPs from kaon decays

Experimentally distinct phenomenology is obtained if ALPs decay promptly within the  instrumented volume of a kaon experiment. This requires relatively low values of the ALP decay constant $f_a$.
For $m_a$ in the MeV to GeV range the ALP decay constant  should be within a few orders of magnitude of the EW scale,  $f_a\sim{\mathcal O}(100~\text{GeV})$. In this range of couplings the ALP is in a thermal equilibrium in the early Universe, and its mass is bounded to be heavier than roughly 1 MeV by BBN constraints on the number of relativistic degrees of freedom <cit.>.

Depending on the specific model, different ALP decays are allowed. For light ALPs, the $a\to\gamma\gamma$ or $a\to e^+e^-$ decays dominate, with at most one loop factor hierarchy between the electron and the photon couplings. At one loop, the ALP-electron and ALP-photon couplings mix, cf. Eq. (<ref>). If kinematically allowed, the $a\to\mu^+\mu^-$ decay should also be considered.

Promptly decaying ALPs can be produced in kaon decays via their couplings to gluons, $W$-bosons, or light quarks. ALP flavor violating couplings can be safely ignored in this context because FCNC constraints from $K-\bar K$ mixing require $f_a\gtrsim$ $10^{7}$ GeV and hence imply an invisible ALP on detector scales (see Section <ref> for further discussion). Flavor diagonal couplings with a small enough decay constant would change appreciably the total kaon decay width.
Summing the measured $K^+$ exclusive branching ratios gives a model-independent bound at 95% CL[The PDG performs a constrained fit to $K^+$ branching ratios that includes the total decay width and also the requirement that the SM channels sum up to one, giving ${\cal B}(K^+\to\text{SM})=(100.00\pm0.15)\%$. Ignoring that a constrained fit was used this would imply  ${\cal B}(K^+\to\text{SM}+\text{new physics})<2.9\times 10^{-3}$ at 95% CL. We expect that repeating the PDG fit, but without imposing the constraint, would give a constraint of a similar size. Until such a fit is performed the constraint in (<ref>) should be viewed as only indicative.]

\begin{equation}
{\cal B}(K^+\to\text{SM}+\text{new physics})\lesssim 3 \times 10^{-3}.
\label{eq:modelind}
\end{equation}

This provides a lower bound on the ALP decay constant that is independent of the ALP decay mode, and only depends on the production mechanism. For instance, for the ALP that couples only to gluons this is denoted with a dark green shaded exclusion region in Fig. <ref> (left).

As we will see below, dedicated experimental searches on exclusive channels already exclude significant parts of the prompt ALP parameter space. However, certain regions of parameter space still remain unexplored and deserve further attention. Typically, the parameter space corresponding to $c\tau_a\gtrsim 10$ cm is constrained by the beam-dump experiments (even though a complete study of the implications of these measurements on the ALP parameter space is still lacking), while for shorter lifetimes the searches for rare $K^+$ and $K_L$ decays become relevant. It is in the latter parameter region that the searches for visible ALP signatures at kaon experiments can make significant improvements, potentially completely covering  the parameter space left unexplored at small decay constant.

Left: Bounds on the ALP coupling only to gluons ($N_3(\Lambda_{\rm UV})\ne 0$ only) scenario. The coupling of ALP  to photons is induced at one loop, cf. Eq. (<ref>), allowing for the $a\to \gamma\gamma$ decays.
The Babar bound was derived assuming $\Lambda_{\rm UV}=1$ TeV.
Right: Bounds on the ALP coupling to only  $SU(2)_L$ gauge bosons ($N_2(\Lambda_{\rm UV})\ne 0$ only)  scenario <cit.>. In this case $a\to \gamma\gamma$ is generated already at tree level. In the parameter space region below the dashed line in both plots, the ALP has proper lifetime longer than a second, and is therefore disfavoured by Big Bang Nucleosynthesis considerations. The red dashed lines on the right-hand plot trace the points at which the ALP has an average lab frame decay length of 100 m and 10 m at NA62. These illustrate why the $K^+\to\pi^+X$ bound where $X$ is invisible in the laboratory does not extend to larger couplings, since the ALP becomes too prompt for the search to be sensitive.

Bounds on flavor-diagonal pseudoscalar quark couplings of ALP to the first generation quarks: coupling to up quarks (left) and down quarks (right). In each case, it is assumed that only the relevant coupling is nonzero in the EFT Lagrangian. The coupling to photons is induced at loop level through Eq. (<ref>), allowing the ALP to decay to two photons. In both plots,  the ALP has proper lifetime longer than a second in the parameter space region below the dashed line,  and is therefore disfavoured by Big Bang Nucleosynthesis.

ALP promptly decaying to $\gamma\gamma$

Light ALPs with couplings to gluons and quarks can be produced in kaon decays,  with the decay amplitudes given in Eqs. (<ref>) and (<ref>). Such an ALP then subsequently decays into photons, electrons or muons. If the ALP couples only to gluons or first-generation quarks at tree level, sizeable couplings to photons are generated at loop level, cf. Eq. (<ref>), such that the ALP decays predominantly to photons. For large enough couplings, for which the ALP decay is prompt, the ALP parameter space can be constrained by the $K\to\pi\gamma \gamma$ measurements (Section <ref>). In contrast, if the couplings are small enough, the ALP can be sufficiently long-lived to escape the detector, in which case $K\to\pi X$ searches with a missing energy signature (Section <ref>) become relevant.

In the following we consider
four distinct limiting cases where a single coupling of ALP to the SM dominates at the scale $\mu=\Lambda_{\rm UV}=4\pi f_a$:

(a) coupling to gluons, $N_3\ne 0$, only <cit.>, shown in Fig. <ref> (left);

(b) coupling to $W^i$-bosons, $N_2\ne 0$ only <cit.>, shown in Fig. <ref> (right);

(c) $C_{uu}^A$ only: tree level coupling to only the first generation up quarks, shown in Fig. <ref> (left);

(d) $C_{dd}^A$ only: tree level coupling to only the first generation down quarks, shown in Fig. <ref> (right).

In each of the above four scenarios the ALP couplings to leptons are generated only through two loop radiative corrections and, if prompt, the ALP decays almost exclusively to photons so that we can take ${\mathcal B}(a\to\gamma\gamma)\simeq 1$. The experimental constraints for the four scenarios in the plane of ALP mass and relevant coupling are shown in Figs. <ref>, <ref>. For large $f_a$ (small ALP couplings), the ALP becomes long-lived, and the bounds on ${\mathcal B}(K\to\pi X_{\rm inv})$ apply (Section <ref>). For smaller $f_a$ (larger ALP couplings), the ALP decays promptly to photons, and thus the bounds on  ${\mathcal B}(K\to\pi\gamma\gamma)$ apply (Section <ref>). If ALP couples to charm or top quarks, the flavor violating $K\to \pi a$ decay is induced at one loop via the one loop generated flavour changing coupling, see Eq. (<ref>) and Fig. <ref>. In this case, a large one-loop coupling to electrons is also induced <cit.>, so that the ALP  decays predominantly to electrons.  Consequently, this scenario is already strongly constrained, as discussed below, in Section <ref>. More exotic scenarios, such as the  lepton flavor violating ALP decays, are discussed in Section <ref>.

In many of the $m_a< m_K$ scenarios discussed here, the dominant constraint at small decay constant $f_a$ is the model independent bound on ${\mathcal B}(K^+\to\text{SM}+\text{ anything})$,  cf. Eq. (<ref>). Complementary bounds from other collider probes such as  LEP <cit.>, CDF <cit.>, GlueX <cit.> and $B$-factories <cit.>[The Babar search for $B\to K a(\to \gamma\gamma)$ excludes the large parameter space for the cases (a) and (b). The bound for (a) $N_3\neq 0$ is derived based on the calculation of <cit.> which assumes $\Lambda_{\rm UV}=1$ TeV, not $4\pi f_a$. The bound will be stronger if $\Lambda_{\rm UV}=4\pi f_a$ because the RGE effect is enhanced. The projection of this channel at Belle II is discussed in Ref. <cit.>.
can also be relevant. For the large decay constant, the ALPs are long-lived, so different bounds associated with the long lifetime should be considered:
i) from beam dump experiments, ii) $K\to \pi +$inv bounds, iii) supernova cooling bounds, iv) BBN bounds.

For an ALP coupled to gluons we present the bound from CHARM <cit.> and NuCal <cit.> derived in Ref. <cit.> by considering only the production mode through meson-ALP mixings.
We expect a similar bound to hold for the first generation up and down quark couplings.
For an ALP coupled to EW gauge bosons we take the results from Ref. <cit.>.

Target branching ratios for $K^+\to\pi^+a(\to \gamma\gamma)$ (left) and $K_L\to \pi^0 a(\to \gamma\gamma)$ (right) decay searches for four different benchmark scenarios: the ALP coupling only to gluons (scenario (a), black dashed line), only to $SU(2)_L$ gauge bosons ((b), blue dashed line),  to up quarks ((c), light blue dashed line) and to down quarks ((d), orange dashed line). In each case the regions below the dashed lines are excluded by the beam-dump constraints. The shaded regions are excluded by prompt searches, as indicated.

In the window of intermediate $f_a$ the ALP decays promptly to photons, and thus the limits on ${\cal B}(K\to\pi\gamma\gamma)$ apply. The experimental status of these searches is discussed in Section <ref>. As shown in Figs. <ref> and <ref>, the visible signature of the ALP in $K\to \pi a (\to \gamma\gamma)$ is essential to close the gap between the model independent bound in Eq. (<ref>) and the invisible searches.
We assume the ALP resonance is narrow for simplicity.[However, the signal is broadened when the ALP decays with a significant displacement. For example, in $K^+\to \pi^+\gamma\gamma$ search, the diphoton invariant mass shifts to a lower value than the ALP mass as $m_{\gamma\gamma}\simeq m_a(1-d/L)$ where $L$ is the distance from the $K^+$ decay point to the ECAL, and $d$ is the displacement of ALP decay. Similar issue was discussed in the context of $B^\pm\to K^\pm a(\to \gamma\gamma)$ <cit.>.]
To set target sensitivities for the future experiments, we evaluate in each scenario the value of the $K\to \pi a(\to\gamma\gamma)$ branching ratio corresponding to the boundary of the region already excluded by the beam dump experiments. These target branching ratios for $K^+$ and $K_L$ decays are shown in Fig. <ref>.
For $m_a\lesssim 100~\text{MeV}$, the promptly decaying ALP that couples predominantly to either gluons or to electroweak gauge bosons (scenarios (a) and (b)) is already excluded by the search at BNL-E949 <cit.>. For $m_a\gtrsim 230$ MeV, an ALP coupling to gluons is excluded by the $K^+\to\pi^+\gamma\gamma$ measurement at NA48/NA62 <cit.>, while a $10^{-10}$ sensitivity in the branching ratio is necessary for a definitive test of the electroweak scenario. Note that the recent Babar analysis excludes the large paramater space in a mass range $175 \lesssim m_a \lesssim 260$ MeV for the gluon  and  electroweak scenarios <cit.>. Similarly, reducing the bounds on ${\cal B}(K_L\to \pi^0\gamma\gamma)$ to $\sim 10^{-10}$ at KOTO or KLEVER would test the scenarios (a) and (b) for $m_a\gtrsim 170$ MeV.

A particularly difficult mass window to cover is the mass region $100~\text{MeV}<m_a<200~\text{MeV}$, where the $K\to\pi a(\to \gamma\gamma)$ searches are affected by the $K\to 2\pi$ background. This mass interval is characterized by large mixing between the ALP and the pion which could be in tension with measurements of the neutral pion lifetime. For this reason a tiny mass window around the neutral pion mass is shaded in the plots. A more detailed theoretical study of this window would be required to properly assess the experimental chances of closing it. Interestingly, $B$ meson decays might close this mass window for an ALP coupled to gluons or electroweak gauge bosons by using the lifetime information to separate the ALP from the pion. Conversely, for an ALP coupled solely to up or down generation quarks the $b\to s$ transition will be highly suppressed and the $100~\text{MeV}<m_a<200~\text{MeV}$ window will have to be closed by kaon experiments.

ALP promptly decaying to di-leptons

Left: Bounds on the coupling of an ALP to gluons for a scenario in which only the couplings to gluons and to leptons are nonzero in the Lagrangian. The coupling to leptons is set by $C^A_{ll}=N_3$ and is lepton universal. The coupling to photons is induced at loop level through equation (<ref>), allowing the ALP to decay to two photons. Right: Bounds on the coupling of an ALP to $SU(2)_L$ gauge bosons for a scenario in which only this coupling and the coupling to leptons are nonzero in the Lagrangian <cit.>. The coupling to leptons is set by $C^A_{ll}=N_2$ and is lepton universal. The ALP decay to photons occurs at tree level through the $N_2$ coupling.

Target branching ratios for $K^+\to\pi^+a(\to e^+e^-)$ (left) and $K_L\to \pi^0 a(\to e^+e^-)$ (right) decay searches for two different benchmark scenarios: the ALP couplings to gauge bosons, either $N_2$ or $N_3$, where the ALP also couples (diagonally and flavour-universally) to leptons. In each case, the regions below the dashed lines are excluded by the beam-dump searches. The shaded regions are excluded by prompt decay searches, as indicated. Due to the displacement of the ALP decay, there is sometimes a discrepancy between the bound assuming prompt decay and the effective bound on each scenario. The dotted lines show the effective bounds for the $N_3=C_{ll}^A$ scenario of $K^+\to\pi^+ e^+e^-$, and the $N_2=C_{ll}^A$ scenario of $K_L\to\pi^0 e^+e^-$.

Measurements of $K^+\to\pi^+e^+e^-$ <cit.> and $K^+\to\pi^+\mu^+\mu^-$ <cit.> decays, dedicated peak searches in the $m_{ee}$ <cit.> and $m_{\mu\mu}$ <cit.> spectra of these decays, and searches for the $K_L\to\pi^0 e^+ e^-$ <cit.> and $K_L\to\pi^0\mu^+\mu^-$ <cit.> decays lead to constraints on the ALP parameter space.
The details and prospects for these measurements are discussed in Section <ref>. To obtain the corresponding constraints on ALP production and decay, we use  the two sigma upper range on the total number of measured events in the measurements and searches of Refs. <cit.>, and use the peak-search results from Refs. <cit.>. The ALP lifetime tends to exceed 10 ps, especially for $m_a<2m_\mu$. To estimate the fraction of ALP decays that pass the event selections we assume that the ALP energy is half of the mean kaon energy, and keep the ALP signal that has a displacement less than $1\,\text{m}$ at KTeV and NA48/2 ($1.4\,\text{cm}$ at BNL E865).

In Fig. <ref> we show two scenarios in which at $\Lambda_{\rm UV}$ in addition to the ALP couplings to gauge bosons, either $N_2$ or $N_3$, the ALP also has a (diagonal and flavour-universal) coupling to leptons, $C_{\ell \ell}^A=N_{2,3}$, large enough to dominate the ALP decay width. In the two scenarios, the ALP branching ratio to photons ($a\to\gamma\gamma$) is reduced compared to Fig. <ref> due to the availability of the leptonic decay channels.

In Fig. <ref> we show the target $K\to\pi a (\to e^+e^-)$ branching ratios for future experiments, required to reach the lower boundary of the parameter space already excluded by the beam dump experiments (while the constraints from kaon invisible searches are weaker).
The sensitivity of the $K\to\pi a(\to e^+e^-)$ searches offered by the currently available samples is potentially sufficient to cover the gap to the region of parameter space in which ALPs have very displaced decays and where ALPs thus appear mainly as missing momentum in the experiments. In previous experimental analyses,  the decay products were assumed to originate from a single prompt decay vertex, which reduces the acceptance to the actual signal events that comprise of displaced $a\to e^+e^-$ decay vertices. The effective bounds for the $N_3=C_{ll}^A$ scenario from $K^+\to\pi^+ e^+e^-$ decays and for the $N_2=C_{ll}^A$ scenario from $K_L\to \pi^0 e^+e^-$ searches (shown as dotted lines in Fig. <ref>, with regions above the lines excluded) take the reduced acceptance into account, and  differ significantly from the bounds that assume prompt decays (shaded regions in Fig. <ref>). The unexcluded parameter space above the beam-dump constraints (dashed lines in Fig. <ref>) corresponds to ALP decays with large vertex displacements. Future searches with rare kaon decays would therefore optimally relax the decay vertex requirement, and perform scans in both ALP mass and lifetime.

Note that beam-dump constraints weaken above the di-muon threshold, thus the $K\to\pi a (\to\mu^+\mu^-)$ searches represent a unique opportunity to cover the unexplored parameter space for $m_a>210$ MeV. In Fig. <ref> we show the present bounds from $B\to K a(\gamma\gamma)$ derived in Ref. <cit.> and we expect searches for $B\to K a(e^+e^-)$ and $B\to K a(\mu^+\mu^-)$ to also have an impact on this parameter space. Currently, we include the $B^0\to K^{0*} a(\to \mu^+\mu^-)$ search at LHCb <cit.>, which is one of the leading constraints for $m_a>2m_\mu$.

§.§.§ A simplified model of prompt ALP decays to $\ell^+\ell^-$ final states

Prompt leptonic ALP decays (with $c\tau_a\lesssim\mathcal{O}(\text{m})$) are only possible for sufficiently low decay constants $f_a/C_{ee}^A\lesssim\mathcal{O}$(TeV), see Eq. (<ref>). Such low $f_a$ values face stringent constraints. Most notably, the flavor off-diagonal couplings $C_{q_i q_j}^A$ ($i\neq j$) must be strongly suppressed to be compatible with bounds on flavor violation, see Fig. <ref>.  In this case, it is natural to express the flavor-diagonal ALP couplings in the Yukawa basis with the ALP appearing as a phase of the fermion masses in the form $e^{i\,Q^{_\text{PQ}}_{f}\,a/f_a}m_f\,ff^c$, where $Q^{_\text{PQ}}_f$
is the PQ charge of the fermion species $f$. Such flavor alignment without flavor universality could be achieved, for instance, in UV completions where only right-handed SM fermions are charged under PQ symmetry. Performing a chiral rotation of the right-handed fermions $f^c$, these couplings can be mapped into the basis in Eqs. (<ref>–<ref>) to yield $C^V_{f_i f_j}=C^A_{f_i f_j}=\delta_{ij}\,Q^{_\text{PQ}}_f$, $N_3 =N_3\big|_{\rm UV}+\sum_{i=1}^6 Q^{_\text{PQ}}_{q_i}$,   $N_2 = N_2\big|_{\rm UV}$,   $N_1 = N_{1}\big|_{\rm UV}-2\,\sum_{i\,=1}^3\big(Q^{_\text{PQ}}_{\ell_i}+\frac{4}{3}\,Q^{_\text{PQ}}_{u_i}+\frac{1}{3}\,Q^{_\text{PQ}}_{d_i}\big)$.

Strong constraints apply to this ALP parameter space for the mass range  $2m_e\lsim m_a \lsim 2m_\mu$.

* A PQ charge for the top quark larger than $\big|Q^{_\text{PQ}}_t\big| \approx f_a/(100\,\text{TeV})$ would lead to experimentally excluded rates for $b\to s a$ and $s\to d a$ transitions via one loop penguin diagrams.

* The PQ charges of charm and bottom quarks, $Q^{_\text{PQ}}_{c}$ and $Q^{_\text{PQ}}_{b}$, are significantly constrained by bounds on axionic decays of $B$, $D$, $\Upsilon$ and $J/\Psi$ mesons, as well as their contributions to ALP-pion mixing.

* ALP-pion mixing is strongly constrained by bounds on ${\cal B}(\pi^{+}\to e^{+}\nu\,(a\to e^{+}e^{-}))$ <cit.> and by the observed value of ${\cal B}(\pi^{0}\to e^{+}e^{-})$ <cit.>, which require the $a-\pi^0$ mixing angle $\big|\theta_{a\pi}\big|\lesssim 10^{-4}$. At LO in ChPT, $\theta_{a\pi}$ is given by:

\begin{equation}
\theta_{a\pi}\,\approx\,
-\frac{(m_u\, Q^{_\text{PQ}}_{u}-m_d\,Q^{_\text{PQ}}_{d})}{(m_u+m_d)}\frac{f_\pi}{f_a}-
\frac{1}{2}\,\bigg(\!Q^{_\text{PQ}}_{s}+N_3\big|_{\rm UV}\!\bigg)\,\frac{(m_u-m_d)}{(m_u+m_d)}\frac{f_\pi}{f_a}\,+\,
\mathcal{O}\bigg(\!\frac{m_{u,d}}{m_s}\!\bigg)\frac{f_\pi}{f_a}\,,
\end{equation}

where $f_\pi=92$ MeV and we assume $m_a\ll m_{\pi^{_0}}$. ALP-pion mixing could then be made compatible with experimental bounds if $Q^{_\text{PQ}}_{u}=2\,Q^{_\text{PQ}}_{d}=\mathcal{O}(1)$ and $N_3\big|_{\rm UV}=Q^{_\text{PQ}}_{s}=0$  <cit.>.

Order-of-magnitude range of ${\cal B}(K^+\to\pi^+a)$ (left purple axes), ${\cal B}(K_L\to\pi^0 a)$ (left cyan axes), ${\cal B}(K^+\to\mu^{+}\nu a)$ (left red axis), and contours of ${\cal B}(K^+\to\pi^+\pi^0 a)$ (orange curves), and ${\cal B}(K_L\to\pi^0\pi^0a)$ (green curves) in the parameter space of ALP mass ($m_a$) vs inverse decay constant ($1/f_a$) for the prompt $a\to e^+e^-$ benchmark scenario (<ref>). Contours of ALP lifetime are shown as dotted black curves. Shaded regions are excluded by the indicated searches. Expected upper limits on ${\cal B}(K^+\to\pi^+(\pi^0)a)$ from the full NA48/2 dataset are shown as dashed dark (light) red curves.
In the left axes we ignore the $\sim 30\%$ (factor of $\sim\!\!3$) variation of ${\cal B}(K\to\pi a)$ (${\cal B}(K^+\to\mu^+\nu a)$) with $m_a$.
For $K\to\pi a$ decays we show the estimated branching ratios with and without the assumption of $\Delta S=1$ octet-enhancement of these amplitudes (see main text and Ref. <cit.>).

The suppression of $\theta_{a\pi}$ combined with poorly determined $\mathcal{O}(1)$ corrections to $a-\eta$ and $a-\eta^{_\prime}$ mixing at second order in the chiral expansion imply that the naive LO ChPT estimation of hadronic ALP couplings is not reliable <cit.>. Therefore, in this subsubsection we consider a simplified model in which we parametrize the low energy hadronic ALP couplings via a generalized Bardeen-Tye current <cit.>, which is defined by the ALP-meson mixing angles $\theta_{a\pi}$, $\theta_{a\eta_{ud}}$, and $\theta_{a\eta_{s}}$.

Axionic kaon decay amplitudes can then be calculated via $\mathcal{A}(K\to aX)=\sum_\phi\theta_{a\phi}\;\mathcal{A}(K\to \phi^* X)\big|_{p^2_{\phi}=m_a^2}$, with $\phi=\pi^0, \eta_{ud}, \eta_s$ (the SM amplitudes can be found in <cit.>).
As illustrative target benchmarks for experimental searches, we take

\begin{equation}\label{paramBench}
\big|Q^{_\text{PQ}}_e\big|=1~,\qquad \big|\theta_{a\pi}\big|\sim\mathcal{O}\left(\frac{50\text{ keV}}{f_a}\right)~,\qquad \big|\thetaUD\big|\sim\big|\thetaS\big|\sim\mathcal{O}\left(\frac{\text{MeV}}{f_a}\right),
\end{equation}

in anticipation that improved ChPT predictions will map these low energy couplings into a nonvanishing region of the UV parameter space.

The resulting branching ratio predictions can still show significant variation due to the possibility of interference between the amplitudes $\mathcal{A}(K\to\phi^*X)$. Therefore in plotting order-of-magnitude estimates for the various axionic kaon decay channels in Fig. <ref> we assume the following branching ratios, consistent with (<ref>):

B(K^+→π^+a)    ≃       (GeV/f_a)^2  |p⃗_a| /|p⃗_π^_+|  ×

3.6×10^-5    (octet-enhanced),

10^-7    (non octet-enhanced),

B(K_L→π^0a)    ≃       (GeV/f_a)^2  |p⃗_a| /|p⃗_π^_0|  ×

10^-9    (octet-enhanced),

10^-11    (non octet-enhanced),

B(K^+→μ^+νa)    ≃        3 ×10^-8 (GeV/f_a)^2 × Φ(m_a)/Φ(m_a=10 MeV),

where $\Phi(m_a)$ is the Dalitz phase-space integral. The predictions for ${\cal B}(K^+\to\pi^+a)$ and ${\cal B}(K_L\to\pi^0a)$ differ in whether $\Delta S=1$ octet-enhancement (a.k.a. $\Delta I=1/2$ enhancement) is assumed to arise at $\mathcal{O}(p^{_{^2}}\!)$ (octet-enhanced) or at $\mathcal{O}(p^{_{^4}}\!)$ (non-octet-enhanced) in ChPT power counting (see <cit.> for details). Fig. <ref> shows branching ratio contours for the three-body decays $K\to\pi\pi a$, where we assume the same relative sign for $\theta_{a\pi}$, $\thetaUD$ and $\thetaS$, and we neglect the contributions of non-octet-enhanced amplitudes.

Sensitivity projections at 90% CL for ${\cal B}(K^+\to\pi^+a)$ and ${\cal B}(K^+\to\pi^+\pi^0a)$ from the NA48/2 analyses of $K^+\to\pi^+e^+e^-$ <cit.> and $K^+\to\pi^+\pi^0e^+e^-$ <cit.> are also shown in Fig. <ref> (see Sections <ref>, <ref>). These projections assume that: (i) the signal acceptance for $K^+\to\pi^+(a\to e^+e^-)$ is the same as in Fig. 4 of Ref.<cit.>; (ii) the signal acceptance for $K^+\to\pi^+\pi^0(a\to e^+e^-)$ is the same as the acceptance for the IB decay shown in Fig. 3 of Ref.<cit.>; and (iii) the expected backgrounds are accounted for using the full background simulations in Refs.<cit.>.
KTeV also placed a bound on the SM decay ${\mathcal B}(K_L\to\pi^0\pi^0e^+e^-)<6.6 \cdot 10^{-9}$ at 90% CL <cit.> (see Section <ref>). However, due to the low signal efficiency resulting from background rejection cuts, recasting this analysis into bounds on ${\mathcal B}(K_L\to \pi^0\pi^0 a)$ would require a dedicated simulation of the ALP signal.

Order-of-magnitude ranges of axionic kaon decays in the parameter space of ALP mass ($m_a$) vs inverse decay constant ($1/f_a$) assuming dominance of the prompt $a\to\mu^+\mu^-$ channel. Shaded regions are excluded. Contours of ALP lifetime are shown as dotted black curves. Contours for ${\cal B}(K^+\to\pi^+a)$, ${\cal B}(K_L\to\pi^0a)$ and ${\cal B}(K^+\to\mu^+\nu a)$ are shown as dashed purple, solid cyan, and solid red curves, respectively (see text for assumptions). For $K^+\to\pi^+a$, we show the estimated branching ratios with and without the assumption of $\Delta S=1$ octet-enhancement of the decay amplitude (see <cit.> for details).

We turn next to the ALP mass range $2\,m_\mu<m_a<m_K-m_\pi$ and consider the case of prompt $a\to\mu^+\mu^-$ decays. If the ALP couples to heavy (top and/or bottom) quarks, constraints from $B\to K(a\to\mu^+\mu^-)$ decays are severe enough to preclude any viable parameter space in the prompt region. Therefore we restrict our discussion to a special region of parameter space for which $Q^{_\text{PQ}}_t=Q^{_\text{PQ}}_b=0$. On the other hand, $\mathcal{O}(1)$ ALP couplings to $G\tilde{G}$ and second generation quarks are still viable as long as constraints on $\theta_{a\pi}$ can be evaded, which can be obtained by setting $Q^{_\text{PQ}}_e=0$. With these assumptions, we have:

\begin{equation}\label{paramBenchmumu}
\big|Q^{_\text{PQ}}_\mu\big|=1~,\qquad \big|\theta_{a\pi}\big| \sim \big|\thetaUD\big|\sim\big|\thetaS\big|\sim\mathcal{O}\left(\frac{f_\pi}{f_a}\right).
\end{equation}

In Fig. <ref> we show existing constraints on the prompt $a\to\mu^+\mu^-$ parameter space derived from NA48/2 measurements (Refs. <cit.>; Section <ref>), as well as $K\to\pi a$ and $K^+\to\mu^+\nu a$ branching ratio contours. The same considerations and caveats for the discussion of $a\to e^+e^-$ apply in this case: the estimated range of branching ratios can vary significantly due to possible constructive or destructive interference between amplitudes. For the contours in Fig. <ref>, we assume the following branching ratios consistent with Eq. (<ref>):

B(K^+→π^+a)    ≃        (f_π/f_a)^2  |p⃗_a| /|p⃗_π^_+|  ×

32.2    (octet-enhanced),

0.072    (non octet-enhanced),

B(K_L→π^0a)    ≃         0.72×10^-3 (f_π/f_a)^2  |p⃗_a| /|p⃗_π^_0|   (octet-enhanced),

B(K^+→μ^+νa)    ≃         3.36 (f_π/f_a)^2 × Φ(m_a)/Φ(m_a=10 MeV),

where $\Phi(m_a)$ in (<ref>) is the $m_a$-dependent integral over the Dalitz phase-space.

§.§ Heavy neutral leptons

Authors: Fox, Hostert,  Kelly

The existence of heavy neutral leptons (HNLs) has been proposed as one possible solution to the smallness of the observed neutrino masses. In these scenarios, the observed neutrino masses are so small due to one of several “seesaw” mechanisms, where a heavier neutrino drives the mass well below that of the charged fermions of the Standard Model. While these heavy neutral leptons can span orders of magnitude in masses, we are interested here in those below the kaon mass to produce observable signatures in kaon factories.

If one or more heavy neutral leptons $N$ exist, they can interact with the other fields of the SM via the Lagrangian

\begin{eqnarray}
- \mathcal{L} \, \supset \,  -y_N \left( L H\right) N\ + \ {\rm h.c.},
\end{eqnarray}

along with a Majorana mass term generating a mass for the heavy neutral lepton. The existence of the Yukawa coupling $y_N$, along with this mass, will generate mixing between the heavy, mostly-sterile state and the light, mostly-active neutrino states. The mixing, in the case where one $N$ is assumed, is parameterized in terms of the elements of an extended, $4\times 4$ leptonic mixing matrix: $\left\lvert U_{eN}\right\rvert^2$, $\left\lvert U_{\mu N}\right\rvert^2$, and $\left\lvert U_{\tau N}\right\rvert^2$ (see, e.g., Refs. <cit.> for constraints on these angles as a function of $N$ mass). In most phenomenological studies, one of these mixing angles is taken to be nonzero at a time, however, combinations of two or three can yield less trivial signatures.

In this minimal “mass plus mixing” scenario, the $N$ can be produced in any process that yields SM neutrinos, as long as its mass is small enough to be kinematically accessible in the process. At kaon factories, this leads to the two most relevant channels for HNL searches <cit.>:

\begin{equation}
K^+\to\ell^+N, \quad \ell=e,\mu.
\end{equation}

The rates of the above decays are proportional to $|U_{\ell N}|^2$, with explicit expressions available e.g. in Ref. <cit.>.
In most of the parameter space of interest, $N$ is long-lived on detector timescales and so the signature is a single charged lepton and missing mass. The HNL lifetime (for $m_N<m_K$) exceeds $10^{-4}/|U_4|^2~\mu$s, where $|U_4|^2$ is the largest of the three mixing parameters $|U_{\ell N}|^2$, $\ell=e,\mu,\tau$ <cit.>.

The HNLs can modify the early Universe evolution, either by contributing excessively to the energy density of the Universe or through decay products. In particular, if HNLs are sufficiently long lived to survive until the time of BBN, the decays to SM particles can modify the precisely-measured abundance of light elements by altering the proton to neutron ratio <cit.>. This leads to an upper bound on the HNL lifetime  $\tau^0 < 0.1$ s (or $\tau^0 < 0.023$ s when the HNL can decay to pions <cit.>), which then sets lower limits on the mixings of HNLs with the SM neutrinos. The BBN limits are complementary to laboratory searches. The combination of the two covers extensively  the parameter space below the kaon mass, with some parameter space  still open in particular due to weak constraints on $|U_{\tau 4}|^2$ <cit.>. It should be noted, however, that the BBN limits are model-dependent, and can be avoided in models where HNLs decay through forces other than the weak force <cit.>. In such models both the decay-in-flight as well as the BBN constraints are modified, depending on the new branching ratios and lifetimes of the HNLs, while the kaon constraints remain largely unchanged (see

Sec. <ref> for some examples).

The NA62 collaboration has performed world-leading searches for HNL production in kaon decays, leading to some of the most stringent limits on $\left\lvert U_{eN}\right\rvert^2$ and $\left\lvert U_{\mu N}\right\rvert^2$ in the $m_N < m_K$ regime <cit.>. These constraints are shown as red excluded regions in Fig. <ref>, compared to other experimental searches for the minimal “mass plus mixing” scenario. The grey bands in each panel correspond to the naïve target for a seesaw solution to the light neutrino masses. These are calculated assuming a simple one-light-one-heavy-neutrino scenario, where the light neutrino acquires a mass equal to $m_\nu = |U_{\alpha 4}|^2 m_N$, where we take $m_\nu = \sqrt{\Delta m_{21}^2}\approx 9\ \mathrm{meV}$ and $\sqrt{|\Delta m_{31}^2|}\approx 50\ \mathrm{meV}$ from neutrino oscillation measurements <cit.> as the minimum mass of two of the three light-neutrino eigenstates in nature. When considering a three-light-two-heavy-neutrino scenario (the minimum required to explain oscillation results in a seesaw mechanism), additional mixing parameters break the simple dependence between light-neutrino masses, heavy-neutrino masses, and mixing, making this target broader. Experimental details and future prospects of the $K^+\to\ell^+N$ measurements are discussed in Section <ref>.

Depending on the HNL mass, the branching ratio of the $K^+\to \pi^0\ell^+ N$ decay can be comparable to that of the $K^+\to\ell^+N$ decay, both of which are proportional to $\left \lvert U_{\ell N}\right\rvert^2$. This is particularly true for the electron-coupled final state with $m_N \lesssim 20$ MeV <cit.>. If experiments can constrain ${\cal B}(K^+\to\pi^0 e^+N)$ at comparable levels as ${\cal B}(K^+\to e^+N$), then smaller $\left \lvert U_{eN}\right\rvert^2$ can be probed. However the currently available datasets only offer a limited sensitivity (Section <ref>).

Experimental constraints on electron-coupled (left) and muon-coupled (right) heavy neutral leptons as functions of their mass and the mixing $|U_{\ell N}|^2$, along with the naïve seesaw target (grey bands).

In the minimal Type-I seesaw model $N$ is a Majorana fermion with a mass term that violates lepton number by two units. After production in kaon decays, the charged-current decays of $N$ would lead to lepton number violating final states <cit.>, such as $K^+\to\ell^+ N$, $N \to \ell^+\pi^-$. Therefore searches for signatures such as $K^+ \to\pi^-(\pi^0)\ell_1^+\ell_2^+$ represent a unique probe of LNV in kaon facilities (Section <ref>).

Beyond the minimal scenarios, there is also the possibility that $N$ has additional interactions, for instance, if it is gauged under a new $U(1)^\prime$ symmetry <cit.>. This can lead to short-lived $N$ which decay within the detector into, for instance, a neutrino and a charged lepton pair:

\begin{equation}
K^+\to\mu^+ N, \quad N\to\nu e^+ e^-.
\end{equation}

With the additional interactions considered, HNL lifetime becomes effectively an additional free parameter. Peak searches in this three-track signature (Section <ref>) could constrain neutrino mixing directly under different HNL lifetime assumptions. More exotic HNL production channels are discussed in Section <ref>.

§.§ Dark photon

Authors: Dobrescu, Fabbrichesi, Gabrielli, Kvedaraitė,  Martin Camalich, Terol-Calvo

The dark photon is a gauge boson of a secluded $U(1)_{\rm d}$ gauge group under which the dark sector states are charged. If the  $U(1)_{\rm d}$ symmetry is exact, the dark photon is massless.  If the symmetry on the other hand is  either spontaneously broken by the Higgs mechanism or explicitly broken, as in the Stueckelberg Lagrangian, the dark photon is massive.

The dark photon interacts with SM particles via kinetic mixing with the photon <cit.>, or via dipole-like operators <cit.>. The most general kinetic part of the $U(1)_{\rm em}\times U(1)_{\rm d}$ Lagrangian is given by
L_0=- 1/4 F_em, μν F^μν_em- 1/4 F_d μν F_d^μν -ε/2 F_d μν F^μν_em   ,
where the term proportional to $\varepsilon$ introduces the mixing. The gauge boson $A^{\mu}_{\rm em}$ couples  to  the ordinary electromagnetic matter current $J_{\mu}$, while  $A_{\rm d}^{\mu}$ couples to the dark-sector matter current $J_{\mu}^\prime $,
L= e   J_μ A^μ_em + e^'J_ μ^'A^μ_d   ,
with $e$ and $e^\prime$ the respective gauge couplings.

For a massless dark photon, the kinetic terms in Eq.( <ref>) can be conveniently diagonalized by a rotation and a field redefinition of the gauge fields such that the  interaction Lagrangian  takes the following form <cit.>
L_int = e^'J^'_μ A^'μ
+ 1/√(1 -ε^2) ( - εe^'J_μ ^'+ e J_μ ) A^μ  .

The redefined dark photon $A^{\prime \mu}$ then does not couple to the SM matter, while $A^{\mu}$ couples to the SM and represents the ordinary photon. As a result, the dark sector particles carry electric charges, which are small (“milli-charges”) for $\varepsilon \ll 1$, being proportional to $\varepsilon e'/e$.

For a massive dark photon, the physical states are obtained through a field redefinition that both diagonalizes the kinetic term and leaves the mass matrix diagonal. This implies that the ordinary photon does not couple to the dark sector, while the dark photon couples to the electrically-charged SM fermions, proportionally to the mixing parameter $\varepsilon$:
ℒ_𝒾𝓃𝓉^' =  -ε  e/ √(1 - ε^2)   J_μ  A^'μ  ≃- ε  e  J_μA^'μ.

Experimental limits on the kinetic mixing  $\varepsilon$ as a function of the dark photon mass for $m_{A'}>1$ MeV  <cit.>.

Although the massless dark photon (also referred to as the “photon-prime" <cit.>) does not interact at tree level with the SM, such interactions can be generated radiatively if the full UV theory contains messengers coupled to the SM and also to states charged under $U(1)_d$. At low energies
the interactions of the dark photon with the SM fermions are described by higher-dimensional operators, starting with the dimension-5 dipole operators <cit.>:

\begin{equation}
\label{dipole}
\mathcal L_{\rm dim-5}= \dfrac{1}{\Lambda}\bar{\psi}_i\sigma^{\mu\nu}\left(\mathbb{D}_M^{ij}+i \,  \mathbb{D}_E^{ij}\gamma_5\right)\psi_j F^{\prime}_{\mu\nu} ~~,
\end{equation}

where $\sigma_{\mu\nu}=i/2 \, [\gamma_{\mu},\gamma_{\nu}]$, $\psi^i$ are the SM fermions with $i=1,2,3$ labeling the generations, and $F^{\prime}_{\mu\nu}$ is the  dark photon field strength.
The effective couplings $\mathbb{D}_M$ ($\mathbb{D}_E$) multiplying the magnetic (electric) dipole moment  depend on the parameters of the underlying UV model and are proportional to $e'$. They are suppressed by the scale $\Lambda$, representing the typical mass of the dark sector states.
Note that for $i\ne j$ the neutral current transitions induced by (<ref>)  are flavor violating. Such flavor violating dipole transitions are possible also for massive dark photon. This is in contrast to the tree-level couplings of the massive dark photon to the SM which are proportional to the electromagnetic SM current, and thus flavor diagonal.
For $s\to d $ transitions the measured kaon mass mixing imposes a strong constraint on the parameters of the UV model and therefore on the scale $\Lambda$ suppressing the dipole operator  (<ref>).

Dark photons can be searched for in various experiments and observational data. The current limits on $\Lambda/\mathbb{D}^{ij}$ for the massless dark photon are shown in Fig. <ref>,
where $\mathbb{D}^{ij}\equiv\big((\mathbb{D}^{ij}_M)^2+(\mathbb{D}^{ij}_E)^2\big)^{1/2}\!$ is the effective dipole coupling to leptons or light quarks.
Even though some of the limits are as stringent as $\Lambda/\mathbb{D}^{ij} \gtrsim 10^{10}$ GeV, the effective dipole coupling may be small enough to allow the new particles charged under $U(1)_{\rm d}$ to have masses at the TeV scale <cit.>. In the next subsection we discuss in details the different probes of dipoles interactions of the massless dark photon.

For the massive  scenario,  the current limits on the photon-dark photon mixing are shown in Fig. <ref>, with $m_{A'}\simeq 1$ MeV delineating the regions where the dark photon decays into either “visible” SM particles (for larger masses). As we see from the figure, rare kaon decays play an important role in constraining the massive dark photon in the region between 10 and 100 MeV. The current best limit are set by the NA48 experiment <cit.>, looking for $\pi^0\to A'\gamma$ with the dark photon $A'$ further decaying into $e^+e^-$ pairs. The pion dataset is collected from identified kaon decays (e.g. $K^{\pm}\to \pi^{\pm}\pi^0$, $K^{\pm}\to \pi^0\mu^+\nu$ and the dark photon yield is controlled by the mixing squared $\sim\varepsilon^2$. Similarly, one could hunt for the dark photon at NA62 in all the channels containing SM photons, performing a bump hunt on the dielectron spectrum of the collected data. For instance one could estimate

\begin{equation}
{\cal B}(K^+\to \mu^+ \nu A' (e^+ e^-))\sim 10^{-3}\times \varepsilon^2,
\end{equation}

which should yield a competitive bound on the dark photon mixing angle.

Experimental limits on the  parameter space of the  massless  dark photon, with $\Lambda$ the new physics scale and $\mathbb{D}$ the effective dipole coupling (adapted from Refs. <cit.>, <cit.>).

§.§.§ Searches for dark photons dipoles in kaon and hyperon decays

Flavor-changing neutral current (FCNC) processes involving the dark photon are induced by off-diagonal magnetic- and electric-like dipole operators in Eq. (<ref>). In particular, the $\Delta S=1$ transitions provide the leading contribution to the gauge-invariant FCNC interactions with SM quarks, for both massless and massive dark photon scenarios.

* Dark photon production in charged kaon decays: The $K^+\to \pi^+\gamma A'$ and $K^+\to\pi^+ \pi^0A'$ decays induced by the dipole operators (<ref>) are analyzed for the massless dark photon in Refs. <cit.>.
The bounds on allowed NP contribution to $\Delta M_{K}$ in kaon mixing, assuming no cancellations with other UV contributions, limit the branching ratios to be at most ${\cal B}(K^+ \to\pi^0\pi^+ A')\lsim 1.6\times 10^{-7}$ and ${\cal B}(K^+ \to \pi^+\gamma A')\lsim 1.1\times 10^{-8}$ for the UV model and hadronic inputs of <cit.> with $\alpha_D\equiv e^{\prime 2}/4 \pi=0.1$. Using the same UV model and $\alpha_D$ but the hadronic inputs of <cit.> one obtains, instead, the upper limits  ${\cal B}(K^+ \to\pi^0\pi^+ A')\lsim 1.3\times 10^{-7}$ and ${\cal B}(K^+ \to \pi^+\gamma A')\lsim 3.2\times 10^{-8}$.
These branching ratios are within the reach of the NA62 experiment (Sections <ref>, <ref>). Similar results are expected to apply for a massive dark photon in the limit of small mass. The $K^+\to\pi^+A'$ decay is not allowed by angular momentum conservation for a massless dark photon, while being allowed for a massive dark photon (Section <ref>), with a branching ratio expected to be of the same order as for $K^+\to \pi^+\pi^0A'$. Production of multiple dark photons in kaon decays is discussed in Section <ref>.

* Dark photon production in neutral kaon decays:
Dark photon emission in FCNC decays of $K_L$ and $K_S$ have been explored for the massless case <cit.>, in particular in the  $K_L\to\gamma A'$, $K_L\to\pi^0 \gamma A'$ and $K_L\to\pi^+\pi^- A'$ processes and their $K_S$ counterparts. The $K_L$ decays have larger branching ratios compared to $K_S$ because  of the longer lifetime of the decaying state. Using upper bounds on $\Delta M_{K}$ from kaon mixing, and assuming a vanishing electric-dipole contribution in the $s\to d A'$ transitions, the upper limits on branching ratios are expected to be <cit.>: ${\cal B}(K_{L(S)}\to \gamma A') < 6.8 (1.2) \times 10^{-5(-7)}$, ${\cal B}(K_{L(S)} \to\pi^0\gamma A') < 5.9 (1.0)\times 10^{-8(-10)}$,  ${\cal B}(K_{L(S)}\to\pi^+\pi^- A') < 5.6(9.7) \times 10^{-7 (-10)}$. These $K_L\to(\pi^0)\gamma A'$ decay rates are within the sensitivity reach of the KOTO experiment (Sections <ref>, <ref>).

* Dark photon production in hyperon decays:
FCNC dipole operators may mediate $s\to dA'$ transitions, which induce hyperon decays into dark photons such as  $\Lambda\to n A^\prime$, $\Sigma^+\to p A^\prime$, $\Xi^{(0,-)}\to \Sigma^{(0,-)} A^\prime $, and $\Omega^{-}\to \Sigma^{-}A^\prime$. These can be used efficiently to search for dark photons <cit.>.

* Dark photon production in $\pi^0$ decays:
Kaon factories serve as the principal tools for studies of rare $\pi^0$ decays. The main source of tagged $\pi^0$ mesons in $K^+$ experiments is the $K^+\to\pi^+\pi^0$ decay with a branching fraction ${\cal B}\approx 21\%$. Searches for the $\pi^0\to\gamma A'$ decay followed by either invisible (Section <ref>) or prompt visible (Section <ref>; see Ref. <cit.> for a review) dark photon decays have been performed in kaon experiments.

* Decays into  invisible states:
The rare FCNC decays $K^+\to\pi^+ \nu\bar\nu$ and $K_L\to\pi^0\nu\bar\nu$ are processes that may be mimicked by dark photon production followed by an invisible decay <cit.>. For instance, a tree-level exchange of a virtual dark photon can  mediate $K^+\to \pi^+ Q\bar{Q}$ and $K_L\to \pi^0 Q\bar{Q}$ decays, where $Q$  are
dark fermions which escape the detector, resulting in the missing energy signature. Due to the long-distance nature of the dark-photon mediated FCNC interactions, the GN bound on ${\cal B}(K_L\to\pi^0 Q\bar{Q})$ set by the constraint on ${\cal B}(K^+\to \pi^+ Q\bar{Q})$, cf. Section <ref>, can be seemingly violated by the asymmetric effect of the different kinematic cuts enforced by the NA62 and KOTO experiments, allowing for a branching ratio for $K_L \to \pi^0 Q \bar{Q}$ within the reach of the KOTO sensitivity <cit.>. The $K_L\to$ invisible decay can be searched for by the NA64 experiment up to the $10^{-6}$ level <cit.>. The $K_L\to Q\bar Q$ decay would be an example of such process, with $K-\bar K$ mixing bounds still allowing branching ratios well above the projected reach of NA64 <cit.>.

* Dark photon production via hyperon decays in supernovae:
The proto-neutron star created right after the supernova explosion contains a thermal population of $\Lambda$ hyperons. The $\Lambda\to n A^\prime$ decays can cool the proto-neutron star via a flux of escaping dark photons, modifying the observed SN1987A neutrino flux.
This sets the most stringent bound on the strength of the $s\to d A'$ dipole operator,
$\Lambda/\big|\mathbb{D}^{sd}\big|\gtrsim1.2\times10^{10}$ GeV <cit.>. The above limit can be taken as a useful benchmark for the sensitivity that one would like to achieve by the next generation of experiments. Note however that the bound does rely on the validity of the standard core-collapse interpretation of the SN1987 A and thus may potentially even not apply at all (see refs. <cit.> for recent discussions). It does not apply also in the case of a massive dark photon with $m_{A^\prime} > 1$ MeV because of its decay into electron pairs and the re-absorption of its energy.

2cKaons   2cHyperons

Decay mode     Max branching ratio      Decay mode      Max branching ratio

$K^+ \to\pi^+\pi^0 A^\prime$     $1.3\times 10^{-7}$   $\; \big(7.8\times 10^{-11}\big)$     $\Lambda\to n A^\prime$         $1.4\times 10^{-5}$  $\; \big(8.0\times 10^{-9}\big)$

$K^+ \to \pi^+ \gamma A^\prime$     $3.2\times 10^{-8}$ $\; \big(1.9\times 10^{-11}\big)$       $\Sigma^+\to p A^\prime$       $8.3\times 10^{-7}$ $\; \big(4.9\times 10^{-10}\big)$

$K_L \to \pi^0 \gamma A^\prime$     $5.9\times 10^{-8}$ $\; \big(3.4\times 10^{-11}\big)$      $\Xi^{0}\to \Sigma^{0} A^\prime $         $5.2\times 10^{-6}$ $\; \big(3.0\times 10^{-9}\big)$

$K_L \to \pi^+\pi^- A^\prime$     $5.6\times 10^{-7}$ $\; \big(3.2\times 10^{-10}\big)$       $ \Xi^{0}\to \Lambda A^\prime$         $2.4\times 10^{-6}$ $\; \big(1.4\times 10^{-9}\big)$

$K_L\to \gamma A^\prime$     $6.8\times 10^{-5}$ $\; \big(4.0\times 10^{-8}\big)$     $\Xi^{-}\to \Sigma^{-} A^\prime $         $6.1\times 10^{-6}$ $\; \big(3.6\times 10^{-9}\big)$

$K_L\to \mbox{inv.}$     $1\times 10^{-4}$ $\;\;\; \big( 6.5 \times 10^{-10}\big) $          $ \Omega^{-}\to \Sigma^{-}A^\prime$      $6.1 \times 10^{-5} $  $\; \big(3.6 \times 10^{-8} \big)$

Upper bounds on the branching ratios for different radiative kaon and hyperon decay channels using the maximal $s\to d A^\prime$ dipole transition allowed by $\Delta M_K$ from kaon mixing. If the stronger SN constraints apply, the maximal branching ratios are those listed in the parentheses. For the hadronic inputs we have used <cit.>.

Kaon and hyperon decays induced by the $d\to s A'$ dipole transition, and upper limits on their branching ratios set by $\Delta M_K$ and the SN1987A cooling constraints, are listed in Table <ref>.

§.§ Leptonic force mediators

Authors: Datta, Krnjaic, Marfatia, Redigolo, Shuve

A leptonic force involves a scalar or vector mediator that couples predominantly to leptons. Due to stringent constraints from experiments with electron initial states, the most interesting viable leptonic mediator scenarios couple predominantly to muons, taus, and/or their neutrinos. Such models lead to signatures that are distinct from canonical hidden sector scenarios for dark photons and dark scalars, and are motivated by:

* Muon $g-2$: One of the leading motivations for leptonic forces is the longstanding discrepancy between the predicted and measured $(g-2)_\mu$ values <cit.>. If this discrepancy is indeed the result of new physics, then we expect new states coupled to muons. The new physics models are of two types <cit.>: if the chirality flip in $(g-2)_\mu$ is due to the external muon leg, the NP needs to be light, while NP can be heavy if the chirality flip is due to the states running in the loop. Rare kaon decays can shed light on the light NP models. Direct searches exclude the simplest hidden sector scenarios, such as dark Higgs bosons and dark photons, from accounting for the observed muon $g-2$ assuming the mediator decays either invisibly or to pairs of SM fermions (for a recent summary, see Ref. <cit.>)[We remark, however, that scenarios where a dark photon undergoes cascade or semi-visible decays could still be viable candidates for the muon $g-2$ anomaly <cit.>.]. If, however, the new particles have an enhanced coupling to muons relative to electrons and quarks, many of the direct searches can be evaded. For a point of reference, if there exists a 100 MeV leptonic mediator, the coupling to muons needed to account for $(g-2)_\mu$ is approximately $6\times10^{-4}$ for a scalar and  $1\times 10^{-3}$ for a vector.

* $B$ anomalies:
Several measurements in charged current and neutral current semileptonic $B$ decays  differ from the standard model predictions. An  observable which has received much attention is $R_K \equiv {\cal B}(B^+ \to K^+ \mu^+ \mu^-)/{\cal B}(B^+ \to K^+ e^+ e^-)$ as this is a clean test of lepton universality with tiny hadronic corrections. There are measurements of similar ratios with a $K^*$ in the final state.
The latest LHCb measurement of $R_K$
is <cit.>

R_K(1.1 < q^2 < 6.0   GeV^2)  = 0.846^+0.044_-0.041,

differing from the SM prediction of by $\sim 3.1\sigma$. This measurement could indicate new interactions coupling preferably to  leptons of the second and third generation.
In many models proposed to address the muon $g-2$ measurement, FCNC (Flavor Changing Neutral Current) decays of the $b$ quark are generated  with the involvement of either SM particles <cit.> or new particles<cit.>. These new FCNC contributions may resolve the $B$ anomalies in the  measurement of the semileptonic $ b \to s \ell \ell$ decays.

* Dark matter: Viable GeV-scale thermal dark matter models predict the existence of new, low-mass mediators that allow for efficient DM annihilations at temperatures  $T\lesssim M_{\rm DM}$. In the simplest models, the mediators have minimal “portal” couplings, such as mixing with the Higgs boson in the case of a scalar or mixing with the photon in the case of a vector. In recent years, there have been extensive tests of the minimal portals at accelerators, colliders, and beam dump experiments, and these tests are expected to constrain dark matter models coupled through such portals <cit.>.

Leptonic force mediators serve equally well as conduits between the DM and SM, but are much less studied. As a result, for a conclusive test of thermal DM models at or below the GeV scale, it is important to test also leptonic mediators <cit.>.

Kaon factories are ideal experiments for probing low-mass leptonic mediators as the dominant charged kaon decay mode is $K^+\to\mu^+\nu$. Thus the dominant leptonic mediator production mode is through final-state radiation off the muon or muon neutrino, $K^+\to\mu^+\nu X$. In this section we argue that kaon factories can be sensitive to promising regions of parameter space that are favored by the above motivations.

§.§.§ Models

*Flavor non-universal gauge boson: Low-mass vectors can couple directly to SM leptons if lepton flavor symmetries are embedded in gauge symmetries. The model is UV complete if the gauge symmetry is anomaly free. A simple example of this is gauged differences of lepton flavor number, such as $L_\mu-L_\tau$ or $L_e-L_\mu$ <cit.>. Since vectors coupled to electron number face stringent constraints from $B$-factories and other experiments with electron initial states, we focus on $L_\mu-L_\tau$ as a representative model with $Z'$ mediator  (Ref. <cit.> showed that only small deviations from the $L_\mu-L_\tau$ model are phenomenologically allowed). Its interaction Lagrangian is

\begin{equation}
\mathcal{L} = -g'Z'_\alpha \left(\bar\mu \gamma^\alpha \mu - \bar\tau\gamma^\alpha\tau+
\bar\nu_{\mu\mathrm{L}} \gamma^\alpha \nu_{\mu\mathrm{L}} - \bar\nu_{\tau\mathrm{L}}\gamma^\alpha\nu_{\tau\mathrm{L}}\right).
\end{equation}

For all $Z'$ masses, there is an appreciable invisible branching fraction due to its coupling to neutrinos. Between the dimuon and ditau thresholds, the $Z'$ has an approximately equal branching fractions to muons and invisible products, while below the dimuon threshold the invisible branching fraction is approximately 100%. For masses $m_{Z'}\gtrsim210$ MeV, the model is constrained by direct searches at BaBar <cit.>, and indirect constraints from neutrino trident processes <cit.>. However the model is relatively unconstrained below the dimuon threshold, which is precisely the region of greatest sensitivity for kaon experiments <cit.>.

*Leptophilic scalar:
A leptophilic scalar, $\phi$, couples predominantly to leptons with suppressed or absent quark couplings <cit.>:

\begin{equation}
\mathcal{L} = -y_{ij}\,\bar\ell_i\phi \ell_j+\mathrm{h.c.}.
\end{equation}

To evade stringent constraints from charged lepton flavor violation experiments, we only consider flavor-diagonal couplings, $y_{ij}=y_i\delta_{ij}$. There are several possible textures for the Yukawa couplings:

* Mass proportional coupling, $y_i\propto m_i/v$. This interaction occurs in models with leptophilic Higgs sectors: for example, in two-Higgs doublet models where one of the doublets gives mass exclusively to the leptons <cit.>. Mixing between the leptophilic doublet and a low-mass singlet scalar, $\phi$, gives rise to this structure of interactions. In this case, the scalar $\phi$ couples predominantly to $\tau$ leptons and gives suppressed coupling to electrons, allowing the model to evade stringent constraints from accelerators and $e^+e^-$ colliders.

* Muonphilic scalar, $y_e=y_\tau=0$. This can arise in a model with a new weak-scale vectorlike lepton, $\Psi$, with the same gauge charges as the SM RH muons <cit.>.  The couplings $\phi\bar\Psi\mu_{\rm R}$ and $H^*\bar\Psi \mu_{\rm L}$ generate the $\phi\bar\mu\mu$ coupling when the vectorlike leptons are integrated out. This scenario can account for $(g-2)_\mu$ with a minimum of couplings to other SM leptons.

* Generic texture. This can arise in models with multiple vectorlike leptons. The principal requirement is that $y_e\ll y_\mu ,y_\tau$ to avoid constraints from experiments with electron initial states.

§.§.§ Branching fraction benchmarks

Because the production of leptonic mediators $X$ predominantly arises from muon FSR in $K^+\to\mu^+\nu X$ decays, the branching fraction is determined by the coupling to muons. The same coupling determines the $X$ contribution to $(g-2)_\mu$, therefore it is possible to map the models that account for $(g-2)_\mu$ to specific kaon branching fractions. The only model dependence arises in the spin of $X$, since the scalar and vector FSR rates are different.

There is more ambiguity in mapping the DM parameter space to the kaon branching fraction. Even assuming that DM thermalizes with the SM and its abundance is determined through thermal freeze-out, the specific process that sets the relic abundance depends on the mediator/DM mass hierarchy. If the mediator is lighter than DM (the “secluded scenario” <cit.>), then DM annihilates predominantly into the mediator and the muon coupling plays no role in establishing the DM abundance. Conversely, if the mediator is heavier than DM, then with specific assumptions about the magnitude of the DM-mediator coupling  we can determine the coupling to muons needed for viable DM[For the mass range probed at kaon factories, the coupling to taus is unimportant because taus are largely absent from the early-universe plasma at temperatures of $\mathcal{O}(100~\mathrm{MeV})$.]. In determining the range of parameters consistent with DM, we consider a vector mediator coupled to three possible types of DM (scalar, Majorana, and pseudo-Dirac).

Using the results of Ref. <cit.>, we derive the target branching fractions motivated by $(g-2)_\mu$ and DM, showing our results in Fig. <ref>. Obtaining a sensitivity to branching fractions $\ge10^{-9}$ would cover the most interesting parameter space for both DM and $(g-2)_\mu$. For masses smaller than 10 MeV, thermal DM models can run into constraints from $\Delta N_{\rm eff}$, the number of additional light degrees of freedom in the early universe,  while for larger mediator masses near the kinematic endpoint, Belle II and other higher-energy experiments can provide complementary sensitivity to kaon factories. Even a null result is interesting: since leptonic mediators are among the last viable models explaining $(g-2)_\mu$ at low masses, branching fraction exclusions of order $10^{-9}$ would decisively exclude these light mediator models.

Model-dependent constraints from a BaBar search for a leptonic scalar produced in association with $\tau$ leptons and decaying to $e^+e^-$ already rules out nearly the entire $(g-2)_\mu$ parameter space in the case that the scalar couples mass proportionally to leptons <cit.>, but these constraints do not apply if the mediator has a reduced coupling to taus.

Branching fraction benchmarks for $K^+\to \mu^+ \nu_\mu X$, where $X$ is a leptonic force mediator for (left) $m_X < 225$ MeV; (right) $m_X\ge 225$ MeV. The red (green) band shows the branching fraction consistent with the observed discrepancy in $(g-2)_\mu$ for a vector (scalar) $X$. The grey shaded range indicates the preferred parameter space for the DM relic abundance for a vector mediator, with a mediator that is three times heavier than DM and with a coupling of 0.03. The band indicates the parameter space spanned by models where DM is a scalar, Majorana fermion, or Dirac fermion <cit.>.

§.§.§ Signatures

The predominant production mode for the leptonic force mediator is $K^+\to\mu^+\nu X$. The experimental signatures are driven by the decay mode of $X$, and the best-motivated decay modes are the following.

* $X\to\mu^+\mu^-$ (Section <ref>): because the $X$ coupling to muons is always present for leptonic mediators, this decay mode is always present at some level. It is the dominant decay mode in models where $X$ is lighter than all other hidden sector particles, and $2m_\mu \le m_X\le m_K-m_\mu$. The dimuon decay is typically prompt for the couplings accessible at kaon factories. The projected sensitivity of NA62 in this channel has been derived in Ref. <cit.>.

* $X\to\rm inv$ (Section <ref>): this decay mode is present if $X$ couples to neutrinos, and dominates if $X$ can decay into invisible dark matter states. This decay mode is also relevant if $X$ does not couple to electrons: while $X\rightarrow\gamma\gamma$ proceeds at one loop, the lifetime can be long enough that $X$ typically decays outside of the detector. Another option is $X$ decaying to neutrinos, discussed in more detail in Section <ref>.

* $X\to \gamma\gamma$ (Section <ref>): this loop-induced decay mode dominates if $X$ does not couple to electrons and is lighter than the dimuon threshold. The partial width depends on the UV completion of the model and can lead to prompt or displaced diphoton decays. If $X$ is a scalar that couples only to muons, the proper decay length to diphotons in the unexcluded parameter space below the dimuon threshold is in the 1 cm–1 m range <cit.> if the photon coupling is generated purely from muon loops. However further UV contributions to the photon coupling could sensibly change this prediction so that a comprehensive range of lifetimes should be searched for experimentally.

* $X\to e^+e^-$ (Section <ref>): while not present in every leptonic mediator model, this decay mode dominates below the dimuon threshold in models with mass-proportional couplings of $X$. This is therefore an important channel for discovering or constraining leptonic mediators. Searches for both prompt and displaced decays are of interest, with proper decay lengths in the 0.1–10 cm range predicted in models with mass-proportional couplings of $X$ <cit.>.

§.§ Strongly interacting neutrinos

Authors: Datta, Marfatia, Zupan

Sterile neutrinos interacting with light particles have been proposed as an explanation for the anomalies in neutrino experiments such as the low energy MiniBoone excess, see e.g. <cit.>.[The goal of the MiniBooNE experiment was to address the 3.3$\sigma$ LSND anomaly in electron-like events seen in the $\bar\nu_e$ channel <cit.>. After 15 years of data taking the MiniBoone dataset itself now shows a 4.8$\sigma$ excess in the low energy part of electron spectra in both the neutrino and antineutrino channels <cit.>, possibly not inconsistent with the LSND anomaly.]
Similar interactions/models have been used to explain the XENON1T electron excess <cit.>. The light mediator would also induce neutrino self-interactions that can potentially explain the Hubble tension <cit.> or allow for sterile neutrino dark matter <cit.>. The mediator that is light enough would be produced in rare kaon decays, $K^+\to \mu^+\nu X$, via bremsstrahlung of $X$ off the final state neutrino, so that rare kaon decays offer a unique opportunity to test for such neutrino interactions (similar signatures where $X$ is radiated off the charged lepton are discussed in Section <ref>).

The phenomenology depends on whether the SM neutrinos are Majorana or Dirac, and on the properties of the mediator. The two popular choices for a light mediator are a light vector  boson, $V_\mu$, or a (pseudo-)scalar, $\phi$. If the SM neutrinos are Majorana, the interactions with the scalar mediator take the form <cit.>

L_S=1/2 g_ϕ^αβν̅_Lα^c ν_Lβ ϕ+h.c.,
with the summation over generational indices $\alpha,\beta=1,2,3$ understood. The couplings $g_\phi^{\alpha\beta}$ are not necessarily related to neutrino masses and can be sizable. This is the case, for instance, if $\phi$ is the Majoron from the inverse see-saw model where the lepton number breaking vev is small. The MiniBoonE low energy excess could be explained for $g_{\phi}^{e\alpha}={\mathcal O}(1)$ or $g_{\phi}^{\mu\alpha}= {\mathcal O}(1)$ and $m_\phi= {\mathcal O}(50~{\rm MeV})$. However such large couplings are already excluded by the upper bounds of the ${\cal B}(K^+\to\ell^+\nu X_{\rm inv})$ decay rates (Section <ref>) <cit.>. Improvements on ${\cal B}(K^+\to\ell^+\nu X_{\rm inv})$, on the other hand, could probe the self-interacting neutrino explanation of the Hubble tension <cit.>. The possibility of $\phi$ coupling to only electron neutrino is excluded (Fig. <ref>, left), while an improvement of a factor of a few on ${\cal B}(K^+\to \mu^+\nu X_{\rm inv})$ would probe fully the preferred region in the case of $\phi$ coupling exclusively to the muon neutrinos (Fig. <ref>, right). The possibility of coupling to tau neutrinos remains wide open, and is not probed by rare kaon decays. The SM neutrinos in (<ref>) can also be replaced by light, MeV or keV, sterile neutrinos, $N_i$, if these couple to light scalars (see also discussion on HNLs in Section <ref>). The couplings of $N_i$ to $\phi$ are in general only very weakly constrained.

Bounds (purple shaded regions for constrains from $\tau$ and rare meson decays, red for BBN constraints) on light neutrino-coupled mediators with flavor- specific couplings to either $\nu_e$ (left) or $\nu_\mu$ (right). The bands labeled MI$\nu$ and SI$\nu$ are the  regions in the $g_\phi-m_\phi$ plane preferred by the Hubble tension.  If neutrinos are Dirac, their right-handed components equilibrate before BBN above the dashed black line. Plots taken from <cit.>.

Similar considerations apply to the Dirac neutrinos. The coupling between $\phi$ and the SM neutrinos can be generated from the $\phi$ interactions with sterile neutrinos, ${\cal L}=- g_D \phi \bar {\nu}_D  \nu_D + \ldots,$ if the sterile neutrinos mixes with the lighter SM neutrinos. The phenomenology can be more involved, if $\phi$ mixes with the Higgs, and resembles the phenomenology of the Higgs portal scalar in Section <ref>, but with the large couplings of $\phi$ to neutrinos allowed. In particular, the ${\cal B}(K\to\ell \nu_\ell \phi)$ can be followed by decays to visible SM particles as well as by decays to neutrinos, with the relative branching ratios controlled by the size of $g_D$ and the mixing angles. A UV complete realization of such a model was discussed in Refs. <cit.>, where a light scalar singlet was added to a two-Higgs-doublet model. In addition decays through either on-shell or off-shell sterile neutrino, such as $K^+\to\ell \bar\nu_D \to\ell \bar\nu\phi$ with $\phi$ decaying to visible/invisible final states are possible.

The phenomenology of the light vector mediator also follows the above considerations in broad strokes. The light $Z'$ can couple to either left-handed or right-handed neutrino components. For instance if the sterile neutrino carries a hidden $U(1)$ gauge charge, with $Z'_\mu$ the corresponding gauge boson, the interaction Lagrangian is given by $\mathcal{L}_{V} \supset  g' Z'_\mu \bar {\nu}_D \gamma^\mu  \nu_D= g' Z'_\mu \left(\bar {\nu }_{DL} \gamma^\mu  \nu_{DL}+\bar{\nu} _{DR} \gamma^\mu  \nu_{DR} \right)$. The sterile-active neutrino mixing then leads to couplings of $Z'$ to the SM neutrinos, with similar phenomenology to the one for the scalar mediator discussed above (see also Section <ref> for more complex decays involving HNLs). An explicit model of a light $Z'$ coupling predominantly to neutrinos can be found in <cit.>.

§.§ Grossman-Nir violating models

Authors: Grossman, Hostert, Kitahara, Pospelov,  Tobioka

There is a robust theoretical relation between
${\cal B}(K_L \to \pi^0 \nu \bar\nu)$ and ${\cal B}(K^+\to \pi^+\nu\bar\nu)$, valid in any model that satisfies some mild assumptions: isospin symmetry between $K_L$ and $K^+$ decays, the dominance of the $\Delta I=1/2$ amplitude, and the smallness of CP violation in the decay <cit.>.
From this one obtains

≃     | p A - q A̅|^2/|√(2) A|^2
= |p|^2/2 | 1 - λ|^2
= |p|^2/2
( 1+|λ|^2 - 2 Reλ)

≃    sin^2 [Arg(λ)/2 ]
≤1 ,

where $K_L = p K^0 - q \bar{K}{}^0$,
$A = \langle \pi^0 \nu \bar \nu |\mathcal{H}|K^0 \rangle$, $\bar A = \langle \pi^0 \nu \bar \nu |\mathcal{H}|\bar{K}{}^0 \rangle$, and $\lambda = (q/p)(\bar{A}/A)$.
In the approximation in the first line
the isospin breaking effects are discarded and  $\Delta I=1/2$ dominance is assumed.
The approximation in the second line assumes that
CP violation is negligible and thus
$|p| =|q|= 1/\sqrt 2$ and $|\lambda|=1$.
The last inequality in (<ref>) translates into a
Grossman-Nir (GN) bound,

\begin{equation}
{\cal B}(K_L\to\pi^0\nu\bar\nu) \leq
4.3\cdot{\cal B}(K^+\to\pi^+\nu \bar\nu),
\label{eq:GNbound}
\end{equation}

where summation over different neutrino species
is understood <cit.>.
The numerical factor accounts
for the difference in the $K_L$ and $K^+$
total decay widths,
the isospin breaking effects, and QED radiative corrections <cit.>.

Since the SM satisfies all of the above assumptions,
the SM branching ratios fulfill the GN bound

The latest NA62 analysis of Run 1 data set results in the definitive measurement of this branching ratio: ${\cal B}(K^+ \to\pi^+\nu \bar\nu) = (10.6^{+4.0}_{
-3.4}|_{\rm stat} \pm 0.9_{\rm syst}) \times 10^{-11}$ <cit.>. The most stringent limit from the KOTO experiment is ${\cal B}(K_L \to \pi^0\nu\bar\nu)
<3.0\times 10^{-9}$ at 90% CL <cit.>, while the standard GN bound (<ref>) indicates ${\cal B}(K_L\to\pi^0 \nu \bar\nu)<7.7\times 10^{-10}$.

Similarly to the $K_L \to \pi^0 \nu \bar\nu$ case,
a generalized version of the GN bound holds for the production of light NP state(s) $X$,

\begin{equation}
{\cal B}(K_L\to\pi^0 X ) \leq
4.3\cdot{\cal B}(K^+ \to \pi^+ X),
\label{eq:generalGNbound}
\end{equation}

if the assumptions are satisfied for the NP decay (for $K\to \pi\pi X$ see <cit.>). The inequality can be saturated
in case $X$ is a CP-even scalar <cit.>.
If $X$ is a single neutral particle, the current NA62 limit is ${\cal B}(K^+ \to \pi^+ X_{\rm inv})<\mathcal{O}(10^{-11})$,
where the exact numerical value has a strong dependence on the particle mass $m_X$, as well as on its lifetime $\tau_X$ when it decays into visible SM particles <cit.>.
The generalized GN bound (<ref>) then leads to ${\cal B}(K_L \to \pi^0 X_{\rm inv})<\mathcal{O}(10^{-10})$.

Because of the inequality (<ref>) rare $K^+$ decays set stronger constraints on the parameter space of most NP models, compared to the rare $K_L$ decays. Nevertheless,  there exist also NP models that violate the assumptions of the (generalized) GN bound, for instance, $K_L$ can mix directly with dark sector states, while the charged $K^+$ cannot. In such cases the $K_L$ decays place stronger constraints. The violation of the generalized GN bound may be a fundamental feature of the NP sector, or it can just be due to the experimental loopholes.

The interest in NP models violating the generalized GN bound increased after the KOTO
experiment observed three events in their signal region <cit.>, where
the initial preliminary background estimate was $0.05$ events, later revised to $1.22\pm 0.26$ after the unblinding. Since KOTO is not able to measure the $\gamma\gamma$ invariant mass and is, therefore, unable to determine if they come from a $\pi^0$ or not, the GN violation can also be just apparent, from enhanced $K_L$ decays to $\gamma\gamma+{\rm inv}$ final states, without intermediate $\pi^0$. Besides the interest in the KOTO excess, which might end up being explained by an underestimated background, we will try to emphasize the different physics mechanisms that were put forward. The latter typically points to new signals, or different analysis to be performed at kaon factories. Also, as we will show below, reasonable scenarios require the existence of new light particles, which is in accord with the concept of this review.

There are five known ways in which the generalized GN violation is achieved in the NP scenarios:

$(i)$ CP violation in the decay;

$(ii)$ $\Delta I=3/2$ dominance;

$(iii)$ exploiting charge conservation in the decay;

$(iv)$ exploiting mass difference;

$(v)$ effective violation due to experimental loopholes.

The first four items in the above list
violate the GN bound fundamentally, in the sense that the corresponding $K^+$ decay is either absent or does not obey the relation in Eq. (<ref>). The effective violation can avoid Eq. (<ref>) by exploiting loopholes in the experimental searches.
We consider both the cases of heavy ($\Lambda_{\rm NP} \gg m_W$) and light
($\Lambda_{\rm NP} < m_K$) NP.
If the GN violation is fundamental, the heavy new physics is implausible, and instead light new states are favored. For the effective violation, we only need one light dark state, which would also be testable with improved $K^+$ analyses.

§.§.§ Fundamental Violation

§.§.§ $(i)$ CP violation in the decay

CP violation in the decay, in particular, that implies $|\lambda |\gg 1$ in Eq. (<ref>), can significantly violate the GN bound. In that case, we also have large CP violation in the charged mode, that is, the $K^+$ and $K^-$ decay rates are different.  Since NA62 measures only one type of charge ($K^+$), this idea is not experimentally ruled out. It can be tested in principle by measuring the charged modes also with $K^-$.

Theoretically, this idea requires both strong and weak phases to be present, and the strong phases from QCD are generically expected to be
suppressed  as the final states involve neutrinos (or other SM
singlets). Thus, this idea cannot work with heavy NP. In principle,
light NP can generate a strong phase via intermediate particles that are on-shell. We are unaware of any realistic model that can do it, but that might be due to the fact that no investigation in this
direction was made, and we cannot conclude that it is impossible.

§.§.§ $(ii)$ $\Delta I=3/2$ dominance

The GN bound assumes $\Delta I=1/2$ dominance. Without new light states, making the $\Delta I=3/2$ transition dominant is the only known way to violate the GN bound <cit.>. The SMEFT operators that can violate the GN bound are required to contain at least all the relevant fermions (four quarks and two neutrinos), and are thus at least of dimension nine. While the numerical coefficient of the GN bound in Eq. (<ref>) can be brought from $4.3\to 20$, the suppression scale of the SMEFT operators needs to be as low as few GeV <cit.> and hence naively incompatible with the SMEFT framework. Therefore, it seems implausible the GN bound can be appreciably violated only using heavy new physics.

The $\Delta I=3/2$ transition may become dominant more easily
if the SM neutrinos in the final state are replaced with a new light singlet state, $S$ <cit.>. The dimension of the $\Delta I=3/2$ operator is then reduced by two, and even for the suppression scale as high as a TeV, the GN bound is still significantly relaxed,

\begin{align}
{\cal B}(K_L\to \pi^0 S ) \leq
51 \cdot {\cal B}(K^+ \to \pi^+ S).
\label{GNboundDim7}
\end{align}

The singlet scalar is expected to be invisible in kaon experiments.
The current NA62 limit on ${\cal B}(K^+\to\pi^+X_{\rm inv})$ is at the ${\mathcal O}(10^{-11})$ level <cit.> (Fig. <ref> and Section <ref>), which implies that ${\cal B}(K_L\to \pi^0 X_{\rm inv}) \lesssim 5 \times 10^{-10}$ for these types of NP models, well within reach of future KOTO and KLEVER sensitivities.

$K_L$ decay     Model

$K_L \to \gamma\gamma +\textrm{inv}$

search strategy

$K^+\to \pi^+ +\textrm{inv}$

Other predictions

3*$ \to \pi^0 X$    (d)  <cit.>

3*current  $\pi^0 X_{\rm inv}$     SM-like
$K^0$–$\bar K^0$ mixing

$X_{1,2}$  at beam-dump

$\Delta I=\frac{3}{2}$ <cit.>

${\cal B}(K_L \to \pi^0 X ) %\leq

$\leq51 {\cal B}(K^+ \to \pi^+ X)$

5*$ \to \pi^0 X_1  X_1 $     3*(c)  <cit.>
5*current $\pi^0 \nu\bar\nu$

${\cal B} (K_L \to \pi^0 X_1 X_1)$

$\sim 4 (4 \pi)^2 {\cal B} (K^+ \to \pi^+ X_i X_j)$

if $m_{X_i}+m_{X_j} < m_{K} - m_{\pi}$

$B\to K (\pi^0) XX$ (in MFV)

$X_2$ at beam dumps

2-24-5     (e)  <cit.>            SM-like

$K^0$–$\bar K^0$ mixing

$X_{2}$  at beam-dump

3*$\to \gamma\gamma X_1$     (a) <cit.>
3* optimize $\gamma\gamma$

${\cal B}(K_L \to   \gamma\gamma X_1)$

$\sim 4 (4 \pi)^2 {\cal B}(K^+ \to \pi^+ \gamma\gamma X_1 ) $

$\sim 4 (4 \pi)^2 {\cal B}(K^+ \to \pi^+ X_1 X_1 ) $

if $m_{X_i}+m_{X_j} < m_{K} - m_{\pi}$

${\cal B}(K_L\to \pi^0 XX)

$\sim 10^{-2}{\cal B}( K_L \to \gamma\gamma X )$

(a) <cit.>

$B\to K (X_2 \to \gamma\gamma) X_1$ (in MFV)

$X_2$ at beam dumps

$ \to \gamma\gamma X_1 X_1  $     (b) <cit.>        optimize $\gamma\gamma$

${\cal B} (K_L \to \gamma\gamma X_1 X_1)$

$\sim 4 (4 \pi)^2 {\cal B} (K^+ \to \pi^+ \gamma\gamma X_1 X_1)$

if $m_{X_1}+m_{X_2} < m_{K} - m_{\pi}$

$B\to K(\gamma \gamma) XX$ (in MFV)

$X_2$ at beam dumps

The NP scenarios that lead to violations of the generalized GN bound at the fundamental level. In all cases $X_1(X)$  is an invisible particle.
Some of the $K^+$ decays are kinematically allowed only for $2m_X<m_K-m_\pi$.

§.§.§  $(iii)$ Exploiting Charge Conservation in the Decay

Without relying on $\Delta I=3/2$ dominance, one can violate the generalized GN bound if there are new light dark sector particles in the sub-GeV mass range. The basic idea is to exploit the charge difference between $K_L$ and $K^+$. The dark sector particles $X_i$ can only couple to electromagnetic neutral operators of the SM, such as $\bar s d$, meaning that operators like  $X_i X_j \bar s d$ would lead to neutral kaon decays to the dark sector without necessarily predicting analogous decays in charged kaons. These types of scenarios are discussed in Refs. <cit.>.

Typically, these models require at least two states in order to generate the GN-violating $K_L\to \pi^0/\gamma\gamma+\textrm{inv}$ decays. After the direct decay of $K_L$ to the dark sector via $K_L\to X_2 X_{1,2}$, the dark states ought to transition back to pion-like final states. This can easily be achieved through subsequent decays such as $X_2\to \gamma\gamma$, $X_2\to X_1 \gamma$, or $X_2\to  X_1 \pi^0$, where $X_1$ is stable and invisible. These scenarios can be categorized as follows,

\begin{align}\label{eq:pair-production-categories}
& \text{(a) $\pi^0$ impostor \cite{Gori:2020xvq, Hostert:2020gou}}
&&K_L\to X_2 X_1\ \&\ X_2\to \gamma\gamma,
\\
&\text{(b) dipole portal \cite{Hostert:2020gou}}
&&K_L\to X_2 X_2 \ \&\ X_2\to X_1\gamma,
\\
&\text{(c) $\pi^0$ production \cite{Hostert:2020gou}}
&&K_L\to X_2 X_1 \ \&\ X_2\to X_1\pi^0.
\end{align}

In Ref. <cit.>, UV completions for all of the above scenarios are provided using either a Higgs or a $Z'$ portal, all the while assuming Minimal-Flavor-Violation (MFV). If this assumption holds, $b\to s,d$ transitions will also be sensitive to potential signals in $K_L$ decay. Of particular relevance is the process $B\to K X_2 X_{1,2}$, which can be searched for at Belle-II, with $B\to K + \textrm{inv}$ constituting the most promising channel when $X_2$ particles decay outside the detector. Another model example is provided by Ref. <cit.>, which introduces a complex scalar with a half strangeness whose real and complex parts correspond to $X_1$ and $X_2$, respectively. This generates the diphoton signal in scenario (a) of Eq. (<ref>) via the $X_2F\tilde F$ coupling, which is responsible for $X_2\to \gamma\gamma$ decays.

Another way to exploit the charge difference is through the mixing of $X_2$ with light neutral mesons. If there is only a linear coupling $X_2  \bar s d$,  both $K^+ \to \pi^+ X_2$ and $K_L \to \pi^0 X_2$ are allowed and the generalized GN bound holds. However, if $X_2$ is heaver than $m_K-m_\pi$  and couples to another light dark sector particle $X_1$, such as $X_1 X_2^2$ ($X_1^2 X_2^2$), the authors of <cit.> pointed out $K_L \to \pi^0 X_1(\pi^0 X_1X_1)$ has a contribution from mixing:

\begin{align}
&\text{(d) Model 1  \cite{Ziegler:2020ize}}
&&K_L\to X_2^* \to  X_1 \pi^0,
\\
&\text{(e) Models 2\&3 (scalar or fermion $X_1$)  \cite{Ziegler:2020ize}}
&&K_L\to X_2^* \to  X_1X_1\pi^0,
\end{align}

while $K^+\to\pi^+ X_1(\pi^+ X_1X_1)$ does not have a similar mixing contribution. The prediction of this scenario is  $X_2\bar d d$ operator for $X_2$–$\pi^0$ mixing, and $K^+ \to \pi^+ X_1(\pi^+ X_1X_1)$ is not kinematically forbidden.
The violation of the GN bound is achieved if a diagram  with mixing between the neutral mesons ($K^0, \pi^0$) and $X_2$ dominates in the neutral kaon decay. The $X_2$ cannot be too heavy because the mixing contribution decouples fast as $m_{X_2}^{-8}$. As long as the intermediate scalar is light enough, the correction to $K_L\to \pi^0 +\textrm{inv}$ can be orders of magnitude above the SM while the correction to $K^+\to \pi^+ +\textrm{inv}$ is order one of the SM branching ratio. Since the mixing between $X_2$ and $K^0$ is essential, $K^0$–$\bar K^0$ oscillation sets some stringent constraint, but viable parameter space still exists. Another essential coupling, $X_2 \bar d d$, predicts signals from $X_2$ production or double $X_1$ production at beam-dump experiments.

So far we only discussed models where the NP couples to flavor-non-diagonal operators. This need not be the case, as pointed out in <cit.>, as $X$ particles could be produced from flavor-diagonal couplings to the first generation of quarks. This can happen via long-distance ($\Delta S = 1$) induced $K_L$ transitions to virtual unflavored mesons ($\pi^0, \eta$, and $\eta^\prime$), which in turn couple to the dark sector in a flavor diagonal fashion ($\pi^0 X_1 X_2$, for instance). The decays $X_2 \to X_1 \pi^0$ are then unavoidable, if kinematically allowed, and will naturally lead to scenario (c) in Eq. (<ref>).

The scenarios discussed above are categorized in the experimental perspectives in Table <ref>. At KOTO/KLEVER, the standard $\gamma\gamma+\text{inv}$ analysis assumes that the two photons are from a $\pi^0$, so the origin of two photons is one criterion. Another criterion would be the number of the invisible particles. The scenarios (c,d,e) and $\Delta =3/2$ dominance scenario will be well tested by the standard $K_L\to \gamma\gamma+\text{inv}$ search, however the standard searches are not optimized for scenarios (a,b). While in scenario (a) it would be possible to adapt the searches to different values of $m_{\gamma \gamma}$, in scenario (b) the two gammas are uncorrelated and will be selected by the standard analyses much less often. Nevertheless, the efficiencies for all cases will be nonzero, depending on the masses, and are shown in the left panel of Fig. 12 of Ref. <cit.> and in Fig. 8 of Ref. <cit.>.

We emphasize that in all of the scenarios discussed above, $K^+\to \pi^+XX$ would be allowed only if $2 m_{X_{1}}<m_K-m_\pi$, and the branching fraction would be smaller when compared to $K_L \to X_1\gamma\gamma$ due to the reduced phase space and the three-body nature of the decay.

The models presented here, with simple modifications, can predict  other neutral kaon decay patterns, while still having new physics effects in $K^+$ decays either suppressed or forbidden. For example, if   $X_2$ in scenario (a) can decay through a dilepton channel, the model then predicts the $K^0 \to \ell^+\ell^-+{\rm inv}$ signature.  Instead, if $X_1$ is unstable and decays to leptons, then the $K^0 \to \ell^+\ell^-\gamma\gamma$ signature is predicted. The experimental searches of these channels are discussed in Section <ref>.
If both $X_1$ and $X_2$ decay to leptons this gives the  $K^0\to \ell_\alpha^+\ell_\alpha^-\ell^+_\beta\ell_\beta^-$ signature, with the  experimental status discussed in Section <ref>. See also Section <ref> for more concrete realizations of such models.

§.§ $(iv)$ Exploiting isospin dependence of light meson masses

Within the SM, radiative corrections produce a mass difference between charged and neutral mesons. As a result, the difference between $K$ and $\pi$ masses is different between the charged and neutral modes. If a new particle $Q$ satisfies the following mass relation:

\begin{equation}
m_{K^+} - m_{\pi^+} =354\,{\rm MeV} < 2 m_{Q}< m_{K_L} - m_{\pi^0} = 363\,{\rm MeV},
\label{eq:kinematicsrange}
\end{equation}

the three-body kinematics maximally violates the isospin symmetry and the generalized GN bound <cit.>.
While $K_L \to \pi^0 Q \bar{Q}$ remains open,
$K^+ \to \pi^+  Q \bar{Q}$ is kinematically forbidden.
However, the predicted $\pi^0$ momentum is  $p_{T}^{\pi^0} < 60\,$MeV, and in most the events the pions are expected to decay down the beam hole and remain unobservable in the current design of the KOTO setup. In addition, for the decays that do get reconstructed, the measured $p_{T}^{\pi^0}$ will lie much below current signal selection cuts, where backgrounds, such as $K_L \to \gamma \gamma$, are expected to be large.

§.§.§ Effective violation due to experimental loopholes

We move on to discuss effective violation of the GN bound due to the specific experimental setups of the NA62 and KOTO experiments. We have identified several possibilities, discussed below.

§.§.§ $(v$-$i)$ $\pi^0$ blind spot

An additional neutral state $X$ is fundamentally constrained by the generalized GN bound in Eq. (<ref>).
However this bound can be effectively violated by resorting to differences in experimental efficiency and background levels.
When one considers the case of $m_{X}\approx m_{\pi^0}$, with $X$ decaying invisibly, constraints
from $K^+\to\pi^+X_{\rm inv}$ searches at the E949 and NA62 experiments are significantly loosened due to the $K^+\to\pi^+\pi^0$ background <cit.>, leading to the $\pi^0$ blind spot. Even under these circumstances, an upper bound on ${\cal B}(\pi^0 \to $ inv) in the $K^+\to\pi^+\pi^0$ channel,

\begin{equation}
{\cal B}
(\pi^0 \to {\rm inv}) < 4.4 \times 10^{-9} \quad (90\%\,{\rm CL}),
\end{equation}

has been adapted to constrain
${\cal B}(K^+\to\pi^+X)$ for the case with $m_{X} = m_{\pi^0}$
as done recently by the NA62 collaboration <cit.>:

\begin{equation}
{\cal B}(K^+ \to \pi^+ X) \lesssim 10^{-9} \quad (\tau_X = \infty).
\end{equation}

Compared to the current NA62 limit of ${\cal B}(K^+ \to \pi^+ X_{\rm inv} ) < \mathcal{O}( 10^{-11})$ for
$\tau_X = \infty$ at 90% CL <cit.>,