1em captionlabel caption titlehead title section headinclude=false,
footinclude=false, twoside=false, parskip=false, draft=false, abstract=true,
numbers=noenddot, DIV=12 DO-TH 23/17

# The Anatomy of $K^{+}\to\pi^{+}\nu\bar{\nu}$ Distributions

Martin Gorbahn<EMAIL_ADDRESS>Department of Mathematical Sciences,
University of Liverpool, Liverpool L69 3BX, UK Ulserik Moldanazarova
<EMAIL_ADDRESS>Department of Mathematical Sciences, University of
Liverpool, Liverpool L69 3BX, UK Faculty of Physics and Technology, Karaganda
Buketov University, 100028 Karaganda, Kazakhstan Kai Henryk Sieja
<EMAIL_ADDRESS>Fakultät für Physik, TU Dortmund, D-44221 Dortmund,
Germany Emmanuel Stamou<EMAIL_ADDRESS>Fakultät für Physik, TU
Dortmund, D-44221 Dortmund, Germany Mustafa Tabet<EMAIL_ADDRESS>Fakultät für Physik, TU Dortmund, D-44221 Dortmund, Germany

(December 10, 2023)

###### Abstract

The excellent experimental prospects to measure the invisible mass spectrum of
the $K^{+}\to\pi^{+}\nu\bar{\nu}$ decay opens a new path to test generalised
quark–neutrino interactions with flavour changing $s\to d$ transitions and as
such to novel probes of Physics beyond the Standard Model. Such signals can be
a consequence of new lepton-number violating or lepton-number conserving
interactions, with their interpretations depending on the Majorana versus
Dirac nature of the neutrinos. Furthermore, the possible existence of new
massive sterile neutrinos can be tested via their distinctive imprints in the
invariant mass spectrum. Within the model-independent framework of the weak
effective theory at dimension-six, we study the New Physics effects of
Majorana and Dirac neutrinos on the differential distribution of
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ allowing for lepton-number violating
interactions and potential new sterile neutrinos. We determine the current and
expected future sensitivity on the corresponding Wilson coefficients using the
distribution measured by the NA62 collaboration and accounting for expected
improvements based on the HIKE experiment. We present single-operator fits and
also determine correlations among different type of operators. Even though we
focus on $s\to d\nu\nu$ transitions, the operator bases for Majorana and Dirac
and the classification of lepton-number-violating/conserving interactions is
applicable also for the study of $b\to s/d\nu\nu$ and $c\to u\nu\nu$
transitions relevant in current phenomenology.

###### Contents

  1. 1 Introduction
  2. 2 Weak effective theories for $K\to\pi\nu\bar{\nu}$
    1. 2.1 Majorana-$\nu$ EFT with lepton-number violation
    2. 2.2 Dirac-$\nu$ EFT with lepton-number conservation
  3. 3 Rates and distributions for $K^{+}\to\pi^{+}\nu\bar{\nu}$ at NA62
    1. 3.1 Decay rates
    2. 3.2 Statistical treatment for EFT fit
  4. 4 New Physics sensitivities
    1. 4.1 Single-operator fits
    2. 4.2 Multi-operator correlations
  5. 5 Discussion and conclusions
  6. A Decay rates

##  1 Introduction

The branching ratio of the rare $K^{+}\to\pi^{+}\nu\bar{\nu}$ decay is
currently being measured at $\mathcal{O}(35\%)$ accuracy at the NA62
experiment [1, 2], while improved upper limits for the corresponding neutral
decay mode, $K_{L}\to\pi^{0}\nu\bar{\nu}$, are provided from the KOTO
experiment [3, 4]. Each event of the charged decay mode has a so-called
missing neutrino mass squared, which leads to a decay spectrum that can be
measured with increased experimental statistics. The existing measurements
already play an important role in current phenomenology. They are sensitive to
modifications of the Standard Model (SM) that are generated via effective
operators at $\mathcal{O}(100\,\text{TeV})$ and also constrain feebly
interacting particles that can provide new final states contributing to the
invisible mass spectrum [5, 6].

This exceptional sensitivity to New Physics (NP) is due to the high
suppression of these decay rates in the SM and the precision reached in their
theory predictions. Both rare $K\to\pi\nu\bar{\nu}$ decay modes are flavour-
changing neutral-current (FCNC) processes, and as such loop-induced in the SM
and additionally suppressed by products of off-diagonal
Cabibbo–Kobayashi–Maskawa Matrix (CKM) elements. The fact that the neutrinos
are only interacting with the weak sector of the SM results in a hard
Glashow–Iliopoulos–Maiani (GIM) mechanism [7], which highly suppresses light-
quark contributions. This generates a severe CKM suppression in the SM and
facilitates a precise theory prediction, where, to a very good approximation,
only a single local charged-current operator contributes below the charm
scale. The evaluation of the operator matrix element and its Wilson
coefficient through an ongoing theoretical effort has lead to a theory
uncertainty of only $3\%$ and $1\%$ for the charged and neutral decay modes,
respectively [8, 9, 10, 11, 12, 13]. This uncertainty includes the
contribution of higher-dimensional operators and light-quark contributions,
which have been estimated in chiral perturbation theory [14] and can be
evaluated on the lattice [15, 16]. The overall uncertainty in the current SM
prediction for the charged decay mode [17] is still dominated by the CKM input
parameters, yet current input values [18] result in an overall theory
uncertainty of $\simeq 6\%$. This has to be compared with the expected
uncertainty of $15\%$ and $5\%$ for the branching-ratio measurement at the end
of NA62 and after the first phase of the proposed HIKE experiment [19],
respectively. These improved measurements, together with the ones for the
$K_{L}$ mode [20], will result in even tighter probes of potential NP
scenarios with flavour violation in the $s-d$ sector.

The present work is motivated by the fact that a more detailed measurement of
the missing-mass spectrum will become available at the end of the NA62
experiment and during the run of HIKE. In anticipation of this data, we will
consider scenarios with heavy NP, including the possibility of additional
massive sterile neutrinos, that can impact the missing-mass spectrum of
$K^{+}\to\pi^{+}\nu\bar{\nu}$. We discuss their effect in a mostly model-
independent manner by considering all the relevant weak effective theory
operators up to dimension six for both Majorana and Dirac neutrinos. Here,
non-zero contributions of (pseudo-)scalar and tensor operators can
distinctively modify the invariant mass spectrum from the SM-like spectrum
generated by vector and axial-vector operators. These generalised neutrino
interactions have recently received considerable attention in the context of
neutrino oscillations [21, 22, 23, 24] and correlations to charged current
decays have for example been studied in Ref. [25] in the context of SMEFT. In
this work, we will study their impact on the FCNC golden Kaon decay. The
impact of sterile neutrinos on (axial-)vector type operators has been studied
in Ref. [26]. Working with the two, physically distinct scenarios of Dirac and
Majorana neutrinos allows us to incorporate mass effects of sterile neutrinos
for all operators, which would not be possible in a basis where neutrinos are
represented as Weyl fields [27]. This distinction has the benefit of providing
a transparent interpretation of a possible NP signal in terms of lepton-number
conserving or lepton-number violating NP interactions. Dirac and Majorana weak
effective theories originate from different UV theories. E.g., in the Standard
Model Effective Field Theory (SMEFT) only three Majorana neutrinos are present
in the weak effective theory. Constraints on scalar-type SMEFT operators on
the $K^{+}\to\pi^{+}\nu\bar{\nu}$ distributions were partially studied in Ref.
[28], but an analysis that incorporates the experimental cuts and backgrounds
is hitherto missing. To incorporate the full experimental information we
derive the differential decay rate as a function of the pion momentum and the
invariant neutrino mass in the lab frame of NA62 for Majorana and Dirac
neutrinos including mass effects and allowing for flavour-violation in the
neutrino sector. We then perform a full statistical analysis of the
sensitivity to probe the different NP Wilson coefficients accounting for the
single-event sensitivity and available background information at the current
and future-planned experimental setup.

The paper is organised as follows. In Section 2 we classify the weak effective
theory relevant for $q\to q^{\prime}\nu\nu$ FCNCs and count the independent
Wilson coefficients that contribute to flavour-changing quark into neutrino
decays. In the following Section 3 we provide the relevant decay rates for
zero neutrino masses, and explain our statistical treatment including the
treatment of theoretical uncertainties. In the next Section 4 we provide the
current and expected future sensitivities for the cases in which a single
operator contributes (single-operator fits) as well as the correlated
constraints on pairs of different operators. After the conclusion we also
provide in Appendix A further formulae for the double differential
distributions $\mathrm{d}^{2}\mathrm{Br}/\mathrm{d}q^{2}\mathrm{d}k^{2}$.

##  2 Weak effective theories for $K\to\pi\nu\bar{\nu}$

The focus of this work is to study the possible impact of heavy New Physics
(NP) scenarios in the distributions of $K^{+}\to\pi^{+}\nu\bar{\nu}$. The goal
is to also include the possibility that the existence of additional heavy
neutrinos could provide further final states with missing energy. To this end,
it is useful and necessary to distinguish the conceptually different cases of
having Majorana-type versus Dirac-type neutrino masses for the three light
neutrinos or the possible extra sterile neutrinos. Both Majorana and Dirac
operator basis have a direct correspondence in terms of the basis written
using chiral (or equivalently two-component Weyl) fermion fields [27].
However, unless the dimension-four mass terms are specified, the chiral basis
with more than three neutrino fields leaves the question of neutrino masses
unspecified: The additional right-handed neutrinos can combine with the left-
handed ones to give Dirac masses without the need for lepton-number-violating
mass terms. Alternatively, the observed neutrinos are Majorana particles in
which case their masses originate from Majorana mass terms and the extra,
sterile neutrinos have their own independent mass.

We will thus separately discuss two cases:

  1. (i)

Three or more Majorana neutrinos with dimension-six interactions with $s-d$
FCNC interactions that can include lepton-number violation.

  2. (ii)

Dirac neutrinos with dimension-six interactions with $s-d$ FCNC interactions
that conserve lepton-number.

Each of these cases will have a distinct low-energy effective theory. We thus
decompose the effective Lagrangian relevant for quark-FCNC transitions with
neutrinos, i.e., $d_{i}\to d_{j}\nu\nu$ or $u_{i}\to u_{j}\nu\nu$, at the
dimension-six level as

$\mathscr{L}^{\nu\text{LEFT}}_{q\to
q^{\prime}\nu\nu}\Big{|}_{\nu\text{-type}}=\mathscr{L}^{(4)}_{\text{QCD}}+\mathscr{L}^{(4)}_{\nu}\Big{|}_{\nu\text{-type}}+\mathscr{L}^{(6),\text{SM}}_{q\to
q^{\prime}\nu\nu}+\mathscr{L}^{(6)}_{q\to
q^{\prime}\nu\nu}\Big{|}_{\nu\text{-type}}\,,$ (1)

where the superscript ‘‘$(n)$’’ indicates the mass-dimension of the operators
and $\nu\text{-type}$ labels the case of Majorana or Dirac neutrinos. Apart
from the usual dimension-four QCD Lagrangian involving quarks and gluons, the
kinetic and mass term of the neutrinos reads

$\displaystyle\mathscr{L}^{(4)}_{\nu}\Big{|}_{\text{Majorana}}=\frac{1}{2}\sum_{a}\bar{\nu}_{Ma}(i\not{\partial}-M_{a})\nu_{Ma}\,,$
$\displaystyle\mathscr{L}^{(4)}_{\nu}\Big{|}_{\text{Dirac}}=\sum_{a}\bar{\nu}_{Da}(i\not{\partial}-M_{a})\nu_{Da}\,,$
(2)

where $a\geq 3$ runs over the different neutrino flavours, and $\nu_{Ma}$ and
$\nu_{Da}$ denote four-component Majorana and Dirac fields, respectively. Eq.
(2) is already written in terms of mass-eigenstates for the neutrinos fields.
Throughout this work we work in the mass-eigenstate basis for the neutrinos.
This choice of basis implicitly also fixes the basis for the definition of the
dimension-six Wilson coefficients.

In the present work we focus on $K\to\pi\nu\nu$. In this case the SM induces a
single dimension-six operator below the charm mass [29, 17]. For the $s\to d$
transition the SM contribution can be written both in terms of Majorana or
Dirac fields

$\ \begin{split}\mathscr{L}^{(6),\text{SM}}_{s\to
d\nu\nu}&=\sum\limits_{\begin{subarray}{c}a,b=\left\\{1,2,3\right\\}\end{subarray}}\left\\{\begin{array}[]{c}C^{A,\,L,\,\text{SM}}_{ab21}~{}O^{A,\,L}_{ab21}\\\\[5.0pt]
C^{V,\,LL,\,\text{SM}}_{ab21}~{}O^{V,\,LL}_{ab21}\end{array}\right\\}+\text{h.c.}\\\
&=-\frac{4G_{F}}{\sqrt{2}}\frac{\alpha}{2\pi
s^{2}_{w}}\sum\limits_{\begin{subarray}{c}\ell=\left\\{e,\mu,\tau\right\\}\\\
a,b=\left\\{1,2,3\right\\}\end{subarray}}U_{\ell a}^{*}U_{\ell
b}\left(\lambda_{c}X^{\ell}+\lambda_{t}X_{t}\right)\times\left\\{\begin{array}[]{r}-O^{A,\,L}_{ab21}\\\\[5.0pt]
O^{V,\,LL}_{ab21}\end{array}\right\\}+\text{h.c.}\,,\end{split}$ (3)

where $G_{F}$ denotes the Fermi constant, $\alpha$ is the fine-structure
constant, and $s_{w}$ is the sine of the weak-mixing angle. The elements of
the quark-mixing matrix are comprised in the parameters
$\lambda_{i}=V^{*}_{is}V_{id}$, while $U_{\ell a}$ is the
Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix. The short-distance physics is
described by the loop functions $X_{t}$ and $X^{\ell}$, where the dependence
on the $\tau$ mass introduces a small off-diagonal neutrino coupling, which we
ignore in the following. The operator induced in the SM case reads

$\displaystyle
O^{A,\,L}_{abij}=\frac{1}{2}(\overline{\nu}_{Ma}\gamma_{\mu}\gamma_{5}\,\,\nu_{Mb})(\overline{d}_{i}\,\gamma^{\mu}\,P_{L}\,{d}_{j})\,,$
$\displaystyle
O^{V,\,LL}_{\,abij}=(\overline{\nu}_{Da}\gamma_{\mu}\,P_{L}\,\nu_{Db})(\overline{d}_{i}\,\gamma^{\mu}\,P_{L}\,{d}_{j})\,.$
(4)

in the Majorana and Dirac case, respectively. We have intentionally decided to
directly separate the SM from NP contributions in Eq. (1) in order to present
the constraints on the NP effects more transparently.

We discuss the dimension-six operator basis for the Majorana and Dirac case in
Section 2.1 and Section 2.2, respectively. There are further possibilities,
e.g., mixed Dirac–Majorana cases or Dirac masses with lepton-number-violating
dimension-six interactions, but as we shall see the main effects can be
illustrated within the two aforementioned cases. The basis of operators can be
equally well be adjusted by trivial flavour renaming to study $b\to s/d\nu\nu$
or $c\to u\nu\nu$ transitions relevant, e.g, for $B\to K\nu\nu$ or rare charm
FCNC with decays to neutrinos.

###  2.1 Majorana-$\nu$ EFT with lepton-number violation

The dimension-six effective Lagrangian relevant for $d_{i}\rightarrow
d_{j}\nu\nu$ transitions for the case in which all neutrinos are Majorana
particles, can be written as

$\mathscr{L}^{(6)}_{d\to
d^{\prime}\nu\nu}\Big{|}_{\text{Majorana}}=\sum_{\begin{subarray}{c}I=\left\\{V,A\right\\}\\\
\tau=\left\\{L,R\right\\}\end{subarray}}\sum_{f}C^{I,\tau}_{f}O^{I,\tau}_{f}+\Biggl{(}\sum_{\begin{subarray}{c}I=\left\\{S,P,T\right\\}\end{subarray}}\\\
\sum_{f}C^{I,L}_{f}O^{I,L}_{f}+\text{h.c.}\Biggr{)}\,,$ (5)

where $I$ denotes the type of interaction, $I=\\{V,A,S,P,T\\}$, which stand
for $V$ for vector-, $A$ for axial-vector-, $S$ for scalar-, $P$ for
pseudoscalar-, and $T$ for tensor-type interactions. $\tau=\\{L,R\\}$
indicates the chirality of the quark currents, and $f$ comprises all neutrino
and quark flavour indices, $f=\\{abij\\}$. The full set of independent
operators read

$\displaystyle\begin{aligned}
O^{V,\,L}_{abij}&=\frac{1}{2}(\overline{\nu}_{Ma}\gamma_{\mu}\,\nu_{Mb})(\overline{d}_{i}\,\gamma^{\mu}\,P_{L}\,{d}_{j})\,,&O^{V,\,R}_{abij}&=\frac{1}{2}(\overline{\nu}_{Ma}\gamma_{\mu}\,\nu_{Mb})(\overline{d}_{i}\,\gamma^{\mu}\,P_{R}\,{d}_{j})\,,\\\
O^{A,\,L}_{abij}&=\tfrac{1}{2}(\overline{\nu}_{Ma}\gamma_{\mu}\gamma_{5}\,\,\nu_{Mb})(\overline{d}_{i}\,\gamma^{\mu}\,P_{L}\,{d}_{j})\,,&O^{A,\,R}_{abij}&=\frac{1}{2}(\overline{\nu}_{Ma}\gamma_{\mu}\,\gamma_{5}\,\nu_{Mb})(\overline{d}_{i}\,\gamma^{\mu}\,P_{R}\,{d}_{j})\,,\\\
O^{S,\,L}_{abij}&=\frac{1}{2}(\overline{\nu}_{Ma}\,\nu_{Mb})(\overline{d}_{i}\,P_{L}\,{d}_{j})\,,&O^{P,\,L}_{abij}&=\frac{1}{2}(\overline{\nu}_{Ma}\,i\gamma_{5}\,\nu_{Mb})(\overline{d}_{i}\,P_{L}\,{d}_{j})\,,\\\
O^{T,\,L}_{abij}&=\frac{1}{2}(\overline{\nu}_{Ma}\,\sigma_{\mu\nu}\,\nu_{Mb})(\overline{d}_{i}\,\sigma^{\mu\nu}\,P_{L}\,{d}_{j})\,,\end{aligned}$
(6)

with $i,j=1,2,3$.

Apart from the quark-flavour dependence, the operators $O^{V/A,\,L/R}$ are
self-adjoint after considering standard relations for four-component Majorana
fields [30]. This is not the case for the $O^{S/P/T,\,L}$ operators. We ensure
the hermiticity of the Lagrangian by imposing additional conditions on the
flavour structure of the Wilson coefficients. The hermiticity condition on the
Wilson coefficients $C^{V/A,\,L/R}_{f}$ is

$C^{V/A,\,L/R}_{abij}=\left(C^{V/A,\,L/R}_{baji}\right)^{*}\,.$ (7)

Additionally, there are the further conditions originating from standard
relations specific to four-component Majorana fields [30]

$C^{I,\,\tau}_{\,abij}=\eta\,C^{I,\,\tau}_{\,baij}\quad\quad\text{where}\quad\eta=\begin{cases}+1,&\text{for}\quad
I=A,\,S,\,P\\\ -1,&\text{for}\quad I=V,\,T\,.\end{cases}$ (8)

The SM contribution to $K^{+}\to\pi^{+}\nu\nu$ ($\bar{s}\to\bar{d}\nu\nu$
transition) in the Majorana case is induced by the $O^{A,\,L}_{abij}$
operators with $i=1$ (down quark) and $j=2$ (strange quark), as discussed in
Eq. (3). Note that imposing lepton-number conservation forbids the neutrinos
mass-term but also the $O^{S/P/T,\,L}$ interactions, which are thus lepton-
number violating.

For $s\to d\nu\nu$ transitions and the case of three Majorana neutrinos, the
number of independent interactions that are symmetric in the flavour of the
three Majorana neutrinos are parametrised by $6\cdot 6$ independent complex
Wilson coefficients contained in $C^{A,\,L/R}_{ab21}$, $C^{S/P,\,L}_{ab21}$,
and $C^{S/P,\,L}_{ab12}$, with $a\geq b$. In addition, the interactions that
are anti-symmetric in the neutrino flavours are parametrised by $3\cdot 4$
complex Wilson coefficients contained in $C^{V,\,L}_{ab21}$,
$C^{V,\,R}_{ab21}$, $C^{T,\,L}_{ab21}$, and $C^{T,\,L}_{ab12}$, with $a>b$.
The same number of $48$ independent complex Wilson coefficients contributes to
$b\to s\nu\nu$ decays, after the quark-flavour indices are appropriately
changed.

###  2.2 Dirac-$\nu$ EFT with lepton-number conservation

For the case of Dirac neutrinos, the dimension-six effective Lagrangian with
lepton-number conservation relevant for $d_{i}\rightarrow d_{j}\nu\bar{\nu}$
transitions can be written as

$\mathscr{L}^{(6)}_{d\to
d^{\prime}\nu\nu}\Big{|}_{\text{Dirac}}\\!\\!=\\!\\!\sum_{\tau,\tau^{\prime}=\left\\{L,R\right\\}}\\!\\!\sum_{f}C^{V,\,\tau\tau^{\prime}}_{f}O^{V,\,\tau\tau^{\prime}}_{f}+\sum_{f}\Bigl{(}C^{S,\,LL}_{f}O^{S,\,LL}_{f}+C^{S,\,LR}_{f}O^{S,\,LR}_{f}+C^{T,\,LL}_{f}O^{T,\,LL}_{f}+\text{h.c.}\Bigr{)},$
(9)

where $\tau,\tau^{\prime}=\\{L,R\\}$ indicate the chirality of the neutrino
and/or quark currents and $f$ comprises all neutrino and quark flavour
indices, $f=\\{abij\\}$. The independent, lepton-number conserving operators
read

$\displaystyle O^{V,\,LL}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\gamma_{\mu}\,P_{L}\,\nu_{Db})(\overline{d}_{i}\,\gamma^{\mu}\,P_{L}\,{d}_{j})\,,$
$\displaystyle O^{V,\,LR}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\gamma_{\mu}\,P_{L}\,\nu_{Db})(\overline{d}_{i}\,\gamma^{\mu}\,P_{R}\,{d}_{j})\,,$
(10) $\displaystyle O^{V,\,RL}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\gamma_{\mu}\,P_{R}\,\nu_{Db})(\overline{d}_{i}\,\gamma^{\mu}\,P_{L}\,{d}_{j})\,,$
$\displaystyle O^{V,\,RR}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\gamma_{\mu}\,P_{R}\,\nu_{Db})(\overline{d}_{i}\,\gamma^{\mu}\,P_{R}\,{d}_{j})\,,$
(11) $\displaystyle O^{S,\,LL}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\,P_{L}\,\nu_{Db})(\overline{d}_{i}\,P_{L}\,{d}_{j})\,,$
$\displaystyle O^{S,\,LR}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\,P_{L}\,\nu_{Db})(\overline{d}_{i}\,P_{R}\,{d}_{j})\,,$
(12) $\displaystyle O^{T,\,LL}_{\,abij}$
$\displaystyle=(\overline{\nu}_{Da}\,\sigma_{\mu\nu}\,P_{L}\,\nu_{Db})(\overline{d}_{i}\,\sigma^{\mu\nu}\,P_{L}\,{d}_{j})\,.$
(13)

In the Dirac case, the $O^{V,\,\tau\tau^{\prime}}$ operators are all self-
adjoint if their flavour structure is neglected. The hermiticity condition on
their Wilson coefficients reads

$C^{V,\,\tau\tau^{\prime}}_{abij}=\left(C^{V,\,\tau\tau^{\prime}}_{baji}\right)^{*}\,.$
(14)

In the Dirac basis, the SM contribution to $K^{+}\to\pi^{+}\bar{\nu}\nu$ is
induced by the operators $O^{V,\,LL}_{abij}$ with $i=1$ (down quark) and $j=2$
(strange quark) and discussed in Eq. (3). In the strict SM case of three
massless neutrinos, neutrinos can be described by three left-handed fields
without right-handed partners. In this case, the scalar and tensor operators
vanish by construction. They do not, however, vanish if the small neutrino
masses are due to Dirac masses for which right-handed neutrinos are required.
In this case, scalar and tensor interactions are possible without the need for
lepton-number violation in contrast to the Majorana case.

For $s\to d\nu\bar{\nu}$ transitions and the case of three Dirac neutrinos,
the interactions are parametrised by a total of $9\cdot 10$ independent
complex Wilson coefficients $C^{V,\,\tau\tau^{\prime}}_{ab21}$,
$C^{S,\,LL/LR}_{ab21}$, $C^{S,\,LL/LR}_{ab12}$, $C^{T,\,LL}_{ab21}$,
$C^{T,\,LL}_{ab12}$, where $\tau,\tau^{\prime}\in\left\\{L,R\right\\}$. The
same number of $90$ complex Wilson coefficients contributes to $b\to s\nu\nu$
decays, after the quark flavour indices are appropriately changed.

##  3 Rates and distributions for $K^{+}\to\pi^{+}\nu\bar{\nu}$ at NA62

In this section, we derive the decay rates and the resulting distributions for
the decay $K^{+}\rightarrow\pi^{+}\nu\nu$ for the cases of Majorana and Dirac
neutrinos. We also discuss the treatment of theory and experimental
uncertainties used throughout the analysis. Finally, we describe the
statistical procedure that we shall use to constrain the parameter-space of
the effective theories based on current and expected data from NA62 and HIKE,
respectively.

###  3.1 Decay rates

To derive the differential decay rates for
$K^{+}\rightarrow\pi^{+}\nu\overset{\scriptscriptstyle(-)}{\nu}$ we evaluate
the matrix element

${\cal
M}=\Braket{\pi^{+}(p_{\pi})\overset{\scriptscriptstyle(-)}{\nu}_{\\!\\!\\!a}\nu_{b}}{H_{\mathrm{eff}}}{K^{+}(p_{K})}\,,$

where the effective Hamiltonian $H_{\mathrm{eff}}$ follows either from the
Lagrangian in Eq. (6) or from the one in Eq. (10). To be precise, for the
Dirac case we compute the amplitude for the transition
$K^{+}\rightarrow\pi^{+}\bar{\nu}_{a}\nu_{b}$, while for the Majorana case the
one for $K^{+}\rightarrow\pi^{+}\nu_{a}\nu_{b}$. We will mention differences
between the Majorana and the Dirac computation when they arise.

The amplitude factorizes into a hadronic part
$\braket{\pi^{+}}{O^{\mathrm{had}}_{sd}}{K^{+}}$ and a leptonic part
$\braket{\overset{\scriptscriptstyle(-)}{\nu}_{\\!\\!\\!a}\nu_{b}}{O^{\mathrm{lep}}_{ab}}{0}$.
The leptonic part is obtained using standard techniques; in particular we
employ Feynman rules for Majorana fields [31] to treat the case of massive
Majorana neutrinos. The hadronic parts of the decay are determined by three
form factors [32, 33, 11]

$\displaystyle\Braket{\pi^{+}}{\bar{d}\gamma^{\mu}s}{K^{+}}$
$\displaystyle=(p_{K}+p_{\pi})^{\mu}f^{V}(q^{2})+\frac{m_{K}^{2}-m_{\pi}^{2}}{q^{2}}(p_{K}-p_{\pi})^{\mu}\left(f^{S}(q^{2})-f^{V}(q^{2})\right)\,,$
(15) $\displaystyle\Braket{\pi^{+}}{\bar{d}s}{K^{+}}$
$\displaystyle=\frac{m_{K}^{2}-m_{\pi}^{2}}{m_{s}-m_{d}}f^{S}(q^{2})\,,$ (16)
$\displaystyle\Braket{\pi^{+}}{\bar{d}\sigma^{\mu\nu}s}{K^{+}}$
$\displaystyle=\left(p^{\mu}_{\pi}p^{\nu}_{K}-p^{\nu}_{\pi}p^{\mu}_{K}\right)\frac{f^{T}(q^{2})}{m_{K}+m_{\pi}}\,,$
(17)

with $m_{p}$ denoting the masses of mesons and quarks, $p_{\pi}$ and $p_{K}$
the $\pi^{+}$ and $K^{+}$ momentum, respectively, and
$q^{2}\equiv(p_{K}-p_{\pi})^{2}$. To obtain the tensor form factor,
$f^{T}(q^{2})$, we use isospin symmetry to relate it to the neutral form
factor, finding $f^{T}(q^{2})=2f^{T}_{0}(q^{2})$, with
$\Braket{\pi^{0}}{\bar{s}\sigma^{\mu\nu}d}{K^{0}}=\left(p^{\mu}_{\pi}p^{\nu}_{K}-p^{\nu}_{\pi}p^{\mu}_{K}\right)\sqrt{2}f^{T}_{0}(q^{2})/(m_{K}+m_{\pi})$
evaluated in Ref. [33]. We drop the subscript ‘‘${K^{+}\pi^{+}}$’’ for
readability. Note, that the parity violating axial parts of the hadronic
matrix elements do not contribute to the decay due to the parity conservation
of QCD. The $q^{2}$-dependence of the remnant scalar form-factor functions
reads [11, 33]

$\displaystyle
f^{V}(q^{2})=f^{+}_{0}\left(1+\lambda^{\prime}_{+}\frac{q^{2}}{m_{\pi}^{2}}+\lambda^{\prime\prime}_{+}\frac{q^{4}}{m_{\pi}^{4}}\right),~{}f^{S}(q^{2})=f^{+}_{0}\left(1+\lambda_{0}\frac{q^{2}}{m_{\pi}^{2}}\right),~{}f^{T}(q^{2})=\frac{2f^{T}_{0}(0)}{1-\lambda_{T}\,q^{2}}\,.$
(18)

Table 1 contains the numerical input for the form factors. For notational
brevity we also define

$\displaystyle F^{S}(q^{2})$
$\displaystyle\equiv\frac{m_{K}^{2}-m_{\pi}^{2}}{m_{s}-m_{d}}f^{S}(q^{2})\,,$
$\displaystyle F^{T}(q^{2})$
$\displaystyle\equiv\frac{f^{T}(q^{2})}{m_{K}+m_{\pi}}\,.$ (19)

Table 1: Numerical values of the form factor parameters used. The values for
the scalar and vector form factor are taken from Ref. [11, 32], while the
parameters of the tensor form factor are taken from Ref. [33].
Form factor | Parameters
---|---
$f_{S}$ | $f^{+}_{0}=$ | $0.9778$
$\lambda_{0}=$ | $(13.38\pm 1.19)\cdot 10^{-3}$
$f_{V}$ | $f^{+}_{0}=$ | $0.9778$
$\lambda^{\prime}_{+}=$ | $(24.82\pm 1.10)\cdot 10^{-3}$
$\lambda^{\prime\prime}_{+}=$ | $(1.64\pm 0.44)\cdot 10^{-3}$
$f_{T}$ | $f^{T}_{0}=$ | $0.417(14_{\text{stat}})(5_{\text{syst}})$
$s_{T}=$ | $1.10(8_{\text{stat}})(11_{\text{syst}})\,\mathrm{GeV}^{-2}$

The differential partial decay rate for the three-body decay is then given by
[18]

$\mathrm{d}\Gamma\left(K^{+}\to\pi^{+}\overset{\scriptscriptstyle(-)}{\nu}_{\\!\\!\\!a}\nu_{b}\right)=\frac{(2\pi)^{4}}{2E_{K}}\left|\overline{{\cal
M}}\right|^{2}\frac{\mathrm{d}^{3}\vec{p}_{\pi}}{(2\pi)^{3}E_{\pi}}\frac{\mathrm{d}^{3}\vec{p}_{\nu,a}}{(2\pi)^{3}E_{\nu,a}}\frac{\mathrm{d}^{3}\vec{p}_{\nu,b}}{(2\pi)^{3}E_{\nu,b}}\,,$
(20)

where $\left|\overline{{\cal M}}\right|$ indicates that we have summed over
the polarisations of the final-state neutrinos. Starting from this equation,
one can derive the differential distribution
$\mathrm{d}^{2}\Gamma_{ab}/\mathrm{d}q^{2}\mathrm{d}k^{2}$ in terms of the
Lorentz invariant parameters $q^{2}\equiv(p_{\nu,a}+p_{\nu,b})^{2}$ and
$k^{2}\equiv(p_{\pi}+p_{\nu,a})^{2}$ to obtain the characteristic for three-
body decays Dalitz plot. The corresponding expressions in the Majorana and
Dirac case are provided in Appendix A.

However, NA62 measures the differential branching ratio
$\mathrm{d}\mathrm{Br}/\mathrm{d}|\vec{p}_{\pi}|$ in the lab frame. To arrive
at this expression, we recursively write the three-body phase-space
integration into a convolution of two two-particle phase-space integrations.
This introduces an additional integration over the variable
$q^{2}=(p_{\nu,a}+p_{\nu,b})^{2}$ but allows one to trivially perform the
integration over both neutrino momenta. After the integration over the angular
part of the pion momentum in the lab frame, we are left with the double
differential decay rate
$\mathrm{d}^{2}\Gamma_{ab}/\mathrm{d}|\vec{p}_{\pi}|\mathrm{d}q^{2}$ in the
lab frame.

For the case of a decay to two massless neutrinos, $\nu_{a}$ and $\nu_{b}$,
the double differential partial decay width of $K^{+}\to\pi^{+}\nu_{a}\nu_{b}$
in the Majorana basis then reads

$\begin{split}\frac{\mathrm{d}^{2}\Gamma_{ab}}{\mathrm{d}|\vec{p}_{\pi}|\mathrm{d}q^{2}}\Bigg{|}_{\text{Majorana}}&=\frac{|\vec{p}_{\pi}|}{6144\pi^{3}E_{K}E_{\pi}|\vec{p}_{K}|\left(1+\delta_{ab}\right)}\Theta(|\vec{p}_{\pi}|,q^{2})\Bigl{|}_{m_{a}=m_{b}=0}\Bigg{[}\\\
&+8\lambda_{K\pi
q^{2}}\left|f^{V}(q^{2})\right|^{2}\left(\left|C^{A,\,L,\,\text{SM}}_{ab12}+C^{A,\,L}_{ab12}+C^{A,\,R}_{ab12}\right|^{2}+\left|C^{V,\,L}_{ab12}+C^{V,\,R}_{ab12}\right|^{2}\right)\\\
&+12q^{2}\left|F^{S}(q^{2})\right|^{2}\Big{(}\left|C^{P,\,L}_{ab12}+C^{P,\,L*}_{ab21}\right|^{2}+\left|C^{S,\,L}_{ab12}+C^{S,\,L*}_{ab21}\right|^{2}\Big{)}\\\
&+3q^{2}\lambda_{K\pi
q^{2}}\left|F^{T}(q^{2})\right|^{2}\Big{(}\left|C^{T,\,L}_{ab12}-C^{T,\,L*}_{ab21}\right|^{2}\Big{)}\Bigg{]}\end{split}$
(21)

with $\lambda_{K\pi q^{2}}\equiv\lambda(m_{K}^{2},m_{\pi}^{2},q^{2})$ and the
Källén function $\lambda(x,y,z)\equiv x^{2}+y^{2}+z^{2}-2xy-2yz-2zx$. The
additional factor $1/(1+\delta_{ab})$ accounts for the reduction of the phase
space in the case of identical particles in the final state. We stress that in
Eq. (21) we have already used the hermiticity and Majorana conditions to
obtain all contributions to the neutrino final state $\nu_{a}\nu_{b}$.
Therefore, the total rate corresponds to the ‘‘constrained’’ sum of channels
$\sum_{\begin{subarray}{c}a,b\\\ a\leq b\end{subarray}}$ of all independent
final states. In Eq. (21), $\Theta(|\vec{p}_{\pi}|,q^{2})$ is a shorthand for
the theta-functions

$\begin{split}\Theta(|\vec{p}_{\pi}|,q^{2})\equiv~{}&\theta(E_{K}-E_{\pi})\theta\left(\sqrt{q^{2}}-m_{a}-m_{b}\right)\\\
\times~{}&\theta\left(1-\frac{m_{K}^{2}+m_{\pi}^{2}-2E_{K}E_{\pi}-q^{2}}{2|\vec{p}_{K}||\vec{p}_{\pi}|}\right)\theta\left(1+\frac{m_{K}^{2}+m_{\pi}^{2}-2E_{K}E_{\pi}-q^{2}}{2|\vec{p}_{K}||\vec{p}_{\pi}|}\right)\,.\end{split}$
(22)

Analogously, we compute the differential partial decay width of
$K^{+}\to\pi^{+}\bar{\nu}_{a}\nu_{b}$ for the massless Dirac-neutrino case

$\begin{split}\frac{\mathrm{d}^{2}\Gamma_{ab}}{\mathrm{d}|\vec{p}_{\pi}|\mathrm{d}q^{2}}\Bigg{|}_{\text{Dirac}}&=\frac{|\vec{p}_{\pi}|}{3072\pi^{3}E_{K}E_{\pi}|\vec{p}_{K}|}\Theta(|\vec{p}_{\pi}|,q^{2})\Bigl{|}_{m_{a}=m_{b}=0}\Bigg{[}\\\
&+2\lambda_{K\pi
q^{2}}|f^{V}(q^{2})|^{2}\Bigl{(}\bigl{|}C^{V,\,LL,\,\text{SM}}_{ab12}+C^{V,\,LL}_{ab12}+C^{V,\,LR}_{ab12}\bigr{|}^{2}+\bigl{|}C^{V,\,RL}_{ab12}+C^{V,\,RR}_{ab12}\bigr{|}^{2}\Bigr{)}\\\
&+3q^{2}|F^{S}(q^{2})|^{2}\Bigl{(}\bigl{|}C^{S,\,LL}_{ab12}+C^{S,\,LR}_{ab12}\bigr{|}^{2}+\bigl{|}C^{S,\,LL}_{ba21}+C^{S,\,LR}_{ba21}\bigr{|}^{2}\Bigr{)}\\\
&+q^{2}\lambda_{K\pi
q^{2}}|F^{T}(q^{2})|^{2}\Bigl{(}\bigl{|}C^{T,\,LL}_{ab12}\bigr{|}^{2}+\bigl{|}C^{T,\,LL}_{ba21}\bigr{|}^{2}\Bigr{)}\Bigg{]}\end{split}$
(23)

Similarly to the Majorana case, also here we have used the hermiticity
conditions on the Wilson coefficients in order to include all contributions to
the $\bar{\nu}_{a}\nu_{b}$ final state. However, contrarily to the Majorana
case, particles are not the same as antiparticles and thus the total rate is
given by the unconstrained sum of channels $\sum_{a,b}$.

The total branching ratio for
$K^{+}\to\pi^{+}\nu\overset{\scriptscriptstyle(-)}{\nu}$ is then the sum of
the different neutrino channels

$\mathrm{Br}\left(K^{+}\to\pi^{+}\nu\overset{\scriptscriptstyle(-)}{\nu}\right)=\Gamma_{\text{tot}}^{-1}\int_{\Omega(|\vec{p}_{\pi}|,q^{2})}\mathrm{d}|\vec{p}_{\pi}|\mathrm{d}q^{2}\overline{\sum_{a,b}}\frac{\mathrm{d}^{2}\Gamma_{ab}}{\mathrm{d}|\vec{p}_{\pi}|\mathrm{d}q^{2}}\,,$
(24)

where $\Omega(|\vec{p}_{\pi}|,q^{2})$ denotes the signal region of the NA62
experiment, $\Gamma_{\text{tot}}$ the total decay width in the lab frame, and
the overline in $\overline{\sum}$ indicates that for the Majorana case the sum
is constrained, i.e., $a\leq b$, while for the Dirac case it is unconstrained.

#### Uncertainties in theory predictions

There are two main sources of uncertainties in the theory predictions for the
$K\to\pi\nu\nu$ rates and distributions: (i) uncertainties associated to the
SM Wilson coefficient and (ii) uncertainties associated to the form-factors.

For (i), i.e., SM Wilson coefficient, we follow Ref. [17]. We include all
state-of-the-art perturbative corrections and include in our analysis both
pure theory uncertainties, associated to residual scale dependencies, and
parametric uncertainties. At the current level of precision, the dominant
uncertainties are parametric, mainly originating from the CKM input. All these
uncertainties are due to short-distance dynamics and thus do not depend on
$q^{2}$.

Figure 1:  The $q^{2}$ dependence of the form-factor functions and their
uncertainties as given in Table 1. The coloured bands correspond to the
$1\sigma$ bands for the vector, scalar and tensor form-factors. The lower
panel shows the dependence of the relative uncertainty $\Delta f/f$ of each
form-factor function.

For (ii), i.e., the form factors, we use their numerical values and
uncertainties as published in the original literature and summarised in Table
1. The numerical evaluation of the form factors shows that their relative
uncertainties are roughly constant within the $q^{2}$ signal regions of NA62.
This is illustrated in Figure 1, where we show the $q^{2}$ dependence of the
three form-factor functions and their uncertainties in the NA62 signal
regions. The lower panel shows the $q^{2}$ dependence of their relative
uncertainties ${\Delta f}/{f}$. Since these uncertainties are roughly flat in
$q^{2}$, we can consider them as global relative uncertainties, which speeds
up the evaluation of the profile likelihood ratio test. This simplifies
accounting for them in the likelihood evaluation. Concretely, we use
throughout the whole signal region an overall relative error of $0.4\%$,
$0.8\%$, and $4\%$ for the vector, scalar, and tensor form factor,
respectively, cf., Figure 1.

###  3.2 Statistical treatment for EFT fit

To quantify the agreement between the experimental data and the different NP
scenarios, we perform a statistical test using a profile likelihood ratio as
test statistic that is sampled in a full frequentist manner. The number of
signal events, $s_{i}$, in a category $i$ is given as

$\displaystyle
s_{i}=\frac{\mathrm{BR}_{K\rightarrow\pi\nu\nu}^{i}(\\{C_{\mathrm{NP}}\\})}{{\mathrm{S}ES}_{i}}\,,$
(25)

where $\mathrm{BR}_{K\rightarrow\pi\nu\nu}^{i}(\\{C_{\mathrm{NP}}\\})$ denotes
the branching ratio of $K^{+}\rightarrow\pi^{+}\nu\nu$ within the category $i$
for a set of Wilson coefficients $\\{C_{\mathrm{NP}}\\}$ corresponding to the
NP scenario under consideration. Further, the variable ${\mathrm{S}ES}_{i}$
denotes the single-event sensitivity of category $i$ taken from Ref. [34].
Note that in Ref. [34, 2], the ${\mathrm{S}ES}_{i}$ are defined such that
multiplied with the total branching ratio one obtains the number of expected
events in bin $i$. This is in contrast to our definition in Eq. (25) where we
use the branching ratio within the corresponding categories. We, thus, remove
the phase-space information from the original definition of the
${\mathrm{S}ES}_{i}$, cf. Table 2.

The number of expected background events, $b_{i}$, and observed data, $n_{i}$,
in each category are taken from the NA62 measurement [2]. The events in the
different categories are uncorrelated, which implies that the total likelihood
is a product of Poisson distributions in the number of signal plus background
events $s_{i}+b_{i}$. We further neglect possible correlations between
experimental parameters like the ${\mathrm{S}ES}_{i}$, which—if non-zero—are
not provided by the NA62 collaboration. We account for theory and experimental
statistical uncertainties in the likelihood by including the corresponding
parameters as global observables associated to auxiliary measurements. The
uncertainties of the parameters for the theory predictions described in
Section 3.1 are all assumed to be Gaussian even though some of them, e.g.,
scale dependencies are of non-statistical nature. However, the numerical
impact of this assumption on the analysis is negligible. The sensitivities,
${\mathrm{S}ES}_{i}$, are constrained to follow uncorrelated Gaussian terms
while the number of expected background events follows a Poisson distribution.

Table 2: The single-event sensitivities, ${\mathrm{S}ES}_{i}$, for the
respective categories provided by the NA62 collaboration [35, 1, 2, 34], as
well as the corresponding sensitivities with the phase-space information
removed. The definition of each category is given in the original NA62
publications.
Category | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
${\mathrm{S}ES}/10^{-10}$ from [35, 1, 2, 34] | $3.15$ | $0.39$ | $0.54$ | $1.48$ | $0.59$ | $0.55$ | $0.63$ | $1.22$ | $1.75$
${\mathrm{S}ES}/10^{-12}$ w/o kinematic cuts | 7$6.01$ | $8.45$ | 1$6.03$ | $8.04$ | $3.20$ | $2.99$ | $3.42$ | $4.86$ | $6.97$

Let us now denote the total likelihood with the shorthand notation
$\mathcal{L}(x|\\{C_{\mathrm{NP}}\\},\nu)$, with $x$ being the outcome, i.e.,
the observed data, $\\{C_{\mathrm{NP}}\\}$ the parameters of interest, i.e.,
the NP Wilson coefficients, and $\nu$ the nuisance parameters. Under the
assumption that the NP scenario under consideration is realised in nature, one
can determine the confidence intervals on the model parameters
$\\{C_{\mathrm{NP}}\\}$ at a given confidence level using the profile
likelihood ratio $\lambda$ as a test statistic

$\displaystyle\lambda_{\\{C_{\mathrm{NP}}\\}}=\frac{\mathcal{L}(x|\\{C_{\mathrm{NP}}\\},\hat{\hat{\nu}}(\\{C_{\mathrm{NP}}\\}))}{\mathcal{L}(x|\\{\hat{C}_{\mathrm{NP}}\\},\hat{\nu})}\,.$
(26)

The parameters with a single hat are the parameter values that maximise the
likelihood. The value for the nuisance parameters $\nu$ that maximise the
likelihood for given values of the parameters of interest
$\\{C_{\mathrm{NP}}\\}$ are denoted with a double-hat. For given parameters of
interest $\\{C_{\mathrm{NP}}\\}$, we calculate the corresponding $p$-value via

$\displaystyle
p_{\\{C_{\mathrm{NP}}\\}}=\\!\\!\\!\int_{\lambda^{\mathrm{obs}}_{\\{C_{\mathrm{NP}}\\}}}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\lambda_{\\{C_{\mathrm{NP}}\\}}\,f(\lambda_{\\{C_{\mathrm{NP}}\\}}|\\{C_{\mathrm{NP}}\\},\nu)\,,$
(27)

where $f$ denotes the test-statistic distribution that we sample by means of a
Monte-Carlo method under the assumption of $\\{C_{\mathrm{NP}}\\}$, and with
$\lambda^{\mathrm{obs}}_{\\{C_{\mathrm{NP}}\\}}$ the observed value of the
test statistic. In order to determine the upper limit on a single Wilson
coefficient $C_{\mathrm{NP}}$ at a confidence level of $100(1-\alpha)\%$, we
solve for $p_{C_{\mathrm{NP}}}=\alpha$. For the procedure on how to estimate
two-dimensional confidence regions from the profile likelihood ratio see,
e.g., Ref. [36]. For further references, see Refs. [37, 18]

##  4 New Physics sensitivities

Figure 2:  Distributions of $\mathrm{Br}(K^{+}\to\pi^{+}\nu\nu)$ for different
NP scenarios containing three massless Majorana or Dirac neutrinos. The left
and middle panel shows differential $q^{2}$ distributions integrated over the
whole physical phase space for $|\vec{p}_{\pi}|$ and over $|\vec{p}_{\pi}|$
signal region of NA62, respectively. The right panel shows the lab frame
$|\vec{p}_{\pi}|$ distribution inside the NA62 $q^{2}$ signal regions as
currently done at NA62. All NP Wilson coefficients that affect the rates are
displayed in the legend. In each panel we use the values for the given Wilson
coefficients that saturate the bound at $90\%$ CL. Therefore, there is no
distinction between $a=b$ and $a\neq b$ cases. The gray, dashed line
corresponds to the SM case. The dashed, vertical lines indicate the different
NA62 signal regions.

In this Section, we use the NP distributions from Section 3.1 together with
the binned likelihood from Section 3.2 to determine NA62’s sensitivity on
constraining the independent NP operators given in Section 2.1 and 2.2 for
Majorana and Dirac neutrinos, respectively. We shall discuss the sensitivity
based on the current experimental data of NA62 [2] and also provide estimates
for the future sensitivity based on the prospects in future experiments like
HIKE. In Section 4.1 we switch on one NP operator at a time, while in Section
4.2 we also show correlations by switching on pairs of operators. The limits
on all Wilson coefficients are always determined in addition to the Standard
Model contribution, which is always included.

###  4.1 Single-operator fits

For the single-operator fits we discuss three different scenarios:

  1. (i)

Dirac effective theory with three light neutrinos, i.e., $a,b=1,2,3$
with $m_{\nu,1}\\!\approx\\!m_{\nu,2}\\!\approx\\!m_{\nu,3}\\!\approx\\!0$
(operator basis of Section 2.2).

  2. (ii)

Majorana effective theory with three light neutrinos, i.e., $a,b=1,2,3$
with $m_{\nu,1}\\!\approx\\!m_{\nu,2}\\!\approx\\!m_{\nu,3}\\!\approx\\!0$
(operator basis of Section 2.1).

  3. (iii)

Majorana effective theory with three light neutrinos and one additional
sterile neutrino with mass $m_{\nu,4}$, i.e., $a,b=1,2,3,4$, (operator basis
of Section 2.1).

Since in this Section we switch on one operator at a time, the bounds are only
sensitive to the absolute value of the respective Wilson coefficients and
insensitive to their CP violating phases. The exceptions are NP contributions
that induce the same operator as the SM, i.e., $O^{A,\,L}_{aa12}$ and
$O^{V,\,LL}_{aa12}$, for which the interference term with the SM depends on
the relative phase between the NP and the SM Wilson coefficient. Since
$K^{+}\to\pi^{+}\nu\nu$ is a CP conserving decay we choose for concreteness to
align the NP with SM phase in the numerical analysis that follows.

Table 3:  Current and estimated future experimental lower limit at $90\%$ CL
for the different NP operators in the case of three massless Majorana
neutrinos. A single operator is turned on at a time. The case labelled by
$\sum a=b$ corresponds to taking neutrino-flavour universal Wilson
coefficients, whereas in the case labelled by $a=b$ a single Wilson
coefficient is switched on. Note that in case of interference with the SM
Wilson coefficient, the CP phase of the NP Wilson coefficient has been aligned
to the SM. See text for the details on the different future scenarios.
$\nu$-Majorana EFT |  |  current |  future |  future[$q^{2}$ bins]
---|---|---|---|---
$1/\sqrt{\left|C^{V,\,L/R}_{ab12}\right|}$ |  $a\neq b$ |  $7.6\cdot 10^{1}\,\text{TeV}$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $1.2\cdot 10^{2}\,\text{TeV}$
$1/\sqrt{\left|C^{A,\,L/R}_{ab12}\right|}$ |  $\sum a=b$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $2.8\cdot 10^{2}\,\text{TeV}$ |  $2.6\cdot 10^{2}\,\text{TeV}$
$a=b$ |  $8.1\cdot 10^{1}\,\text{TeV}$ |  $1.7\cdot 10^{2}\,\text{TeV}$ |  $1.6\cdot 10^{2}\,\text{TeV}$
$a\neq b$ |  $7.6\cdot 10^{1}\,\text{TeV}$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $1.2\cdot 10^{2}\,\text{TeV}$
$1/\sqrt{\left|C^{S,\,L}_{ab12/21}\right|}$ |  $\sum a=b$ |  $1.6\cdot 10^{2}\,\text{TeV}$ |  $2.4\cdot 10^{2}\,\text{TeV}$ |  $2.5\cdot 10^{2}\,\text{TeV}$
$a=b$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$
$a\neq b$ |  $1.4\cdot 10^{2}\,\text{TeV}$ |  $2.2\cdot 10^{2}\,\text{TeV}$ |  $2.3\cdot 10^{2}\,\text{TeV}$
$1/\sqrt{\left|C^{P,\,L}_{ab12/21}\right|}$ |  $\sum a=b$ |  $1.6\cdot 10^{2}\,\text{TeV}$ |  $2.4\cdot 10^{2}\,\text{TeV}$ |  $2.5\cdot 10^{2}\,\text{TeV}$
$a=b$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$
$a\neq b$ |  $1.4\cdot 10^{2}\,\text{TeV}$ |  $2.2\cdot 10^{2}\,\text{TeV}$ |  $2.3\cdot 10^{2}\,\text{TeV}$
$1/\sqrt{\left|C^{T,\,L}_{ab12/21}\right|}$ |  $a\neq b$ |  $3.0\cdot 10^{1}\,\text{TeV}$ |  $4.7\cdot 10^{1}\,\text{TeV}$ |  $4.7\cdot 10^{1}\,\text{TeV}$
Table 4:  Current and estimated future experimental lower limit at $90\%$ CL
for the different NP operators in the case of three massless Dirac neutrinos.
A single operator is turned on at a time. The case labelled by $\sum a=b$
corresponds to taking neutrino-flavour universal Wilson coefficients, whereas
in the case labelled by $a=b$ a single Wilson coefficient is switched on. Note
that in case of interference with the SM Wilson coefficient, the CP phase of
the NP Wilson coefficient has been aligned to the SM. See text for the details
on the different future scenarios.
$\nu$-Dirac EFT |  |  current |  future |  future[$q^{2}$ bins]
---|---|---|---|---
$1/\sqrt{\left|C^{V,L\tau}_{ab12}\right|}$ |  $\sum a=b$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $2.8\cdot 10^{2}\,\text{TeV}$ |  $2.6\cdot 10^{2}\,\text{TeV}$
$a=b$ |  $8.2\cdot 10^{1}\,\text{TeV}$ |  $1.7\cdot 10^{2}\,\text{TeV}$ |  $1.6\cdot 10^{2}\,\text{TeV}$
$a\neq b$ |  $6.3\cdot 10^{1}\,\text{TeV}$ |  $1.0\cdot 10^{2}\,\text{TeV}$ |  $9.6\cdot 10^{1}\,\text{TeV}$
$1/\sqrt{\left|C^{V,R\tau}_{ab12}\right|}$ |  $\sum a=b$ |  $8.4\cdot 10^{1}\,\text{TeV}$ |  $1.3\cdot 10^{2}\,\text{TeV}$ |  $1.3\cdot 10^{2}\,\text{TeV}$
$a=b$ |  $6.3\cdot 10^{1}\,\text{TeV}$ |  $1.0\cdot 10^{2}\,\text{TeV}$ |  $9.6\cdot 10^{1}\,\text{TeV}$
$a\neq b$ |  $6.3\cdot 10^{1}\,\text{TeV}$ |  $1.0\cdot 10^{2}\,\text{TeV}$ |  $9.6\cdot 10^{1}\,\text{TeV}$
$1/\sqrt{\left|C^{S,L\tau}_{ab12/21}\right|}$ |  $\sum a=b$ |  $1.6\cdot 10^{2}\,\text{TeV}$ |  $2.5\cdot 10^{2}\,\text{TeV}$ |  $2.5\cdot 10^{2}\,\text{TeV}$
$a=b$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$
$a\neq b$ |  $1.2\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$ |  $1.9\cdot 10^{2}\,\text{TeV}$
$1/\sqrt{\left|C^{T,LL}_{ab12/21}\right|}$ |  $\sum a=b$ |  $3.5\cdot 10^{1}\,\text{TeV}$ |  $5.7\cdot 10^{1}\,\text{TeV}$ |  $5.7\cdot 10^{1}\,\text{TeV}$
$a=b$ |  $2.7\cdot 10^{1}\,\text{TeV}$ |  $4.2\cdot 10^{1}\,\text{TeV}$ |  $4.2\cdot 10^{1}\,\text{TeV}$
$a\neq b$ |  $2.7\cdot 10^{1}\,\text{TeV}$ |  $4.2\cdot 10^{1}\,\text{TeV}$ |  $4.2\cdot 10^{1}\,\text{TeV}$

Even though the cases (i) and (ii) are physically distinguishable due to the
different number of degrees of freedom there will always be a scenario in (i)
that leads to the same NP distribution as a scenario in (ii) in the limit of
$m_{\nu}\to 0$ if a single NP operator is switched on. This is best seen by
comparing the Majorana to the Dirac distributions for massless neutrinos in
Eqs. (21) and (23) and is illustrated also in Figure 2, where we present all
the independent differential distributions possible for the case (i) and (ii).
Throughout Figure 2, the value of the Wilson coefficients is chosen to
saturate the current experimental limit at $90\%$ Confidence Level (CL). The
limits are determined from the prescription of Section 3.2.

Figure 3:  Distributions of $\mathrm{Br}(K^{+}\to\pi^{+}\nu\nu)$ for different
NP scenarios containing an additional fourth massive Majorana neutrino with
mass, $m_{\nu,4}=50\,\text{MeV}$. Similarly to Figure 2, the left and middle
panel shows differential $q^{2}$ distributions integrated over the whole
physical phase space for $|\vec{p}_{\pi}|$ and over the $|\vec{p}_{\pi}|$
signal region of NA62, respectively. The right panel shows the lab frame
$|\vec{p}_{\pi}|$ distribution inside the NA62 $q^{2}$ signal regions as
currently done at NA62. All NP Wilson coefficients that affect the rates are
displayed in the legend. In each panel we use the values for the given Wilson
coefficients that saturate the bound at $90\%$ CL. The gray, dashed line
corresponds to the SM case. The dashed, vertical lines indicate the different
NA62 signal regions.

In the left panel and middle panel of Figure 2, we show the different
$\mathrm{dBr}(K^{+}\rightarrow\pi^{+}\nu\nu)/\mathrm{d}q^{2}$ distributions.
The difference between left and middle panel is that the former shows the
distributions integrated over the whole pion-momentum, $|\vec{p}_{\pi}|$,
phase space while in the latter we have integrated over the pion-momentum
signal region of NA62. A comparison of the different shapes shows how a
differential measurement can be used to distinguish different operators. The
right panel shows the distributions
$\mathrm{dBr}(K^{+}\rightarrow\pi^{+}\nu\nu)/\mathrm{d}|\vec{p}_{\pi}|$ in the
NA62 lab frame after integrating over the $q^{2}$ signal region of NA62. This
corresponds exactly to the categories 3–8 currently employed by NA62 [2]. This
plot also makes it evident that a binning in the missing-momentum square,
$q^{2}$, is able to distinguish between the different NP scenarios within the
Majorana or the Dirac basis but not between the Majorana or Dirac nature of
the neutrinos. In Table 3 and 4 we collect the constraints at $90\%$ CL on the
different Wilson coefficients for the case of three Majorana neutrinos, i.e.,
scenario (i), and three Dirac neutrinos, i.e, scenario (ii), respectively.

Additionally, we perform first estimates of the future sensitivity reach based
on the proposed HIKE experiment. To this end we use the expected improvements
highlighted in the letter of intent for HIKE [38, 19] and perform a naive
upscaling of the luminosity such that the expected number of signal events
reach 400–500 while also reducing the number of background events to ensure
that the intended signal to background ratio of $B/S\sim 0.5\text{--}0.7$ is
achieved.111Ref. [19] and private communication with NA62. Apart from this, we
assume the same experimental setup as NA62. We take the expected data to be
the sum of the expected SM contribution plus the upscaled background of NA62
after the aforementioned reduction. The resulting limits are shown in the
column labeled ‘‘future’’ in Table 3 and 4 for the Majorana and Dirac case,
respectively.

For the NP scenarios (i) and (ii), we also estimate the expected limits on the
Wilson coefficients by choosing an alternative categorisation of the events.
For this we integrate over pion-momentum signal region, but choose to bin the
data in the current low-$q^{2}$ and high-$q^{2}$ signal region NA62 in two and
four equidistant $q^{2}$-bins, respectively. To achieve this, we rescale the
${\mathrm{S}ES}$ of the 2018-NEWCOL sample (category 3-8) under the assumption
that the efficiencies entering the ${\mathrm{S}ES}$ and the background
distributions are flat in $q^{2}$. The expected data is taken to be the sum of
the expected SM contribution plus the upscaled background. The estimated
sensitivities to the NP operators are shown in the right column of Table 3 and
4 for the Majorana and Dirac case, respectively. In both future scenarios, we
observe an increase of approximately $50\%$–$60\%$ for all NP operators apart
from the diagonal vector ones that interfere with the Standard Model. For them
we observe an increase in sensitivity of $100\%$–$130\%$.

Figure 4:  Current (left) and future (right) lower limits of $1/\sqrt{C}$ at
$90\%$ CL for the scenarios with an additional massive fourth neutrino
interacting solely via one of the given NP Wilson coefficients as a function
of its mass $m_{\nu,4}$.

Finally, in the context of single-operator fits we also study scenario (iii).
Namely, the case of having a fourth massive Majorana neutrino, $\nu_{4}$, with
a mass such that the decays $K^{+}\to\pi^{+}\nu\nu_{4}$ or
$K^{+}\to\pi^{+}\nu_{4}\nu_{4}$ are kinematically possible. To illustrate the
possible impact of a fourth massive neutrino on the distributions we choose
the benchmark mass $m_{\nu,4}=50$ MeV and show in Figure 3 the analogous
distributions as for the three massless neutrino cases.

Figure 4 shows the current constraints at $90\%$ CL as a function of the
fourth neutrino mass $m_{\nu,4}$. The left panel shows the current sensitivity
while the right one the projected sensitivity. All limits are given at $90\%$
CL. Here we observe an increase of $50\%$–$60\%$ in the sensitivities
depending on the NP operator. Note that the kinematic endpoint of the
contribution of the diagonal operators, i.e., $O^{S/P,\,L}_{4412/21}$ and
$O^{A,\,L/R}_{4412}$ (dashed lines), lies at around $\sim 130\,\text{MeV}$ as
the final state necessarily contains two massive neutrinos.

Massless Dirac Neutrinos

Figure 5: Allowed parameter space at $68\%$ CL (dark blue) and $90\%$ CL
(light blue) in the plane of two NP Wilson coefficients corresponding to
various combinations of scenario (i), i.e. the case of three massless Dirac
neutrinos. Note that here $a\neq b$. The dashed line corresponds to the future
projection for the HIKE experiment at $90\%$ CL.

Massless Majorana Neutrinos

Figure 6: Allowed parameter space at $68\%$ CL (dark blue) and $90\%$ CL
(light blue) in the plane of two NP Wilson coefficients corresponding to
various combinations of scenario (ii), i.e. the case of three massless
Majorana neutrinos. Note that here $a\neq b$. The dashed line corresponds to
the future projection for the HIKE experiment at $90\%$ CL.

###  4.2 Multi-operator correlations

Next we study the correlations of probing different operators by considering
the effect of switching on pairs of different operators simultaneously. For
concreteness we cover (almost) all different non-trivial combinations for the
case of three massless Dirac or Majorana neutrinos, scenario (i) and (ii),
respectively. For these correlations it can be important to specify the phase
we chose for the NP Wilson coefficients, since relative phases can be relevant
for interferences. We assume there is no new CP-violating phases beyond the
SM, i.e., all NP phases are aligned to the SM one. Additionally, in the
discussion of neutrino-flavour diagonal operators, e.g., $O^{V,\,LL}_{aa12}$,
we consider the flavour-universal case in which we switch on all operators of
the category simultaneously but with universal Wilson coefficients, e.g.,
$C^{V,\,LL}_{1112}=C^{V,\,LL}_{2212}=C^{V,\,LL}_{3312}$. This additional
assumption of flavour universality reduces the number of different
combinations.

The results are shown in Figure 5 and 6. For each case we determine the
allowed region at $68\%$ CL (dark blue region) and $90\,\%$ CL (light blue
region) in the plane of two Wilson coefficients assuming that the NP scenario
is realised in nature. The dashed, black lines correspond to the $90\%$ CL
future projection.

##  5 Discussion and conclusions

The golden, rare kaon decays $K\to\pi\nu\nu$ provide an excellent test of the
Standard Model and are already placing stringent constraints on the parameter
space of New Physics. The expected and planned experimental activities will
test modifications of the SM operators at energies of
$\mathcal{O}(300\,\text{TeV})$ and will also allow to extract additional
information from the invariant mass spectrum of the neutrinos. In the weak
effective theory, scalar and tensor operators induce a modified invisible mass
distribution, c.f. Figure 1. We we have classified all independent Majorana
and Dirac operators that contribute to these FCNC decays to neutrinos. We have
derived the sensitivity to each operator for the current and expected future
experimental situation including experimental and theoretical uncertainties,
see Tables 3 and 4. The expected constraints on the NP scale range from
$\mathcal{O}(40\,\text{TeV})$ for tensor operators to
$\mathcal{O}(300\,\text{TeV})$ for axial-vector Majorana and vector Dirac
operators. We have studied the impact of a massive sterile neutrino on the
invisible mass spectrum and have discussed the correlation of various
operators by analysing the allowed parameter space for pairs of different NP
operators.

Majorana and Dirac neutrinos, which can be generated from a lepton-number
violating SMEFT or a lepton-number conserving $\nu$SMEFT theory with
additional right-handed neutrinos, have a completely different UV origin.
Therefore, the relevant weak effective theories also contain different number
of independent parameters. This implies that concrete UV models could be
tested from the data of the spectrum in future searches. The NP sensitivity of
the invisible spectrum provides a strong motivation for the future Kaon
program and we plan to study the discriminating power of the charged decay
mode spectrum on different realisations of the SMEFT and $\nu$SMEFT in the
future [39].

#### Acknowledgments

We thank Radoslav Marchevski for the discussion regarding the future prospects
of the HIKE experiment. Moldanazarova Ulserik was supported by the grant
AP14972893 from the Committee of Science of the Ministry of Science and Higher
Education of the Republic of Kazakhstan. The work of MG is partially supported
by the UK Science and Technology Facilities Council grant ST/X000699/1. UM and
MG thanks the particle theory group at the TU Dortmund for the hospitality
during the time in which this work was in part done.

##  Appendix A Decay rates

The full expressions for the double differential decay widths
$\mathrm{d}^{2}\Gamma_{ab}/\mathrm{d}q^{2}\mathrm{d}|\vec{p}_{\pi}|$ and
$\mathrm{d}^{2}\Gamma_{ab}/\mathrm{d}q^{2}\mathrm{d}k^{2}$ for massive
Majorana and Dirac neutrinos are openly available in a git repository on
GitLab [40]. Here we give expressions for the double differential partial
decay width in terms of the Lorentz invariant Dalitz variables
$q^{2}=(p_{K}-p_{\pi})^{2}$ and $k^{2}=(p_{\pi}+p_{\nu,a})^{2}$ for massless
neutrinos. In the Majorana case we find

$\displaystyle\frac{\mathrm{d}^{2}\Gamma_{ab}}{\mathrm{d}q^{2}\mathrm{d}k^{2}}\Bigg{|}_{\text{Majorana}}=$
$\displaystyle\frac{1}{512\pi^{3}m_{K}^{3}(1+\delta_{ab})}\Theta_{\text{L.
inv.}}\left(q^{2},k^{2}\right)\Bigl{|}_{m_{a}=m_{b}=0}\Bigg{[}$ (28)
$\displaystyle-4\left|f^{V}(q^{2})\right|^{2}(k^{4}-k^{2}(m_{K}^{2}+m_{\pi}^{2}-q^{2})+m_{K}^{2}m_{\pi}^{2})$
$\displaystyle\qquad\times\left(\left|C^{A,\,L,\,\text{SM}}_{ab12}+C^{A,\,L}_{ab12}+C^{A,\,R}_{ab12}\right|^{2}+\left|C^{V,\,L}_{ab12}+C^{V,\,R}_{ab12}\right|^{2}\right)$
$\displaystyle+q^{2}\left|F^{S}(q^{2})\right|^{2}\Big{(}\left|C^{P,\,L}_{ab12}+C^{P,\,L*}_{ab21}\right|^{2}+\left|C^{S,\,L}_{ab12}+C^{S,\,L*}_{ab21}\right|^{2}\Big{)}$
$\displaystyle+q^{2}\left|F^{T}(q^{2})\right|^{2}(-2k^{2}+m_{K}^{2}+m_{\pi}^{2}-q^{2})^{2}\Big{(}\left|C^{T,\,L}_{ab12}-C^{T,\,L*}_{ab21}\right|^{2}\Big{)}$
$\displaystyle+2q^{2}(2k^{2}-m_{K}^{2}-m_{\pi}^{2}+q^{2})\Big{(}$
$\displaystyle\qquad-F^{S}(q^{2})F^{T}(q^{2})\mathrm{Im}\Big{(}C^{S,\,L}_{ab12}C^{T,\,L}_{ab21}+C^{S,\,L}_{ab21}C^{T,\,L}_{ab12}\Big{)}$
$\displaystyle\qquad+F^{S}(q^{2})\left(F^{T}(q^{2})\right)^{*}\mathrm{Im}\Big{(}C^{S,\,L}_{ab21}C^{T,\,L*}_{ab21}+C^{S,\,L}_{ab12}C^{T,\,L*}_{ab12}\Big{)}\Big{)}\Bigg{]}\,,$

with the theta function in terms of $q^{2}$ and $k^{2}$ defined as

$\begin{split}\Theta_{\text{L.
inv.}}\left(q^{2},k^{2}\right)\Bigl{|}_{m_{a}=m_{b}=0}&=\theta(q^{2})\theta\bigl{(}(m_{K}-m_{\pi})^{2}-q^{2}\bigr{)}\\\
&\times\theta\left(k^{2}-\frac{1}{2}\left(m_{K}^{2}+m_{\pi}^{2}-q^{2}-\lambda^{1/2}_{K\pi
q^{2}}\right)\right)\\\
&\times\theta\left(\frac{1}{2}\left(m_{K}^{2}+m_{\pi}^{2}-q^{2}+\lambda^{1/2}_{K\pi
q^{2}}\right)-k^{2}\right)\,.\end{split}$ (29)

For massless Dirac neutrinos we similarly find

$\displaystyle\frac{\mathrm{d}^{2}\Gamma_{ab}}{\mathrm{d}q^{2}\mathrm{d}k^{2}}$
$\displaystyle=\frac{1}{256\pi^{3}m_{K}^{3}}\Theta_{\text{L.
inv.}}\left(q^{2},k^{2}\right)\Bigl{|}_{m_{a}=m_{b}=0}$ (30)
$\displaystyle\times\Bigg{(}-(f^{V}(q^{2}))^{2}(k^{4}-k^{2}(m_{K}^{2}+m_{\pi}^{2}-q^{2})+m_{K}^{2}m_{\pi}^{2})$
$\displaystyle\qquad\times\left(\left|C^{V,\,LL,\,\text{SM}}_{ab12}+C^{V,\,LL}_{ab12}+C^{V,\,LR}_{ab12}\right|^{2}+\left|C^{V,\,RL}_{ab12}+C^{V,\,RR}_{ab12}\right|^{2}\right)$
$\displaystyle+\frac{q^{2}}{4}\left((F^{S}(q^{2}))^{2}\left(\left|C^{S,\,LL}_{ab12}+C^{S,\,LR}_{ab12}\right|^{2}+\left|C^{S,\,LR}_{ba21}+C^{S,\,LL}_{ba21}\right|^{2}\right)\right.$
$\displaystyle\qquad\left.+(F^{T}(q^{2}))^{2}(-2k^{2}+m_{K}^{2}+m_{\pi}^{2}-q^{2})^{2}\left(\left|C^{T,\,LL}_{ab12}\right|^{2}+\left|C^{T,\,LL}_{ba21}\right|^{2}\right)\right)$
$\displaystyle-\frac{q^{2}}{2}(-2k^{2}+m_{K}^{2}+m_{\pi}^{2}-q^{2})\Big{(}$
$\displaystyle\qquad+F^{S}(q^{2})F^{T}(q^{2})\mathrm{Im}\left(C^{T,\,LL}_{ab12}\left(C^{S,\,LL}_{ba21}+C^{S,\,LR}_{ba21}\right)\right)$
$\displaystyle\qquad+F^{S}(q^{2})(F^{T}(q^{2}))^{*}\mathrm{Im}\left(C^{T,\,LL*}_{ba21}\left(C^{S,\,LL}_{ba21}+C^{S,\,LR}_{ba21}\right)\right)$
$\displaystyle\qquad+(F^{S}(q^{2}))^{*}F^{T}(q^{2})\mathrm{Im}\left(C^{T,\,LL}_{ab12}\left(C^{S,\,LL}_{ab12}+C^{S,\,LR}_{ab12}\right)^{*}\right)$
$\displaystyle\qquad+(F^{S}(q^{2}))^{*}(F^{T}(q^{2}))^{*}\mathrm{Im}\left(C^{T,\,LL*}_{ba21}\left(C^{S,\,LL}_{ab12}+C^{S,\,LR}_{ab12}\right)^{*}\right)\Big{)}\Bigg{)}\,.$

Note, that for massless neutrinos the interference term between the scalar and
tensor operators vanishes for both Majorana, and Dirac neutrinos, after
integration over the full $k^{2}$ phase space, which is why such terms do not
appear in Eqs. (21) and (23).

## References

  * [1] E. Cortina Gil, A. Kleimenova, E. Minucci, S. Padolski, P. Petrov, A. Shaikhiev et al., _An investigation of the very rare ${K}^{+}\to{\pi}^{+}\nu\overline{\nu}$ decay_, _Journal of High Energy Physics_ 2020 (Nov., 2020) .
  * [2] NA62 collaboration, E. Cortina Gil et al., _Measurement of the very rare K +→${\pi}^{+}\nu\overline{\nu}$ decay_, _JHEP_ 06 (2021) 093, [2103.15389].
  * [3] KOTO collaboration, J. K. Ahn et al., _Search for the $K_{L}\\!\to\\!\pi^{0}\nu\overline{\nu}$ and $K_{L}\\!\to\\!\pi^{0}X^{0}$ decays at the J-PARC KOTO experiment_, _Phys. Rev. Lett._ 122 (2019) 021802, [1810.09655].
  * [4] KOTO collaboration, J. K. Ahn et al., _Study of the $K_{L}\to\pi^{0}\nu\bar{\nu}$ Decay at the J-PARC KOTO Experiment_, _Phys. Rev. Lett._ 126 (2021) 121801, [2012.07571].
  * [5] E. Goudzovski et al., _New physics searches at kaon and hyperon factories_ , _Rept. Prog. Phys._ 86 (2023) 016201, [2201.07805].
  * [6] G. Anzivino et al., _Workshop summary – Kaons@CERN 2023_ , in _Kaons@CERN 2023_ , 11, 2023, 2311.02923.
  * [7] S. L. Glashow, J. Iliopoulos and L. Maiani, _Weak Interactions with Lepton-Hadron Symmetry_ , _Phys. Rev. D_ 2 (1970) 1285–1292.
  * [8] G. Buchalla and A. J. Buras, _The rare decays $K\to\pi\nu\bar{\nu}$, $B\to X\nu\bar{\nu}$ and $B\to l^{+}l^{-}$: An Update_, _Nucl. Phys. B_ 548 (1999) 309–327, [hep-ph/9901288].
  * [9] M. Misiak and J. Urban, _QCD corrections to FCNC decays mediated by Z penguins and W boxes_ , _Phys. Lett. B_ 451 (1999) 161–169, [hep-ph/9901278].
  * [10] A. J. Buras, M. Gorbahn, U. Haisch and U. Nierste, _Charm quark contribution to $K^{+}\to\pi^{+}\nu\bar{\nu}$ at next-to-next-to-leading order_, _JHEP_ 11 (2006) 002, [hep-ph/0603079].
  * [11] F. Mescia and C. Smith, _Improved estimates of rare K decay matrix-elements from Kl3 decays_ , _Phys. Rev. D_ 76 (2007) 034017, [0705.2025].
  * [12] J. Brod and M. Gorbahn, _Electroweak Corrections to the Charm Quark Contribution to $K^{+}\to\pi^{+}\nu\bar{\nu}$_, _Phys. Rev. D_ 78 (2008) 034006, [0805.4119].
  * [13] J. Brod, M. Gorbahn and E. Stamou, _Two-Loop Electroweak Corrections for the $K\to\pi\nu\bar{\nu}$ Decays_, _Phys. Rev. D_ 83 (2011) 034030, [1009.0947].
  * [14] G. Isidori, F. Mescia and C. Smith, _Light-quark loops in K — $>$ pi nu anti-nu_, _Nucl. Phys. B_ 718 (2005) 319–338, [hep-ph/0503107].
  * [15] Z. Bai, N. H. Christ, X. Feng, A. Lawson, A. Portelli and C. T. Sachrajda, _$K^{+}\to\pi^{+}\nu\bar{\nu}$ decay amplitude from lattice QCD_, _Phys. Rev. D_ 98 (2018) 074509, [1806.11520].
  * [16] RBC, UKQCD collaboration, N. H. Christ, X. Feng, A. Portelli and C. T. Sachrajda, _Lattice QCD study of the rare kaon decay $K^{+}\to\pi^{+}\nu\bar{\nu}$ at a near-physical pion mass_, _Phys. Rev. D_ 100 (2019) 114506, [1910.10644].
  * [17] J. Brod, M. Gorbahn and E. Stamou, _Updated standard model prediction for $k\to\pi\nu\bar{\nu}$ and $\epsilon_{k}$_, 2021.
  * [18] Particle Data Group collaboration, R. L. Workman, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * [19] HIKE collaboration, M. U. Ashraf et al., _High Intensity Kaon Experiments (HIKE) at the CERN SPS Proposal for Phases 1 and 2_ , 2311.08231.
  * [20] KOTO collaboration, H. Nanjo, _KOTO II at J-PARC : toward measurement of the branching ratio of $K_{L}\to\pi^{0}\nu\bar{\nu}$_, _J. Phys. Conf. Ser._ 2446 (2023) 012037.
  * [21] D. Aristizabal Sierra, V. De Romeri and N. Rojas, _COHERENT analysis of neutrino generalized interactions_ , _Phys. Rev. D_ 98 (2018) 075018, [1806.07424].
  * [22] W. Altmannshofer, M. Tammaro and J. Zupan, _Non-standard neutrino interactions and low energy experiments_ , _JHEP_ 09 (2019) 083, [1812.02778].
  * [23] A. Falkowski, M. González-Alonso and Z. Tabrizi, _Reactor neutrino oscillations as constraints on Effective Field Theory_ , _JHEP_ 05 (2019) 173, [1901.04553].
  * [24] I. Bischer and W. Rodejohann, _General neutrino interactions from an effective field theory perspective_ , _Nucl. Phys. B_ 947 (2019) 114746, [1905.08699].
  * [25] T. Han, J. Liao, H. Liu and D. Marfatia, _Scalar and tensor neutrino interactions_ , _JHEP_ 07 (2020) 207, [2004.13869].
  * [26] A. Abada, D. Bečirević, O. Sumensari, C. Weiland and R. Zukanovich Funchal, _Sterile neutrinos facing kaon physics experiments_ , _Phys. Rev. D_ 95 (2017) 075023, [1612.04737].
  * [27] E. E. Jenkins, A. V. Manohar and P. Stoffer, _Low-Energy Effective Field Theory below the Electroweak Scale: Operators and Matching_ , _JHEP_ 03 (2018) 016, [1709.04486].
  * [28] F. F. Deppisch, K. Fridell and J. Harz, _Constraining lepton number violating interactions in rare kaon decays_ , _JHEP_ 12 (2020) 186, [2009.04494].
  * [29] G. Buchalla, A. J. Buras and M. E. Lautenbacher, _Weak decays beyond leading logarithms_ , _Rev. Mod. Phys._ 68 (1996) 1125–1144, [hep-ph/9512380].
  * [30] H. K. Dreiner, H. E. Haber and S. P. Martin, _Two-component spinor techniques and Feynman rules for quantum field theory and supersymmetry_ , _Phys. Rept._ 494 (2010) 1–196, [0812.1594].
  * [31] A. Denner, H. Eck, O. Hahn and J. Kublbeck, _Compact Feynman rules for Majorana fermions_ , _Phys. Lett. B_ 291 (1992) 278–280.
  * [32] R.-X. Shi et al., _Revisiting the new-physics interpretation of the b $\rightarrow$ c$\tau$$\nu$ data_, _Journal of High Energy Physics_ 2019 (12, 2019) .
  * [33] I. Baum et al., _Matrix elements of the electromagnetic operator between kaon and pion states_ , _Physical Review D_ 84 (10, 2011) .
  * [34] F. Brizioli, _Measurement of $Br(K^{+}\to\pi^{+}\nu\bar{\nu})$ with the NA62 experiment at CERN_, Ph.D. thesis, Perugia U., 2021.
  * [35] E. Cortina Gil, E. Minucci, S. Padolski, P. Petrov, B. Velghe, G. Georgiev et al., _First search for k +→${\pi}^{+}\nu\overline{\nu}$ using the decay-in-flight technique_, _Physics Letters B_ 791 (Apr., 2019) 156–166.
  * [36] W. A. Rolke, A. M. Lopez and J. Conrad, _Limits and confidence intervals in the presence of nuisance parameters_ , _Nucl. Instrum. Meth. A_ 551 (2005) 493–503, [physics/0403059].
  * [37] K. Cranmer, _Practical Statistics for the LHC_ , in _2011 European School of High-Energy Physics_ , pp. 267–308, 2014, 1503.07622, DOI.
  * [38] HIKE collaboration, E. Cortina Gil et al., _HIKE, High Intensity Kaon Experiments at the CERN SPS: Letter of Intent_ , 2211.16586.
  * [39] M. Gorbahn, U. Moldanazarova, K. H. Sieja, E. Stamou and M. Tabet, _in preparation_ , .
  * [40] https://gitlab.com/manstam/kpinunu_distributions_in_left, 2023.