aainstitutetext: AHEP Group, Instituto de Física Corpuscular –
CSIC/Universitat de València, Apartado 22085, E–46071 València,
Spainbbinstitutetext: Departamento de Ciencias, Facultad de Artes Liberales,
Universidad Adolfo Ibáñez, Diagonal Las Torres 2640, Santiago,
Chileccinstitutetext: Millennium Institute for Subatomic Physics at the High
Energy Frontier (SAPHIR), Fernández Concha 700, Santiago,
Chileddinstitutetext: Departamento de Física, Facultad de Ciencias,
Universidad de La Serena, Avenida Cisternas 1200, La Serena, Chile
eeinstitutetext: Departament de Física Teòrica, Universitat de València and
Instituto de Física Corpuscular – CSIC/Universitat de València, Dr. Moliner
50, E–46100 Burjassot, Spainffinstitutetext: Department of Physics, National
Tsing Hua University, Hsinchu 300, Taiwangginstitutetext: Center for Theory
and Computation, National Tsing Hua University, Hsinchu 300, Taiwan

# Long-lived heavy neutral leptons from mesons in effective field theory

Rebeca Beltrán<EMAIL_ADDRESS>b,c Giovanna Cottin
<EMAIL_ADDRESS>d,c Juan Carlos Helo<EMAIL_ADDRESS>a Martin Hirsch
<EMAIL_ADDRESS>e Arsenii Titov<EMAIL_ADDRESS>f,g Zeren Simon
Wang<EMAIL_ADDRESS>
###### Abstract

In the framework of the low-energy effective field theory of the Standard
Model extended with heavy neutral leptons (HNLs), we calculate the production
rates of HNLs from meson decays triggered by dimension-six operators. We
consider both lepton-number-conserving and lepton-number-violating four-
fermion operators involving either a pair of HNLs or a single HNL. Assuming
that HNLs are long-lived, we perform simulations and investigate the reach of
the proposed far detectors at the high-luminosity LHC to (i) active-heavy
neutrino mixing and (ii) the Wilson coefficients associated with the effective
operators, for HNL masses below the mass of the $B$-meson. We further convert
the latter to the associated new-physics scales. Our results show that scales
in excess of hundreds of TeV and the active-heavy mixing squared as small as
$10^{-15}$ can be probed by these experiments.

††preprint:  FTUV-22-1005.8958 IFIC/22-27

##  1 Introduction

The next 10–15 years will see a notable increase in experimental sensitivity
to the search for heavy neutral leptons (HNLs), mostly, but not exclusively,
thanks to the diverse long-lived particle program that has been initiated at
CERN for the future high-luminosity LHC (HL-LHC) Alimena:2019zri ;
Curtin:2018mvb . Many new and dedicated “far-detector” proposals have been put
forward in this context: ANUBIS Bauer:2019vqk , AL3X Gligorov:2018vkc ,
CODEX-b Gligorov:2017nwh , FACET Cerci:2021nlb , FASER and FASER2 Feng:2017uoz
; FASER:2018eoc , MoEDAL-MAPP1 and MAPP2 Pinfold:2019nqj ; Pinfold:2019zwp ,
and MATHUSLA Chou:2016lxi ; Curtin:2018mvb ; MATHUSLA:2020uve . Some of them
have already been approved and even deployed. In particular, FASER and MoEDAL-
MAPP1 will be taking data during recently started Run 3 of the LHC, whereas
the rest of the program is envisioned for the HL phase.

Consequently, a number of papers have studied the discovery prospects of these
(and other) future experiments for minimal HNLs, i.e. HNLs coupled to the
Standard Model (SM) only through active-heavy mixing, see Refs. Kling:2018wct
; Helo:2018qej ; Dercks:2018wum ; Hirsch:2020klk ; DeVries:2020jbs ;
Ovchynnikov:2022its . Already in the minimal scenario one can expect up to
four orders of magnitude improvement in mixing angles squared relative to
existing constraints, for HNLs with masses up to roughly 5 GeV. Similarly,
also the main detectors at the LHC can search for HNLs. However, ATLAS/CMS
will be more sensitive at larger HNL masses, say, roughly 3–30 GeV; see e.g.
Refs. Cottin:2018kmq ; Cottin:2018nms ; Cheung:2020buy ; Jones-Perez:2019plk ;
Liu:2019ayx . Note that ATLAS Aad:2019kiz  as well as CMS CMS:2022fut  have
already published first results of searches for HNLs using displaced vertices.

HNLs could appear also in SM extensions with additional particles, a prominent
example being a $Z^{\prime}$. For sensitivity estimates for this scenario, see
e.g. Refs. Chiang:2019ajm ; Deppisch:2019kvs . However, since the LHC has so
far not found any unambiguous signal of beyond the Standard Model (BSM)
physics, attention has shifted recently towards effective field theories
(EFTs). At the center-of-mass energy of the LHC, one expects that the correct
description is the Standard Model effective field theory (SMEFT) (see Ref.
Brivio:2017vri  for a review) or, if we allow also for the existence of light
(with masses below or around the electroweak scale) right-handed (RH) gauge
singlet fermions $N_{R}$, the $N_{R}$SMEFT delAguila:2008ir ; Aparici:2009fh ;
Liao:2016qyd ; Li:2021tsq . Recently, we have studied the expectation for HNL
searches in the $N_{R}$SMEFT for pair-$N_{R}$ operators Cottin:2021lzz  and
for single-$N_{R}$ operators Beltran:2021hpq .

However, below the electroweak scale, and in particular, at the scale of meson
masses, the massive gauge and Higgs bosons, as well as the top quark, have to
be integrated out. The correct EFT in this lower energy regime is then usually
called LEFT (low-energy effective field theory) Jenkins:2017jig , or in case
we add RH singlet states, $N_{R}$LEFT Bischer:2019ttk ; Chala:2020vqp ;
Li:2020lba ; Li:2020wxi . The sensitivities of future far detectors at the LHC
for the four-fermion $N_{R}$LEFT operators containing a single $N_{R}$, a
charged lepton $\ell$, and a quark bilinear, and thus, triggering decays of
mesons (produced in $pp$ collisions at the LHC) into $N_{R}$ and $\ell$, have
already been studied in Ref. DeVries:2020jbs . Furthermore, the projected
sensitivities of Belle II to the same set of operators have been derived in
Ref. Zhou:2021ylt  (see also Ref. Han:2022uho ). Therefore, in the current
work, we will concentrate on pair-$N_{R}$ operators and operators with one
$N_{R}$ and one SM neutrino, exploring the potential of future far detectors
at the LHC. We note that Ref. Li:2020lba  has studied a similar set of
operators and derived bounds from current data on the coherent elastic
neutrino-nucleus scattering (CE${\nu}$NS) process and invisible decays of
light unflavored mesons.

The rest of this paper is organized as follows. In the next section, we
discuss $N_{R}$LEFT operators at $d=6$ and their matching to the $N_{R}$SMEFT
at $d=6$ and $d=7$. Then in section 3, we elaborate on HNL production from
meson decays in the framework of the $N_{R}$LEFT. Section 4 gives a brief
summary of the experiments we consider and provides the corresponding
simulation details. In section 5, we present the results of our simulation. We
then close with a brief summary. Technical details of the computations of
meson and HNL decay widths in the $N_{R}$LEFT are relegated to the appendices.

##  2 Effective field theory approach

###  2.1 $N_{R}$LEFT: low-energy effective field theory with right-handed
neutrinos

In the present work, we are interested in the HNL production from meson decays
in an EFT approach. Meson decays take place at energy scales
$E\sim\mathcal{O}(1)$ GeV, much smaller than the electroweak scale $v=246$
GeV. At such energies, the $W$ and $Z$ gauge bosons, the Higgs boson and the
top quark are not propagating degrees of freedom, and the effects of potential
new physics (NP) associated with a scale $\Lambda\gg E$ can be described in
terms of the LEFT extended with RH neutrinos (or HNLs), $N_{R}$. We will refer
to such minimal extension of the LEFT as $N_{R}$LEFT. This EFT respects
$SU(3)_{C}\times U(1)_{\mathrm{em}}$ gauge symmetry and in addition to the
gluons and the photon contains all SM leptons, $u$, $d$, $s$, $c$, $b$ quarks
and, in general, an arbitrary number $n_{N}$ of HNLs.

In addition to the renormalizable terms and higher-dimensional operators
constructed from the SM fields that have been systematically classified up to
$d=6$ in Ref. Jenkins:2017jig , the $N_{R}$LEFT Lagrangian includes operators
with $N_{R}$, which at $d\leq 6$ have been constructed in Refs. Chala:2020vqp
; Li:2020lba . The renormalizable part of the Lagrangian is given by

$\mathcal{L}_{\mathrm{ren}}=\mathcal{L}_{\mathrm{QCD+QED}}+\overline{N_{R}}i\not{\partial}N_{R}-\left[\frac{1}{2}\overline{\nu_{L}}M_{\nu}\nu_{L}^{c}+\frac{1}{2}\overline{N_{R}^{c}}M_{N}N_{R}+\overline{\nu_{L}}M_{D}N_{R}+\text{h.c.}\right],$
(1)

where $\mathcal{L}_{\mathrm{QCD+QED}}$ is the QCD and QED Lagrangian built out
of the SM leptons (including left-handed (LH) neutrinos $\nu_{L}$ without any
mass term) and the five aforementioned flavors of quarks (for its form, see
e.g. Eq. (5.1) in Ref. Jenkins:2017jig ). The superscript $c$ denotes a charge
conjugate field, $\psi^{c}\equiv C\overline{\psi}^{T}$. The Majorana mass
matrices $M_{\nu}$ and $M_{N}$ are complex symmetric $3\times 3$ and
$n_{N}\times n_{N}$ matrices, respectively. $M_{D}$ is a generic $3\times
n_{N}$ Dirac mass matrix. If lepton number is conserved, the Majorana mass
terms are not present. The full $N_{R}$LEFT Lagrangian reads

$\mathcal{L}_{N_{R}\mathrm{LEFT}}=\mathcal{L}_{\mathrm{ren}}+\sum_{d\geq
5}\sum_{i}c_{i}^{(d)}\mathcal{O}_{i}^{(d)}\,,$ (2)

where $c_{i}^{(d)}$ are the Wilson coefficients of higher-dimensional
operators $\mathcal{O}_{i}^{(d)}$, and the second sum goes over all
independent operators of a given mass dimension $d$. The dimensionful
coefficients $c_{i}^{(d)}$ scale as $\Lambda^{4-d}$.

At $d=5$, there are two dipole operators with $N_{R}$:

$\mathcal{O}_{NN\gamma}=\overline{N_{R}^{c}}\sigma^{\mu\nu}N_{R}F_{\mu\nu}\qquad\text{and}\qquad\mathcal{O}_{\nu
N\gamma}=\overline{\nu_{L}}\sigma^{\mu\nu}N_{R}F_{\mu\nu}\,,$ (3)

where $F_{\mu\nu}$ is the photon field strength tensor. The first of these
operators is antisymmetric in the HNL generation indices, and thus, it
vanishes identically for $n_{N}=1$. $\mathcal{O}_{NN\gamma}$ is lepton-number-
violating (LNV), whereas $\mathcal{O}_{\nu N\gamma}$ is lepton-number-
conserving (LNC). The phenomenology of the RH neutrino dipole operator has
been studied in Ref. Aparici:2009fh  and more recently in Ref.
Barducci:2022gdv  in the context of the proposed far detectors at the LHC.

At $d=6$, there are 55 four-fermion interactions, of which 23 LNC
Chala:2020vqp , 26 LNV, and six operators that violate both lepton and baryon
numbers Li:2020lba .111Here, we do not count possible flavor structures and
Hermitian conjugates. The contact operators triggering meson decays into HNLs
involve a pair of quarks and a leptonic part, which in turn contains either
(i) a pair of HNLs, or (ii) $N_{R}$ and $\nu_{L}$, or else (iii) $N_{R}$ and a
charged lepton. The last possibility has been investigated in detail in Refs.
DeVries:2020jbs  and Zhou:2021ylt , in the context of long-lived HNLs that
could be produced in meson decays at the high-luminosity LHC and Belle II,
respectively. In this work, we focus on the operators with neutral leptons. In
tables 1 and 2, we summarize the pair-$N_{R}$ and single-$N_{R}$ operators of
interest.

| Name | Structure | $n_{N}=1$ | $n_{N}=3$
---|---|---|---|---
LNC | ${\cal O}_{dN}^{V,RR}$ | $\left(\overline{d_{R}}\gamma_{\mu}d_{R}\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 9 | 81
${\cal O}_{uN}^{V,RR}$ | $\left(\overline{u_{R}}\gamma_{\mu}u_{R}\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 4 | 36
${\cal O}_{dN}^{V,LR}$ | $\left(\overline{d_{L}}\gamma_{\mu}d_{L}\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 9 | 81
${\cal O}_{uN}^{V,LR}$ | $\left(\overline{u_{L}}\gamma_{\mu}u_{L}\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 4 | 36
LNV | ${\cal O}_{dN}^{S,RR}$ | $\left(\overline{d_{L}}d_{R}\right)\left(\overline{N_{R}^{c}}N_{R}\right)$ | 18 | 108
${\cal O}_{dN}^{T,RR}$ | $\left(\overline{d_{L}}\sigma_{\mu\nu}d_{R}\right)\left(\overline{N_{R}^{c}}\sigma^{\mu\nu}N_{R}\right)$ | 0 | 54
${\cal O}_{uN}^{S,RR}$ | $\left(\overline{u_{L}}u_{R}\right)\left(\overline{N_{R}^{c}}N_{R}\right)$ | 8 | 48
${\cal O}_{uN}^{T,RR}$ | $\left(\overline{u_{L}}\sigma_{\mu\nu}u_{R}\right)\left(\overline{N_{R}^{c}}\sigma^{\mu\nu}N_{R}\right)$ | 0 | 24
${\cal O}_{dN}^{S,LR}$ | $\left(\overline{d_{R}}d_{L}\right)\left(\overline{N_{R}^{c}}N_{R}\right)$ | 18 | 108
${\cal O}_{uN}^{S,LR}$ | $\left(\overline{u_{R}}u_{L}\right)\left(\overline{N_{R}^{c}}N_{R}\right)$ | 8 | 48
Table 1: Dimension-six operators in the $N_{R}$LEFT, involving two quarks and two HNLs. For each operator structure, we give the number of independent real parameters for $n_{N}=1$ and $n_{N}=3$ generations of $N_{R}$. We recall that in the $N_{R}$LEFT, $n_{d}=3$ and $n_{u}=2$. The LNV operator structures require “+h.c.”. | Name | Structure | $n_{N}=1$ | $n_{N}=3$
---|---|---|---|---
LNC | ${\cal O}_{d\nu N}^{S,RR}$ | $\left(\overline{d_{L}}d_{R}\right)\left(\overline{\nu_{L}}N_{R}\right)$ | 54 | 162
${\cal O}_{d\nu N}^{T,RR}$ | $\left(\overline{d_{L}}\sigma_{\mu\nu}d_{R}\right)\left(\overline{\nu_{L}}\sigma^{\mu\nu}N_{R}\right)$ | 54 | 162
${\cal O}_{u\nu N}^{S,RR}$ | $\left(\overline{u_{L}}u_{R}\right)\left(\overline{\nu_{L}}N_{R}\right)$ | 24 | 72
${\cal O}_{u\nu N}^{T,RR}$ | $\left(\overline{u_{L}}\sigma_{\mu\nu}u_{R}\right)\left(\overline{\nu_{L}}\sigma^{\mu\nu}N_{R}\right)$ | 24 | 72
${\cal O}_{d\nu N}^{S,LR}$ | $\left(\overline{d_{R}}d_{L}\right)\left(\overline{\nu_{L}}N_{R}\right)$ | 54 | 162
${\cal O}_{u\nu N}^{S,LR}$ | $\left(\overline{u_{R}}u_{L}\right)\left(\overline{\nu_{L}}N_{R}\right)$ | 24 | 72
LNV | ${\cal O}_{d\nu N}^{V,RR}$ | $\left(\overline{d_{R}}\gamma_{\mu}d_{R}\right)\left(\overline{\nu_{L}^{c}}\gamma^{\mu}N_{R}\right)$ | 54 | 162
${\cal O}_{u\nu N}^{V,RR}$ | $\left(\overline{u_{R}}\gamma_{\mu}u_{R}\right)\left(\overline{\nu_{L}^{c}}\gamma^{\mu}N_{R}\right)$ | 24 | 72
${\cal O}_{d\nu N}^{V,LR}$ | $\left(\overline{d_{L}}\gamma_{\mu}d_{L}\right)\left(\overline{\nu_{L}^{c}}\gamma^{\mu}N_{R}\right)$ | 54 | 162
${\cal O}_{u\nu N}^{V,LR}$ | $\left(\overline{u_{L}}\gamma_{\mu}u_{L}\right)\left(\overline{\nu_{L}^{c}}\gamma^{\mu}N_{R}\right)$ | 24 | 72
Table 2: Dimension-six operators in the $N_{R}$LEFT, involving two quarks, one
active neutrino and one HNL. For each operator structure, we give the number
of independent real parameters for $n_{N}=1$ and $n_{N}=3$ generations of
$N_{R}$. We recall that in the $N_{R}$LEFT, $n_{\nu}=n_{d}=3$ and $n_{u}=2$.
All operator structures require “+h.c.”.

Though in our numerical analysis we will focus on the case of a single HNL
generation ($n_{N}=1$), in the tables we provide for completeness the numbers
of independent real parameters associated with each operator for $n_{N}=1$ and
$n_{N}=3$. We have checked these numbers using the Sym2Int package
Fonseca:2017lem ; Fonseca:2019yya . In the group of pair-$N_{R}$ operators,
there are four LNC structures and six LNV ones (of which two vanish
identically for $n_{N}=1$). Among the single-$N_{R}$ operators, six structures
conserve lepton number, whereas four violate it.

###  2.2 Matching to the $N_{R}$SMEFT

In this subsection, we assume that the $N_{R}$LEFT originates from the
$N_{R}$SMEFT,222It is worth noting that ($N_{R}$)LEFT is the correct effective
description at low energies even if the high-energy EFT is not given by the
($N_{R}$)SMEFT but by the HEFT (with $N_{R}$), in which the Higgs is not part
of a fundamental weak doublet Jenkins:2017jig . i.e. the EFT of the SM
extended with RH neutrinos delAguila:2008ir ; Aparici:2009fh ; Liao:2016qyd .
The $N_{R}$SMEFT respects the SM gauge symmetry $SU(3)_{C}\times
SU(2)_{L}\times U(1)_{Y}$, and it is valid at energies above the electroweak
scale, set by the vacuum expectation value (VEV) of the Higgs, $v$. In the
considered case of the $N_{R}$LEFT as the low-energy limit of the
$N_{R}$SMEFT, there are no new particles at energies between $E$ and $v$, and
the $N_{R}$LEFT scale $\Lambda$ can be identified with $v$. At this scale, the
two EFTs must be matched.

At tree level, the matching of the $d=6$ LNC operators has been worked out in
Ref. Chala:2020vqp  and that of the $d=6$ LNV operators in Ref. Li:2020lba
.333For the $N_{R}$LEFT operators with two quarks, a charged lepton and a
neutrino, the tree-level matching has been performed also in Refs.
Dekens:2020ttz ; DeVries:2020jbs . For the first-generation charged fermions,
such operators contribute, among other processes, to neutrinoless double beta
decay Dekens:2020ttz ; Dekens:2021qch ; deVries:2022nyh ; Cirigliano:2022oqy .
Below we summarize the results of this matching relevant to the operators of
interest. The $N_{R}$SMEFT operators contributing to the matching relations
for the coefficients of the $N_{R}$LEFT operators from tables 1 and 2 are
given in tables 3 and 4, respectively.

| Name | Structure | $n_{N}=1$ | $n_{N}=3$
---|---|---|---|---
$d=6$ (LNC) | ${\cal O}_{dN}$ | $\left(\overline{d_{R}}\gamma_{\mu}d_{R}\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 9 | 81
${\cal O}_{uN}$ | $\left(\overline{u_{R}}\gamma_{\mu}u_{R}\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 9 | 81
${\cal O}_{QN}$ | $\left(\overline{Q}\gamma_{\mu}Q\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 9 | 81
${\cal O}_{HN}$ | $\left(H^{\dagger}i\overleftrightarrow{D_{\mu}}H\right)\left(\overline{N_{R}}\gamma^{\mu}N_{R}\right)$ | 1 | 9
$d=7$ (LNV) | ${\cal O}_{QNdH}$ | $\left(\overline{Q}N_{R}\right)\left(\overline{N_{R}^{c}}d_{R}\right)H$ | 18 | 162
${\cal O}_{dQNH}$ | $H^{\dagger}\left(\overline{d_{R}}Q\right)\left(\overline{N_{R}^{c}}N_{R}\right)$ | 18 | 108
${\cal O}_{QNuH}$ | $\left(\overline{Q}N_{R}\right)\left(\overline{N_{R}^{c}}u_{R}\right)\widetilde{H}$ | 18 | 162
${\cal O}_{uQNH}$ | $\widetilde{H}^{\dagger}\left(\overline{u_{R}}Q\right)\left(\overline{N_{R}^{c}}N_{R}\right)$ | 18 | 108
Table 3: Dimension-six and dimension-seven operators in the $N_{R}$SMEFT contributing to the tree-level matching relations for the Wilson coefficients of the dimension-six $N_{R}$LEFT operators from table 1. For each operator structure, we give the number of independent real parameters for $n_{N}=1$ and $n_{N}=3$ generations of $N_{R}$. The LNV operator structures require “+h.c.”. | Name | Structure | $n_{N}=1$ | $n_{N}=3$
---|---|---|---|---
$d=6$ (LNC) | ${\cal O}_{LNQd}$ | $\epsilon_{ab}\left(\overline{L^{a}}N_{R}\right)\left(\overline{Q^{b}}d_{R}\right)$ | 54 | 162
${\cal O}_{LdQN}$ | $\epsilon_{ab}\left(\overline{L^{a}}d_{R}\right)\left(\overline{Q^{b}}N_{R}\right)$ | 54 | 162
${\cal O}_{QuNL}$ | $\left(\overline{Q}u_{R}\right)\left(\overline{N_{R}}L\right)$ | 54 | 162
$d=7$ (LNV) | ${\cal O}_{dNLH}$ | $\epsilon_{ab}\left(\overline{d_{R}}\gamma_{\mu}d_{R}\right)\left(\overline{N_{R}^{c}}\gamma^{\mu}L^{a}\right)H^{b}$ | 54 | 162
${\cal O}_{uNLH}$ | $\epsilon_{ab}\left(\overline{u_{R}}\gamma_{\mu}u_{R}\right)\left(\overline{N_{R}^{c}}\gamma^{\mu}L^{a}\right)H^{b}$ | 54 | 162
${\cal O}_{QNLH1}$ | $\epsilon_{ab}\left(\overline{Q}\gamma_{\mu}Q\right)\left(\overline{N_{R}^{c}}\gamma^{\mu}L^{a}\right)H^{b}$ | 54 | 162
${\cal O}_{QNLH2}$ | $\epsilon_{ab}\left(\overline{Q}\gamma_{\mu}Q^{a}\right)\left(\overline{N_{R}^{c}}\gamma^{\mu}L^{b}\right)H$ | 54 | 162
${\cal O}_{NL1}$ | $\epsilon_{ab}\left(\overline{N_{R}^{c}}\gamma_{\mu}L^{a}\right)\left(iD^{\mu}H^{b}\right)\left(H^{\dagger}H\right)$ | 6 | 18
${\cal O}_{NL2}$ | $\epsilon_{ab}\left(\overline{N_{R}^{c}}\gamma_{\mu}L^{a}\right)H^{b}\left(H^{\dagger}i\overleftrightarrow{D^{\mu}}H\right)$ | 6 | 18
Table 4: Dimension-six and dimension-seven operators in the $N_{R}$SMEFT
contributing to the tree-level matching relations for the Wilson coefficients
of the dimension-six $N_{R}$LEFT operators from table 2. For each operator
structure, we give the number of independent real parameters for $n_{N}=1$ and
$n_{N}=3$ generations of $N_{R}$. All operator structures require “+h.c.”.

For the dimensionful Wilson coefficients of the $N_{R}$SMEFT operators, we use
capital $C_{i}^{(d)}$. They scale as $\Lambda_{\mathrm{NP}}^{4-d}$, where
$\Lambda_{\mathrm{NP}}\gg v$ is the scale of new physics. For simplicity, we
assume only one generation of $N_{R}$, while considering a general structure
of quark, charged lepton and LH neutrino flavor indices. We denote the quark
flavor indices by $i\,,j$ and the lepton flavor index by $\alpha$.
Furthermore, we assume that in the $N_{R}$LEFT, quarks are in the mass basis,
whereas the $N_{R}$SMEFT operators are written in the weak interaction basis.

The matching conditions for the coefficients of the LNC pair-$N_{R}$ operators
holding at the electroweak scale read:

$\displaystyle c_{dN,ij}^{V,RR}$
$\displaystyle=C_{dN}^{ij}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{d_{R}}^{ij}Z_{N}\,,\hskip
49.79231ptc_{uN,ij}^{V,RR}=C_{uN}^{ij}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{u_{R}}^{ij}Z_{N}\,,$
(4) $\displaystyle c_{dN,ij}^{V,LR}$
$\displaystyle=V_{ki}^{\ast}V_{lj}C_{QN}^{kl}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{d_{L}}^{ij}Z_{N}\,,\qquad
c_{uN,ij}^{V,LR}=C_{QN}^{ij}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{u_{L}}^{ij}Z_{N}\,.$
(5)

Here, $V$ is the CKM matrix defined through $d^{\prime}_{L}=Vd_{L}$, where
$d^{\prime}$ and $d$ are flavor and mass eigenstates, respectively. In
addition,

$g_{Z}\equiv\frac{e}{s_{W}c_{W}}\,,\qquad
Z_{\psi}^{ij}\equiv\left(T_{\psi}^{3}-Q_{\psi}s_{W}^{2}\right)\delta^{ij}\,,\qquad
Z_{N}\equiv-\frac{v^{2}}{2}C_{HN}\,,$ (6)

with $\psi$ denoting a SM fermion, and $T_{\psi}^{3}$ and $Q_{\psi}$ being its
weak isospin and electric charge, respectively. Further,
$s_{W}\equiv\sin\theta_{W}$ and $c_{W}\equiv\cos\theta_{W}$, where
$\theta_{W}$ is the weak mixing angle, and $m_{Z}$ is the mass of the $Z$
boson. Note that in $Z_{\psi}^{ij}$ we neglect the contributions of $d=6$
$N_{R}$SMEFT operators (cf. Eq. (2.18) in Ref. Li:2020lba ), since when
multiplied by $Z_{N}$ these would lead to $d=8$ effects.

For the LNC single-$N_{R}$ operators, we have:

$\displaystyle c_{d\nu N,ij\alpha}^{S,RR}$
$\displaystyle=V_{ki}^{\ast}\left(C_{LNQd}^{\alpha
kj}-\frac{1}{2}C_{LdQN}^{\alpha jk}\right),\qquad c_{d\nu
N,ij\alpha}^{T,RR}=-\frac{1}{8}V_{ki}^{\ast}C_{LdQN}^{\alpha jk}\,,$ (7)
$\displaystyle c_{u\nu N,ij\alpha}^{S,RR}$ $\displaystyle=0\,,\hskip
136.5733ptc_{u\nu N,ij\alpha}^{T,RR}=0\,,$ (8) $\displaystyle c_{d\nu
N,ij\alpha}^{S,LR}$ $\displaystyle=0\,,\hskip 136.5733ptc_{u\nu
N,ij\alpha}^{S,LR}=C_{QuNL}^{ji\alpha\,\ast}\,.$ (9)

We note that the Higgs exchange contributions (which a priori seem relevant)
are parametrically of the same order as $d=8$ terms, and thus, can be dropped
when working to $d=6$ Jenkins:2017jig .

In the $N_{R}$SMEFT at $d=6$, there is only one operator which violates lepton
number, but conserves baryon number, namely,
$\mathcal{O}_{NNNN}=(\overline{N_{R}^{c}}N_{R})(\overline{N_{R}^{c}}N_{R})$.
It vanishes identically in the case of one generation of $N_{R}$. Thus, the
$d=6$ LNV operators in the $N_{R}$LEFT are expected to get contributions from
$d=7$ LNV operators in the $N_{R}$SMEFT. A basis of the latter can be found in
Ref. Liao:2016qyd  (see also Ref. Bhattacharya:2015vja ). A subset of the
$d=7$ operators relevant for the matching we are interested in is collected in
the bottom part of tables 3 and 4.

For the LNV pair-$N_{R}$ operators, the matching conditions read:

$\displaystyle c_{dN,ij}^{S,RR}$
$\displaystyle=-\frac{v}{2\sqrt{2}}V_{ki}^{\ast}C_{QNdH}^{kj}\,,\qquad
c_{uN,ij}^{S,RR}=-\frac{v}{2\sqrt{2}}C_{QNuH}^{ij}\,,$ (10) $\displaystyle
c_{dN,ij}^{S,LR}$
$\displaystyle=\frac{v}{\sqrt{2}}V_{kj}C_{dQNH}^{ik}\,,\hskip
36.98866ptc_{uN,ij}^{S,LR}=\frac{v}{\sqrt{2}}C_{uQNH}^{ij}\,.$ (11)

Finally, for the LNV single-$N_{R}$ operators, we have:

$\displaystyle c_{d\nu N,ij\alpha}^{V,RR}$
$\displaystyle=-\frac{v}{\sqrt{2}}C_{dNLH}^{ij\alpha}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{d_{R}}^{ij}Z_{\nu
N}^{\alpha}\,,$ (12) $\displaystyle c_{u\nu N,ij\alpha}^{V,RR}$
$\displaystyle=-\frac{v}{\sqrt{2}}C_{uNLH}^{ij\alpha}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{u_{R}}^{ij}Z_{\nu
N}^{\alpha}\,,$ (13) $\displaystyle c_{d\nu N,ij\alpha}^{V,LR}$
$\displaystyle=-\frac{v}{\sqrt{2}}V_{ki}^{\ast}V_{lj}\left(C_{QNLH1}^{kl\alpha}-C_{QNLH2}^{kl\alpha}\right)-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{d_{L}}^{ij}Z_{\nu
N}^{\alpha}\,,$ (14) $\displaystyle c_{u\nu N,ij\alpha}^{V,LR}$
$\displaystyle=-\frac{v}{\sqrt{2}}C_{QNLH1}^{ij\alpha}-\frac{g_{Z}^{2}}{m_{Z}^{2}}Z_{u_{L}}^{ij}Z_{\nu
N}^{\alpha}\,.$ (15)

Here, the couplings $Z_{\nu N}^{\alpha}$ are defined as

$Z_{\nu
N}^{\alpha}\equiv\frac{v^{3}}{4\sqrt{2}}\left(C_{NL1}^{\alpha}+2C_{NL2}^{\alpha}\right).$
(16)

###  2.3 Running of the $N_{R}$LEFT operators of interest

Between the energy scale at which meson decays take place and the weak scale,
the Wilson coefficients of the $N_{R}$LEFT operators evolve according to
renormalization group equations (RGEs). These are analogous to the RGEs in the
LEFT (without $N_{R}$) Jenkins:2017dyc , and some partial results for the
operators with $N_{R}$ can be found in Refs. Chala:2020vqp ; DeVries:2020jbs ;
Li:2020lba . Below we provide the one-loop RGEs for the pair-$N_{R}$ and
single-$N_{R}$ operators summarized in tables 1 and 2. The operators given by
the product of vector currents are renormalized by QED corrections:444In these
equations, we ignore the contributions from the four-fermion leptonic
operators $\mathcal{O}_{eN}^{V,RR}$ and $\mathcal{O}_{eN}^{V,LR}$ as well as
from the double insertions of the $d=5$ neutrino dipole operators.

$\displaystyle\dot{c}_{qN,ij}^{V,RR}$
$\displaystyle=\dot{c}_{qN,ij}^{V,LR}=\frac{4}{3}e^{2}Q_{q}N_{c}\delta_{ij}\left[Q_{u}\left(c_{uN,kk}^{V,RR}+c_{uN,kk}^{V,LR}\right)+Q_{d}\left(c_{dN,kk}^{V,RR}+c_{dN,kk}^{V,LR}\right)\right],$
(17) $\displaystyle\dot{c}_{q\nu N,ij\alpha}^{V,RR}$
$\displaystyle=\dot{c}_{q\nu
N,ij\alpha}^{V,LR}=\frac{4}{3}e^{2}Q_{q}N_{c}\delta_{ij}\left[Q_{u}\left(c_{u\nu
N,kk\alpha}^{V,RR}+c_{u\nu N,kk\alpha}^{V,LR}\right)+Q_{d}\left(c_{d\nu
N,kk\alpha}^{V,RR}+c_{d\nu N,kk\alpha}^{V,LR}\right)\right].$ (18)

Here, $\dot{c}\equiv 16\pi^{2}\mu\dfrac{\mathrm{d}c}{\mathrm{d}\mu}$, with
$\mu$ being the renormalization scale; $q=d$ or $u$; $e$ is the QED coupling
constant, $Q_{q}$ is the electric charge of the quark $q$; and $N_{c}=3$ is
the number of colors. The flavor index $k$ runs over $u$ and $c$ for the up-
type quarks, and over $d,s,b$ for the down-type quarks. We see that only the
coefficients with $i=j$ run, whereas those with $i\neq j$, which will be of
interest to us in what follows, do not.

The operators of the scalar and tensor types are renormalized due to both QCD
and QED corrections:

$\displaystyle\dot{c}_{q}^{S}$
$\displaystyle=-6\left(g^{2}C_{F}+e^{2}Q_{q}^{2}\right)c_{q}^{S}\,,\qquad
c_{q}^{S}=c_{qN,ij}^{S,RR}\,,~{}c_{qN,ij}^{S,LR}\,,~{}c_{q\nu
N,ij\alpha}^{S,RR}\,,~{}c_{q\nu N,ij\alpha}^{S,LR}\,,$ (19)
$\displaystyle\dot{c}_{q}^{T}$
$\displaystyle=2\left(g^{2}C_{F}+e^{2}Q_{q}^{2}\right)c_{q}^{T}\,,\qquad
c_{q}^{T}=c_{qN,ij}^{T,RR}\,,~{}c_{q\nu N,ij\alpha}^{T,RR}\,,$ (20)

where $g$ is the strong coupling constant, and
$C_{F}=(N_{c}^{2}-1)/(2N_{c})=4/3$.

Solving these RGEs to the leading-logarithm approximation, we find:

$\displaystyle c_{q}^{S}(\mu)$
$\displaystyle=\left[\frac{g(\mu)}{g(\mu_{1})}\right]^{\frac{6C_{F}}{b_{g}}}\left[\frac{e(\mu)}{e(\mu_{1})}\right]^{\frac{6Q_{q}^{2}}{b_{e}}}c_{q}^{S}(\mu_{1})\,,$
(21) $\displaystyle c_{q}^{T}(\mu)$
$\displaystyle=\left[\frac{g(\mu)}{g(\mu_{1})}\right]^{-\frac{2C_{F}}{b_{g}}}\left[\frac{e(\mu)}{e(\mu_{1})}\right]^{-\frac{2Q_{q}^{2}}{b_{e}}}c_{q}^{T}(\mu_{1})\,,$
(22)

where $\mu_{1}<\mu$ are two energy scales, and $b_{g}$ and $b_{e}$ are the
coefficients of the one-loop QCD and QED beta functions, respectively:555For
the beta functions, we use $\dot{g}=-b_{g}g^{3}$ and $\dot{e}=-b_{e}e^{3}$.

$\displaystyle b_{g}$
$\displaystyle=\frac{11}{3}N_{c}-\frac{2}{3}\left(n_{u}+n_{d}\right),$ (23)
$\displaystyle b_{e}$
$\displaystyle=-\frac{4}{3}\left(Q_{e}^{2}n_{e}+Q_{d}^{2}N_{c}n_{d}+Q_{u}^{2}N_{c}n_{u}\right).$
(24)

Setting $\mu_{1}=m_{b}$ and $\mu=v$, we find

$\displaystyle c_{u}^{S}(v)$ $\displaystyle=0.68\,c_{u}^{S}(m_{b})\,,\qquad
c_{d}^{S}(v)=0.69\,c_{d}^{S}(m_{b})\,,$ (25) $\displaystyle c_{u}^{T}(v)$
$\displaystyle=1.14\,c_{u}^{T}(m_{b})\,,\qquad
c_{d}^{T}(v)=1.13\,c_{d}^{T}(m_{b})\,.$ (26)

Thus, when running from $m_{b}$ to $v$, the scalar type operators get
suppressed, whereas the tensor type operator get enhanced, cf. Ref. Li:2020lba
. Finally, we note that the contribution to the running from QED corrections
is practically negligible.

Concerning the running in the $N_{R}$SMEFT, i.e. between the EW scale and the
scale of new physics, $\Lambda_{\mathrm{NP}}$, the one-loop RGEs for $d=6$
operators with $N_{R}$ have been derived in Refs. Chala:2020pbn ;
Datta:2020ocb ; Datta:2021akg , while partial results for some $d=7$
$N_{R}$-operators (relevant for neutrinoless double beta decay) have been
obtained in Ref. Dekens:2020ttz .

###  2.4 An example UV completion for $d=7$ $N_{R}$SMEFT

The operators including one or two $N_{R}$ in the $N_{R}$SMEFT can be
generated in the ultra-violet (UV) in a variety of ways. For example, Ref.
Bischer:2019ttk  discusses all possible leptoquark states that can induce
$d=6$ $N_{R}$SMEFT operators. Also Refs. Cottin:2021lzz  and Beltran:2021hpq
study some example models for pair-$N_{R}$ and single-$N_{R}$ operators at
$d=6$, respectively, including leptoquark and $Z^{\prime}$ models. However,
neither of these references gives example models for the $d=7$ $N_{R}$SMEFT
operators, to which our low-energy LNV $N_{R}$ operators must be matched.

Here, we do not attempt to give a systematic list of such models for the $d=7$
operators. Instead, we will briefly discuss one concrete, but very typical UV
model as an example. This example model generates the operator ${\cal
O}_{QNdH}$. For the other $d=7$ operators containing one or two $N_{R}$, other
UV models can be constructed in an analogous manner.

Our example model contains two scalar leptoquark states,
$S_{1}\sim(\mathbf{3},\mathbf{1})_{-1/3}$ and
$S_{2}\sim(\mathbf{3},\mathbf{2})_{1/6}$, with the numbers inside parentheses
denoting the $(SU(3)_{C},SU(2)_{L})$ representations and the subscript
standing for the hypercharge. Apart from the mass and kinetic terms for the
leptoquark states, the model allows to write down the following interactions,
in addition to the SM Lagrangian:

$\mathcal{L}\supset-\left[Y_{QN}\overline{Q}N_{R}S_{2}+Y_{Nd}\overline{N_{R}^{c}}d_{R}S_{1}^{\dagger}+\mu
HS_{2}^{\dagger}S_{1}+{\rm h.c.}\right]+\cdots,$ (27)

where, for simplicity, we have suppressed generation indices. Note that this
Lagrangian violates lepton number, but lepton number is recovered in the limit
$\mu\rightarrow 0$ Hirsch:1996qy . Integrating out both leptoquarks will
produce ${\cal O}_{QNdH}$ with $C_{QNdH}=(\mu
Y_{QN}Y_{Nd})/(m_{S_{1}}^{2}m_{S_{2}}^{2})$. In the limit
$\mu=m_{S_{1}}=m_{S_{2}}=\Lambda_{\mathrm{NP}}$, this gives the correct
$C_{QNdH}\propto\Lambda_{\mathrm{NP}}^{-3}$ dependence for a $d=7$ operator,
as expected. Clearly, this model will also generate $d=6$ operators with
$N_{R}$. For example, the term proportional to $Y_{QN}$ ($Y_{Nd}$) times its
hermitian conjugate will generate ${\cal O}_{QN}$ (${\cal O}_{dN}$). However,
these are only examples, as actually more $d=6$ operators can be generated in
this model; note the dots in Eq. (27).

##  3 HNL production in meson decays

We are interested in HNLs produced in meson decays. In the framework of the
$N_{R}$LEFT, such decays are triggered by both active-heavy neutrino mixing
and higher-dimensional operators. The HNL production via mixing has been
studied in detail in the literature, see e.g. Refs. Bondarenko:2018ptm ;
Ballett:2019bgd ; Coloma:2020lgy  for recent reappraisals. The HNL production
via effective interactions has been considered in Ref. DeVries:2020jbs  for
the four-fermion operators with two quarks, one charged lepton and $N_{R}$. In
this section, we focus on the meson decay modes triggered by the operators
with two quarks and two neutral leptons (either a pair of $N_{R}$, or $N_{R}$
and $\nu_{L}$) discussed in section 2.

Given their large production rates at the LHC, we focus on pseudoscalar mesons
containing either a $c$ or $b$ quark and a light quark, i.e. $D$\- and
$B$-mesons. Vector mesons with the same quark content have much smaller
lifetimes, and thus, their contribution in the HNL production is sub-dominant
DeVries:2020jbs .

In figure 1, we depict two- and three-body pseudoscalar meson decays triggered
by the effective operators of interest.

Figure 1: Two- and three-body pseudoscalar meson decays triggered by the
effective operators with a pair of HNLs, or one HNL and one SM neutrino. The
blue blob denotes the insertion of an effective operator. In the right
diagram, the quark $q$ can be $u$, $d$, $s$, $c$ or $b$, whereas the quarks
$q_{i}$ and $q_{j}$ are both of either the up- or down-type, with
$m_{q_{j}}>m_{q_{i}}$. $P^{(\prime)}$ and $V$ denote pseudoscalar and vector
mesons, respectively.

The same $N_{R}$LEFT operator induces a leptonic two-body decay as well as
contributes to a number of semi-leptonic three-body decays. In the diagram
corresponding to the three-body decay, $q$ can be any quark but top, since the
latter does not form bound states and is not present in the spectrum of the
$N_{R}$LEFT. For the transition $q_{j}\to q_{i}$ with $i\neq j$, both quarks
have to be of either the up- or down-type, since the operators considered
involve quark bilinears with zero electric charge. In addition,
$m_{q_{j}}>m_{q_{i}}$ for the decay to be kinematically allowed.666This
condition is necessary, but not sufficient, since a meson mass is not a simple
sum of the masses of its constituent quarks. These considerations lead to the
following four possibilities: $c\to u$ (for the up-quark sector), and $b\to
d$, $b\to s$ and $s\to d$ (for the down-quark sector). Below we discuss each
of them in detail.

  * •

The operators with $q_{j}=c$ and $q_{i}=u$ lead to the leptonic two-body
decays $D^{0}\to NN$ and $D^{0}\to\nu N$.777In the case of Dirac $\nu$ and
$N$, these decays should read $D^{0}\to N\overline{N}$ and
$D^{0}\to\nu\overline{N}$ ($D^{0}\to\overline{\nu}N$). In addition, they
generate the following transitions resulting in semi-leptonic three-body
decays:
$D^{0}\to\pi^{0}\,,~{}\eta\,,~{}\eta^{\prime}~{}(\rho^{0}\,,~{}\omega)$ for
$q=u$; $D^{+}\to\pi^{+}~{}(\rho^{+})$ for $q=d$; $D_{s}^{+}\to
K^{+}~{}(K^{\ast+})$ for $q=s$;
$\eta_{c}\to\overline{D^{0}}~{}(\overline{D^{\ast 0}})$ for $q=c$; and
$B_{c}^{+}\to B^{+}~{}(B^{\ast+})$ for $q=b$. Here, the final states outside
(inside) parentheses are pseudoscalar (vector) mesons. In what follows, we
will not take into account the decays of $\eta_{c}$ and $B_{c}^{+}$. The
$\eta_{c}$ meson has a lifetime of $\mathcal{O}(10^{-23})$ s Workman:2022ynf ,
which is much smaller than that of other pseudoscalar mesons and comparable
with the lifetimes of vector mesons. Concerning $B_{c}^{+}$, its production
rate at the LHC is much smaller than that of $D$\- and $B$-mesons, see e.g.
Ref. LHCb:2019tea .

  * •

The operators with $q_{j}=b$ and $q_{i}=d$ 888In figure 1, we need to send
$\overline{q}\to q$ and $q_{j,i}\to\overline{q_{j,i}}$ to match the particle
content of $B$-mesons. induce $B^{0}\to NN$ and $B^{0}\to\nu N$ decays, as
well as the following transitions: $B^{+}\to\pi^{+}~{}(\rho^{+})$ for $q=u$;
$B^{0}\to\pi^{0}\,,~{}\eta\,,~{}\eta^{\prime}~{}(\rho^{0}\,,~{}\omega)$ for
$q=d$; $B_{s}^{0}\to\overline{K^{0}}~{}(\overline{K^{\ast 0}})$ for $q=s$;
$B_{c}^{+}\to D^{+}~{}(D^{\ast+})$ for $q=c$; and
$\eta_{b}\to\overline{B^{0}}~{}(\overline{B^{\ast 0}})$ for $q=b$. Again, we
will not consider the decays of $B_{c}^{+}$ as well as those of $\eta_{b}$.

  * •

The operators with $q_{j}=b$ and $q_{i}=s$ give rise to $B_{s}^{0}\to NN$ and
$B_{s}^{0}\to\nu N$ decays. They also contribute to $B^{+}\to
K^{+}~{}(K^{\ast+})$ for $q=u$; $B^{0}\to K^{0}~{}(K^{\ast 0})$ for $q=d$;
$B_{s}^{0}\to\eta\,,~{}\eta^{\prime}~{}(\phi)$ for $q=s$; $B_{c}^{+}\to
D_{s}^{+}~{}(D_{s}^{\ast+})$ for $q=c$; and
$\eta_{b}\to\overline{B_{s}^{0}}~{}(\overline{B_{s}^{\ast 0}})$ for $q=b$.

  * •

The operators with $q_{j}=s$ and $q_{i}=d$ trigger $K^{0}\to NN$ and
$K^{0}\to\nu N$ decays, as well as the following transitions:
$K^{+}\to\pi^{+}$ for $q=u$; $K^{0}\to\pi^{0}$ for $q=d$;
$\eta^{\prime}\to\overline{K^{0}}~{}(\overline{K^{\ast 0}})$ for $q=s$;
$D_{s}^{+}\to D^{+}$ for $q=c$; and
$\overline{B_{s}^{0}}\to\overline{B^{0}}~{}(\overline{B^{\ast 0}})$ for $q=b$.
We will not consider these decay modes in what follows. This is because for
the listed $D_{s}^{+}$ and $\overline{B_{s}^{0}}$ decays, the available mass
window is very small, and for the light mesons’ decays (kaons and
$\eta^{\prime}$), the numerical predictions from Pythia8 Sjostrand:2014zea ;
Bierlich:2022pfr , which we will use for numerical simulation in this work, do
not agree well with the LHCf experiment pythiaForward .

We have computed the corresponding partial decay widths assuming both $N$ and
$\nu$ are (i) Dirac and (ii) Majorana particles. In appendix A, we provide the
formulae for the two-body decay widths and the three-body decay amplitudes in
both cases. To compute the three-body decay widths, we have used the procedure
described in Ref. DeVries:2020jbs . We summarize it in appendix A.2 for
completeness.

If the products of HNL decays are not detected, the modes we consider here
contribute to the $P\to~{}\text{inv.}$ or $P\to P^{\prime}/V\nu\overline{\nu}$
decays. The existing upper limits on the branching ratios of such decays are
listed in table 5. In what follows, we will take these constraints into
account.

Decay | Limit on BR | Decay | Limit on BR | Decay | Limit on BR
---|---|---|---|---|---
$D^{0}\to\text{inv.}$ | $9.4\times 10^{-5}$ | $B^{0}\to\text{inv.}$ | $2.4\times 10^{-5}$ | $B_{s}^{0}\to\phi\nu\overline{\nu}$ | $5.4\times 10^{-3}$
|  | $B^{0}\to\pi^{0}\nu\overline{\nu}$ | $9.0\times 10^{-6}$ | $B^{0}\to K^{0}\nu\overline{\nu}$ | $2.6\times 10^{-5}$
|  | $B^{0}\to\rho^{0}\nu\overline{\nu}$ | $4.0\times 10^{-5}$ | $B^{0}\to K^{\ast 0}\nu\overline{\nu}$ | $1.8\times 10^{-5}$
|  | $B^{+}\to\pi^{+}\nu\overline{\nu}$ | $1.4\times 10^{-5}$ | $B^{+}\to K^{+}\nu\overline{\nu}$ | $1.6\times 10^{-5}$
|  | $B^{+}\to\rho^{+}\nu\overline{\nu}$ | $3.0\times 10^{-5}$ | $B^{+}\to K^{\ast+}\nu\overline{\nu}$ | $4.0\times 10^{-5}$
Table 5: Upper limits on the branching ratios (BRs) of (semi-)invisible $D$\-
and $B$-meson decays extracted from Ref. Workman:2022ynf . The limits
highlighted in bold provide the most stringent constraints on the
corresponding Wilson coefficients, see section 5.

In figure 2, we show the branching ratios of $D$\- and $B$-meson decays
induced by the LNC pair-$N_{R}$ operators $\mathcal{O}_{uN,12}^{V,RR}$ (top),
$\mathcal{O}_{dN,31}^{V,RR}$ (middle) and $\mathcal{O}_{dN,32}^{V,RR}$
(bottom), assuming $N$ is Majorana.

Figure 2: Branching ratios of $D$\- and $B$-meson decays triggered by the LNC
pair-$N_{R}$ operators $\mathcal{O}_{uN,12}^{V,RR}$ (top),
$\mathcal{O}_{dN,31}^{V,RR}$ (middle) and $\mathcal{O}_{dN,32}^{V,RR}$
(bottom). The corresponding Wilson coefficient has been set to
$10^{-3}v^{-2}$.

We turn on one operator at a time, setting the corresponding Wilson
coefficient to $10^{-3}v^{-2}$. This choice allows us to avoid the bounds from
(semi-)invisible meson decays summarized in table 5. We see that the two-body
decay dominates over the three-body ones for $m_{N}\gtrsim 325$ MeV in the
case of $\mathcal{O}_{uN,12}^{V,RR}$, for $m_{N}\gtrsim 1.2$ GeV in the case
of $\mathcal{O}_{dN,31}^{V,RR}$, and for $m_{N}\gtrsim 1.05$ GeV in the case
of $\mathcal{O}_{dN,32}^{V,RR}$. For small $m_{N}$, the two-body decay is
suppressed, since the corresponding decay width is proportional to
$m_{N}^{2}$, see Eq. (42) in appendix A. Turning on the LR operators
$\mathcal{O}_{uN,12}^{V,LR}$, $\mathcal{O}_{dN,31}^{V,LR}$, and
$\mathcal{O}_{dN,32}^{V,LR}$ (one at a time) would lead to the same results as
can be seen from Eqs. (42), (53) and (54).

In figure 3, we display the branching ratios of the same meson decays, but
this time triggered by the LNV pair-$N_{R}$ operators
$\mathcal{O}_{uN,12}^{S,RR}$ (top), $\mathcal{O}_{dN,31}^{S,RR}$ (middle), and
$\mathcal{O}_{dN,32}^{S,RR}$ (bottom).

Figure 3: Branching ratios of $D$\- and $B$-meson decays triggered by the LNV
pair-$N_{R}$ operators $\mathcal{O}_{uN,12}^{S,RR}$ (top),
$\mathcal{O}_{dN,31}^{S,RR}$ (middle) and $\mathcal{O}_{dN,32}^{S,RR}$
(bottom). The corresponding Wilson coefficient has been set to
$10^{-3}v^{-2}$.

Again, we turn on one operator at a time and set the corresponding Wilson
coefficient to $10^{-3}v^{-2}$.999In the case of non-hermitian operators, as
e.g. $\mathcal{O}_{qN}^{S,RR}$, $q=u$ or $d$, “one operator at a time” means
$c_{qN,ij}^{S,RR}\neq 0$, whereas $c_{qN,ji}^{S,RR}=0$. We note also that
$c_{dN,31}^{S,RR}=10^{-3}v^{-2}$ violates slightly the limit on
$\mathrm{BR}(B^{0}\to\text{inv.})$ from table 5. However, we opt for this
value of the Wilson coefficient to be uniform with the other plots. Contrary
to the LNC operators, now the leptonic two-body decays dominate for all
kinematically available values of $m_{N}$, since there is no helicity
suppression from $m_{N}$, cf. Eq. (42). Switching on the LR operators
$\mathcal{O}_{uN,12}^{S,LR}$, $\mathcal{O}_{dN,31}^{S,LR}$, and
$\mathcal{O}_{dN,32}^{S,LR}$ (one at a time) would lead to the same results as
can be inferred from Eqs. (42), (53), and (54). The tensor operators
$\mathcal{O}_{dN}^{T,RR}$ and $\mathcal{O}_{uN}^{T,RR}$ vanish identically for
one generation of $N_{R}$. For this reason, we do not consider them in what
follows.

Concerning single-$N_{R}$ operators, we have also studied their contribution
to two- and three-body $D$\- and $B$\- meson decays, showing some results in
appendix B. In particular, we present in figure 8 the branching ratios of
$B$-meson decays induced by the LNC operators $\mathcal{O}_{d\nu
N,31\alpha}^{S,RR}$ and $\mathcal{O}_{d\nu N,31\alpha}^{T,RR}$ and the LNV
operator $\mathcal{O}_{d\nu N,31\alpha}^{V,RR}$. We observe that our results
are analogous to the ones obtained in Ref. DeVries:2020jbs , cf. figure 2
therein, for single-$N_{R}$ operators including a charged lepton. Therefore,
we expect the constraints derived from the far detectors to be essentially the
same for the neutral-current single-$N_{R}$ operators, and we hence do not
include the latter in our numerical simulation.

Before closing the section, we mention that with only one $N_{R}$ generation,
pair-$N_{R}$ operators cannot make the HNL decay by themselves, in contrast to
single-$N_{R}$ operators. In appendix C, we have computed the partial decay
widths of HNLs into a light meson and an active neutrino triggered by the
neutral single-$N_{R}$ operators. Since we only simulate HNLs produced by
pair-$N_{R}$ operators and only turn on one $N_{R}$LEFT operator at a time, we
consider that the HNL decay is exclusively due to active-heavy neutrino
mixing. The partial decay widths of HNL via mixing are computed according to
the formulae given in Refs. Atre:2009rg ; Bondarenko:2018ptm ; DeVries:2020jbs
.

##  4 Experiments and details of simulation

In this work, we focus on a series of far-detector experiments proposed or
even approved at the LHC. These programs were mostly proposed in recent years
and include ANUBIS Bauer:2019vqk , AL3X Gligorov:2018vkc , CODEX-b
Gligorov:2017nwh , FACET Cerci:2021nlb , FASER and FASER2 Feng:2017uoz ;
FASER:2018eoc , MAPP1 and MAPP2 Pinfold:2019nqj ; Pinfold:2019zwp , and
MATHUSLA Chou:2016lxi ; Curtin:2018mvb ; MATHUSLA:2020uve . Since they are
supposed to be operated at different interaction points (IPs) and time scales
of the LHC, the projected integrated luminosities also vary, ranging from as
low as 30 fb-1 for MAPP1, up to the full HL-LHC target, 3 ab-1 for ANUBIS,
FACET, FASER2, and MATHUSLA. Moreover, typically these far detectors are
located with a distance of $\mathcal{O}(1)-\mathcal{O}(100)$ m to their
corresponding IP, allowing for sufficient amount of rock, lead, or other
shielding material to effectively remove the potential SM background events
stemming from the IP. As a result, vanishing background is usually assumed for
phenomenological studies on these experiments.

Here, we study HNLs pair-produced at the LHC with the center-of-mass energy
$\sqrt{s}=14$ TeV from charm and bottom mesons’ rare decays mediated by a
$N_{R}$LEFT operator. These HNLs then travel a certain distance before
decaying via their mixing with the SM active neutrino. If the boosted decay
length of the HNLs falls into a suitable range, the HNLs have a high
probability of decaying inside some of the far detectors. Since the HNLs
mostly decay to channels that are not fully invisible (three active
neutrinos), the decay vertex can be reconstructed by these experiments.

While we analytically compute the production and decay rates of the HNLs, we
rely on a Monte-Carlo simulation tool, Pythia8.3 Bierlich:2022pfr ;
Sjostrand:2014zea , to estimate the decay probabilities of the HNLs inside
each of the far detectors. Depending on the flavor indices of the EFT operator
we choose to turn on, Pythia8 can generate either the $pp\to c\bar{c}$ or
$pp\to b\bar{b}$ process inclusively for $D$\- or $B$-meson production at the
LHC, respectively. These mesons then decay to different channels including at
least two HNLs, with relative decay branching ratios implemented according to
the analytical computation detailed in section 3. This allows for maximal
usage of all the simulated data while respecting the relative importance of
each decay channel. Since Pythia8 provides the kinematics of each simulated
HNL, we can compute its decay probability in each detector according to the
geometry and position with respect to the IP. We then take the average of the
decay probabilities of all the simulated HNLs, to obtain the expected
acceptance rate at each detector, $\langle P[N\text{ decay}]\rangle$. This
allows us to compute the projected number of signal events, $N_{S}$, with the
following equation:

$\displaystyle N_{S}=\sum_{i}2\cdot N_{M_{i}}\cdot\text{BR}(M_{i}\to
NN+\text{anything})\cdot\langle P[N\text{
decay}]\rangle\cdot\text{BR}(N\to\text{vis.}),$ (28)

where the factor 2 accounts for the fact that each decay of the mesons into
the HNLs via the $N_{R}$LEFT operator includes two HNLs, $N_{M_{i}}$ is the
number of the mother mesons, and $\text{BR}(N\to\text{vis.})$ denotes the
decay branching ratio of the HNLs into visible final states. For our study, we
consider only the decay channel into three active neutrinos to be invisible.

In table 6, we reproduce the inclusive production rates of all the relevant
types of heavy mesons from Ref. DeVries:2020jbs , which were obtained therein
combining both LHCb measurements and extrapolation with the state-of-the-art
simulation tools FONLL Cacciari:1998it ; Cacciari:2001td  and Pythia8. Note
that in table 6 we have in addition multiplied the production rates of these
heavy mesons given in Ref. DeVries:2020jbs  by a factor of 1.08 and 1.06 for
the $B$\- and $D$-mesons, respectively, in order to take into account the fact
that we consider $\sqrt{s}=14$ TeV instead of 13 TeV Hirsch:2020klk .

$D^{0}$ | $D^{\pm}$ | $D^{\pm}_{s}$ | $B^{0}$ | $B^{\pm}$ | $B_{s}^{0}$
---|---|---|---|---|---
$4.12\times 10^{16}$ | $2.16\times 10^{16}$ | $7.02\times 10^{15}$ | $1.58\times 10^{15}$ | $1.58\times 10^{15}$ | $2.73\times 10^{14}$
Table 6: Inclusive production numbers of $D$\- and $B$-mesons for an
integrated luminosity of 3 ab-1, extracted from Ref. DeVries:2020jbs , but
multiplied by a factor of 1.08 (1.06) for the $B$-mesons ($D$-mesons) obtained
with Pythia8 to account for the increase in production cross section of the
heavy mesons from the center-of-mass energy 13 TeV to 14 TeV Hirsch:2020klk .

In principle, the computation of the decay probability of each HNL can be
summarized with a simple formula, if the flight path of the HNL traverses the
detector:

$\displaystyle P[N\text{ decay}]=\text{e}^{-L_{1}/\beta\gamma
c\tau}-\text{e}^{-L_{2}/\beta\gamma c\tau},$ (29)

where $\beta$, $\gamma$, and $\tau$ are the speed, Lorentz boost factor, and
the proper lifetime of the HNL, $c$ is the speed of light, and $L_{1}$ and
$L_{2}$ label respectively the distance from the IP (where we assume the heavy
mesons decay into the HNLs instantly) to the incoming and outgoing sides of
the detector. This formula holds perfectly for a far detector with a
spherical-shell shape facing towards the IP. However, since the far-detector
experiments mostly have a more complicated geometrical shape and relative
orientation with respect to the IP, more sophisticated formulae are required
for more accurate estimate of the decay probability. In this paper, we
strictly follow Ref. DeVries:2020jbs  for computing the decay probabilities in
each far detector with the corresponding geometrical configuration. As
exception, FACET Cerci:2021nlb  was not included in Ref. DeVries:2020jbs , so
here, we give a brief description of this detector. Its concept is essentially
rather similar to that of FASER(2) and AL3X. It is a proposed cylindrical
detector in the forward direction of the CMS experiment, with a distance of
101 m between the CMS IP and its near end. Further, it has a radius of 0.5 m
and length 18 m. However, it does not cover the pseudorapidity up to infinity
as FASER(2) does; instead, it covers the polar angle between 1 and 4 mrad, or
equivalently, the pseudorapidity between 6.2 and 7.6. We require that the HNLs
travel towards the far end of the FACET decay chamber (but not to the sides).

Before closing, we should also mention that since Pythia8 is not well
validated for heavy meson production in the very forward direction, we use
FONLL to correct the differential production cross section of mesons for
various ranges of small transverse momentum and large pseudorapidity, for
three experiments: FACET, FASER, and FASER2.

##  5 Results of numerical analysis

We present and discuss the numerical results in this section. We restrict
ourselves to Majorana HNLs, and study the LNC and LNV pair-$N_{R}$ operators
in the $N_{R}$LEFT triggering the meson decays described in section 3. More
specifically, the set of the LNC operators we consider is
$\mathcal{O}_{uN,12}^{V,RR}$, $\mathcal{O}_{dN,31}^{V,RR}$, and
$\mathcal{O}_{dN,32}^{V,RR}$, and the set of the LNV operators is
$\mathcal{O}_{uN,12}^{S,RR}$, $\mathcal{O}_{dN,31}^{S,RR}$, and
$\mathcal{O}_{dN,32}^{S,RR}$.101010We recall that the corresponding LR
operators, if switched on one at a time, would lead to the same results. We
first show in figures 4 and 5 sensitivity reaches of the various far detectors
we consider, in the plane $|V_{lN}|^{2}$ vs. $m_{N}$, for the LNC and LNV
operators, respectively. These isocurves are for three signal events,
corresponding to 95% confidence level (C.L.) limits, given zero background
events. We switch on one operator only each time and fix the corresponding
Wilson coefficient to $10^{-3}v^{-2}$, where $v=246$ GeV is the SM Higgs VEV.
This is approximately 700 times smaller than the Fermi constant and, for the
LNC operators, corresponds to a NP scale $\Lambda_{\mathrm{NP}}$ of 7.78 TeV
for $c\sim 1/\Lambda_{\mathrm{NP}}^{2}$. For this value of the Wilson
coefficient, the limits coming from the $D^{0}$ and $B^{0}$ invisible decays,
as well as from $P\to P^{\prime}/V\nu\overline{\nu}$, are avoided (see section
3). In our numerical simulation, we consider the HNLs only dominantly mixing
with either the electron neutrino or the muon neutrino: $l=e,\mu$. We do not
consider HNLs decaying to $\tau$’s, because we expect the sensitivities to be
much reduced, compared to $e,\mu$.

Figure 4: Exclusion limits in the plane $|V_{lN}|^{2}$ vs. $m_{N}$, for the
three LNC pair-$N_{R}$ operators we consider, with the corresponding Wilson
coefficients fixed to $10^{-3}v^{-2}$. Figure 5: The same plots as Fig. 4, but
for the three LNV pair-$N_{R}$ operators.

In figure 4, we observe that in the top-right and bottom-left plots, there is
stronger reach in $m_{N}$ than in the top-left plot. This is due to the larger
masses of $B$-mesons than those of $D$-mesons and hence higher kinematic
thresholds. Moreover, the $B$-meson scenarios show better exclusion limits in
$|V_{lN}|^{2}$ as well, compared to the $D$-meson one, despite the larger
production rates of $D$-mesons than those of $B$-mesons by about one order of
magnitude, cf. table 6. The main reason is that the decay branching ratios of
$B$-mesons into the HNLs are in general larger than those of $D$-mesons by
about two orders of magnitude, cf. figure 2, assuming equal values of the
corresponding Wilson coefficients.

Among these far-detector experiments, we find the best limits come from
MATHUSLA, which can probe the values of $|V_{lN}|^{2}$ down to $10^{-12}$
($10^{-15}$) in the charm (bottom) scenario(s), followed by ANUBIS and AL3X,
while FASER shows the weakest sensitivities, mainly because of its small
volume. We note that the approved experiment MAPP1 with only 30 fb-1 of
integrated luminosity can already compete with the future experiment FASER2
with the full HL-LHC data, especially in the bottom scenarios. The reason is
its larger size and smaller distance to the corresponding IP. Also, CODEX-b,
MAPP2, and FACET can exclude almost the same parts of the parameter space.

The plots in figure 5 are very similar to those in figure 4, with the main
difference being that the exclusion limits in $|V_{lN}|^{2}$ are slightly
stronger because the decay branching ratios of the mesons into a pair of HNLs
are larger in the LNV case than in the LNC case, given the same value of the
Wilson coefficients, see figures 2 and 3.

The best present limits for HNL masses below 5 GeV that can be found in the
literature come from different experiments such as NA62 NA62:2020mcv , T2K
T2K:2019jwa , CHARM Bergsma:1985is , PS191 Bernardi:1987ek , JINR
Baranov:1992vq , Belle Liventsev:2013zz , and DELPHI Abreu:1996pa . These
limits are around $|V_{lN}|^{2}\sim 10^{-9}$–$10^{-5}$, depending on the HNL
mass, and correspond to a minimal HNL scenario in which both the production
and the decay of the HNL are determined by their mixing with the active
neutrinos. In our theoretical setup, however, the production of the HNLs
receives additional contributions from effective operators affecting the HNL
sensitivity of these experiments. A detailed reinterpretation of the current
HNL limits coming from different experiments in the $N_{R}$LEFT depends on
specific details of each particular experiment, and it is beyond the scope of
this current work. For this reason, we do not show these bounds in figures 4
and 5.

Figure 6: Exclusion limits in the plane $c$ vs. $m_{N}$, for the three LNC
pair-$N_{R}$ operators we consider, with $|V_{lN}|^{2}$ fixed at $10^{-7}$ or
$10^{-10}$.

Figure 7: The same plots as Fig. 6, but for the three LNV pair-$N_{R}$
operators.

We then present in figures 6 and 7 the sensitivity curves in the plane $c$ vs.
$m_{N}$, for two fixed values of active-heavy mixing, $|V_{lN}|^{2}=10^{-7}$
and $|V_{lN}|^{2}=10^{-10}$. In order to obtain these plots, we select the
chosen benchmark mixing parameters and find the value of the Wilson
coefficient in each scenario such that the predicted number of signal events
is equal to three. We find the same hierarchy in sensitivities of these
experiments as in figures 4 and 5, but the final reaches in $c$ differ between
the two chosen $|V_{lN}|^{2}$ values. This is because for any fixed HNL mass,
varying the mixing parameter squared would change the proper lifetime of the
HNL and hence the detectors’ acceptance to these HNLs. $|V_{lN}|^{2}=10^{-7}$
corresponds to an acceptance rate roughly three orders of magnitude better
than that given by $|V_{lN}|^{2}=10^{-10}$, in the large decay length limit
($\beta\gamma c\tau\gg L_{1},L_{2}$, as in Eq. (29)). Therefore, the lower
sensitivity reach in $c$ for $|V_{lN}|^{2}=10^{-10}$ is roughly weaker than
that for $|V_{lN}|^{2}=10^{-7}$, by a factor of $\sqrt{1000}\approx 32$, since
the HNL production rates are proportional to $c^{2}$. In the charm scenario,
for the larger mixing parameter case, we find the best experiments (MATHUSLA,
ANUBIS, and AL3X) can reach $c$ as low as about $10^{-4}$ TeV-2, while the
other weaker experiments can still probe $c$ down to the order of $10^{-3}$
TeV-2 except FASER, or even better by about one order of magnitude in the two
bottom scenarios. In the LNC case, $c\sim 10^{-4}~{}(10^{-5})$ TeV-2
corresponds roughly to a new physics scale $\Lambda_{\mathrm{NP}}$ of 100
(316) TeV, if we take $c\sim 1/\Lambda_{\mathrm{NP}}^{2}$ and ignore the
relatively minor effects of QCD running.

Again as in figures 4 and 5, we find the LNV results in figure 7 are slightly
stronger than the LNC ones in figure 6, for the same reason as stated above.
Here, the conversion from Wilson coefficient to $\Lambda_{\mathrm{NP}}$ is
more subtle than the LNC operators’ case, because the $d=6$ LNV operators at
the scale below $v=246$ GeV are matched to the $d=7$ operators that include a
Higgs boson in the $N_{R}$SMEFT, cf. Eq. (10). Therefore, we estimate that for
$c\sim 10^{-4}~{}(10^{-5})$ TeV-2 in the low energy scale, the corresponding
$\Lambda_{\mathrm{NP}}$ is about 10 (21) TeV, with $c\sim
v/(2\sqrt{2}\Lambda_{\mathrm{NP}}^{3})$, if we again ignore the QCD running
effects. This is weaker than the LNC case by about one order of magnitude.

The gray shaded regions in figures 6 and 7 correspond to the limits coming
from invisible and semi-invisible $D$\- and $B$-meson decays collected in
table 5. More specifically, in the $c\to u$ scenario, it is
$D^{0}\to\text{inv.}$ which provides the only constraint. In the $b\to d$
scenario and for the LNC operator (middle-row plots in figure 6), there is an
interplay between $B^{0}\to\pi^{0}\nu\overline{\nu}$, which dominates for
$m_{N}\lesssim 1.3$ GeV, and $B^{0}\to\text{inv.}$, which takes over for
larger HNL masses. For the LNV operator (figure 7), $B^{0}\to\text{inv.}$
provides the most stringent constraint for any kinematically allowed $m_{N}$.
Finally, in the $b\to s$ scenario, the leading bound originates from the limit
on the branching ratio of $B^{+}\to K^{+}\nu\overline{\nu}$. All these bounds
have been derived by equating $\mathrm{BR}(P\to NN)$ ($\mathrm{BR}(P\to
P^{\prime}NN)$) obtained using Eq. (42) (Eq. (53)) to the corresponding upper
limit on $\mathrm{BR}(P\to\text{inv.})$ ($\mathrm{BR}(P\to
P^{\prime}\nu\overline{\nu})$) from table 5.

##  6 Summary and conclusions

In this work, we have focused on dimension-6 operators involving a pair of
heavy neutral leptons (HNLs), in the framework of the low-energy effective
field theory of the Standard Model extended with HNLs ($N_{R}$LEFT). We have
considered both lepton-number-conserving (LNC) and lepton-number-violating
(LNV) operators. Such operators with certain quark-flavor combinations can
lead to charm and bottom meson decays into a pair of HNLs, in either two-body
or three-body decays. We have considered three combinations of the quark
flavor indices, leading to the $c\to u$, $b\to d$, and $b\to s$ transitions,
and computed the corresponding production rates of the HNL pair from various
$D$\- and $B$-meson decays for (i) the LNC operators
$\mathcal{O}_{uN,12}^{V,RR}$, $\mathcal{O}_{dN,31}^{V,RR}$, and
$\mathcal{O}_{dN,32}^{V,RR}$, and (ii) the LNV operators
$\mathcal{O}_{uN,12}^{S,RR}$, $\mathcal{O}_{dN,31}^{S,RR}$, and
$\mathcal{O}_{dN,32}^{S,RR}$. The corresponding decay branching ratios for
both LNC and LNV cases are shown in figures 2 and 3, respectively. Switching
on the corresponding LR operators (one at a time) would lead to the same
results. It is interesting to observe that for the LNC case, the mesons’
three-body and two-body decays are dominant, respectively, for relatively
lower and higher masses in the kinematically allowed range of the HNL mass,
for all three quark-flavor scenarios. On the other hand, for the LNV case, the
two-body decays are the most important channel for the whole kinematically
relevant mass range.

Since these operators entail a pair of HNLs and in the absence of other
effective interactions, HNLs can decay only via an additional HNL mixing with
the active neutrinos in the weak current.111111Actually, an HNL can also decay
through a pair-$N_{R}$ operator in conjunction with mixing, but such decays
are suppressed compared to the decays through mixing alone. Therefore, we
study HNL production from the EFT operators, and HNL decays through active-
heavy neutrino mixing. Since we focus on GeV-scale HNLs, small neutrino mixing
parameters naturally lead these HNLs to become long-lived. Here, we choose
therefore to consider a series of far-detector proposals at the LHC. These
include AL3X, ANUBIS, CODEX-b, FACET, FASER and FASER2, MAPP1 and MAPP2, and
MATHUSLA. These experiments have been proposed (some even approved and
operated) mainly for searching for long-lived particles (LLPs) and would be
constructed with a distance of about $5$–$500$ meters away from different
interaction points (IPs). With shielding materials on site between a far
detector and its corresponding IP, such experiments are usually assumed to be
background-free for LLP searches.

We then perform Monte-Carlo simulation with the tool Pythia8, to estimate the
acceptance rates of these far detectors for HNLs pair-produced from $D$\- and
$B$-meson decays. For small mixing and mass, these HNLs are long-lived and
decay outside the LHC main detectors. With the total production rates of the
mesons at the LHC in knowledge, we can estimate the number of HNL decays in
the far detectors, using decay width formulae for (i) the HNL production from
meson decays in the $N_{R}$LEFT, and (ii) the HNL decays via their mixing with
the active neutrinos. For simplicity, we consider the HNLs to be mixed with
only either the electron neutrino or the muon neutrino. Further, for the
signal, we consider only the HNL decays into three active neutrinos to be
invisible; all the other decay channels are included as signature.

We show our numerical results in two types of planes: $|V_{lN}|^{2}$ vs.
$m_{N}$ in figures 4 and 5, and $c$ vs. $m_{N}$ in figures 6 and 7. For all
these results, we consider both LNC and LNV dimension-six operators, for all
three quark-flavor scenarios. In general, we find the constraints on the LNV
operators to be slightly stronger than those on the LNC ones, because the
meson decay branching ratios into HNLs in the LNV case are larger by about one
order of magnitude than those in the LNC case (for the same value of the
corresponding Wilson coefficients); see figures 2 and 3.

Among all the studied far detectors, we find MATHUSLA, ANUBIS, and AL3X show
the best sensitivities, while FASER is usually the weakest. Compared to
FASER2, MAPP1 achieves very similar results with 100 times less data, mainly
as a result of its much larger volume and closer distance to the IP.

In the $|V_{lN}|^{2}$ vs. $m_{N}$ plots, we find these experiments can probe
$|V_{lN}|^{2}$ as low as $10^{-15}$–$10^{-6}$, for GeV-scale HNLs. Then, for
two fixed benchmark values of $|V_{lN}|^{2}=10^{-7}$ and
$|V_{lN}|^{2}=10^{-10}$, we obtain the sensitivity reach in the Wilson
coefficients $c$ as functions of the HNL mass, for these far-detector
experiments. We further convert $c$ to the new-physics (NP) scale
$\Lambda_{\mathrm{NP}}$. We find for $|V_{lN}|^{2}=10^{-7}$, the considered
experiments can probe the NP scale up to between approximately $10$ and
several hundred TeV depending on the type of the operator considered, while
for $|V_{lN}|^{2}=10^{-10}$, the reach in $\Lambda_{\mathrm{NP}}$ is reduced
by slightly less than one order of magnitude, as a result of the larger decay
length and hence smaller acceptance.

## Acknowledgements

This work is supported by the Spanish grants PID2020-113775GB-I00
(AEI/10.13039/ 501100011033) and CIPROM/2021/054 (Generalitat Valenciana).
R.B. acknowledges financial support from the Generalitat Valenciana (grant
ACIF/2021/052). G.C. acknowledges support from ANID FONDECYT grant No.
11220237. G.C. and J.C.H. also acknowledge support from ANID FONDECYT grant
No. 1201673 and ANID – Millennium Science Initiative Program ICN2019_044. The
work of A.T. is supported by the “Generalitat Valenciana” under the grant
PROMETEO/2019/087 and by the MICINN-AEI (10.13039/501100011033) grant
PID2020-113334GB-I00. Z.S.W. is supported by the Ministry of Science and
Technology (MoST) of Taiwan with grant number MoST-110-2811-M-007-542-MY3.

##  Appendix A Pseudoscalar meson decays into HNLs

The four-fermion interactions described by the $N_{R}$LEFT operators in tables
1 and 2 lead to two-body and three-body meson decays, if the HNL mass is small
enough. The former involve two neutral leptons in the final state, whereas the
latter correspond to semileptonic decays with an additional lighter meson
within the final products. This appendix gathers the necessary formulae to
determine meson decay rates into HNLs, assuming a pseudoscalar meson (of mass
$m_{P}$ and momentum $p$) in the initial state.

In pure leptonic decays we define the meson decay constants through the
transition matrix elements (see e.g. Ref. DeVries:2020jbs ):

$\displaystyle\langle
0|\overline{q_{i}}\gamma^{\mu}\gamma_{5}q_{j}|P(p)\rangle$
$\displaystyle=if_{P}p^{\mu}\,,$ $\displaystyle\langle
0|\overline{q_{i}}\gamma_{5}q_{j}|P(p)\rangle$
$\displaystyle=i\frac{m_{P}^{2}}{m_{q_{i}}+m_{q_{j}}}f_{P}\equiv
if_{P}^{S}\,,$ (30)

where $i,j$ denote quark flavor indices, and $|P(p)\rangle$ is the decaying
pseudoscalar meson.

In semileptonic meson decays, the hadronic matrix elements are described by
the corresponding form factors. We distinguish two scenarios depending on the
nature of the outgoing meson. If it is a pseudoscalar (of mass $m$ and
momentum $p^{\prime}$), the non-vanishing matrix elements are:

$\displaystyle\langle
P^{\prime}(p^{\prime})|\overline{q_{i}}\gamma^{\mu}q_{j}|P(p)\rangle$
$\displaystyle=f_{+}(q^{2})\left[(p+p^{\prime})^{\mu}-\frac{m_{P}^{2}-m^{2}}{q^{2}}q^{\mu}\right]+f_{0}(q^{2})\frac{m_{P}^{2}-m^{2}}{q^{2}}q^{\mu}\,,$
$\displaystyle\langle
P^{\prime}(p^{\prime})|\overline{q_{i}}q_{j}|P(p)\rangle$
$\displaystyle=f_{S}(q^{2})\,,$ $\displaystyle\langle
P^{\prime}(p^{\prime})|\overline{q_{i}}\sigma^{\mu\nu}q_{j}|P(p)\rangle$
$\displaystyle=\frac{2i}{m_{P}+m}\left[p^{\mu}p^{\prime\nu}-p^{\nu}p^{\prime\mu}\right]f_{T}(q^{2})\,,$
(31)

where $q^{\mu}\equiv p^{\mu}-p^{\prime\mu}$. Of the four form factors
($f_{+}$, $f_{0}$, $f_{S}$ and $f_{T}$), the scalar one can be related to the
others using the equations of motion, through

$f_{S}(q^{2})=\frac{m_{P}^{2}-m^{2}}{m_{q_{j}}-m_{q_{i}}}f_{0}(q^{2})\,.$ (32)

When the final meson is a vector (of mass $m$, outgoing momentum $p^{\prime}$
and polarization vector $\epsilon$), the non-zero matrix elements are:

$\displaystyle\langle
V(p^{\prime},\epsilon)|\overline{q_{i}}\gamma^{\mu}q_{j}|P(p)\rangle$
$\displaystyle=ig(q^{2})\epsilon^{\mu\nu\alpha\beta}\epsilon_{\nu}^{\ast}P_{\alpha}q_{\beta}\,,$
$\displaystyle\langle
V(p^{\prime},\epsilon)|\overline{q_{i}}\gamma^{\mu}\gamma_{5}q_{j}|P(p)\rangle$
$\displaystyle=f(q^{2})\epsilon^{\ast\mu}+a_{+}(q^{2})P^{\mu}\epsilon^{\ast}\cdot
p+a_{-}(q^{2})q^{\mu}\epsilon^{\ast}\cdot p\,,$ $\displaystyle\langle
V(p^{\prime},\epsilon)|\overline{q_{i}}\sigma^{\mu\nu}q_{j}|P(p)\rangle$
$\displaystyle=\epsilon^{\mu\nu\alpha\beta}\left(g_{+}(q^{2})\epsilon_{\alpha}^{\ast}P_{\beta}+g_{-}(q^{2})\epsilon_{\alpha}^{\ast}q_{\beta}+g_{0}(q^{2})p_{\alpha}p^{\prime}_{\beta}\,\epsilon^{\ast}\cdot
p\right)\,,$ $\displaystyle\langle
V(p^{\prime},\epsilon)|\overline{q_{i}}\gamma_{5}q_{j}|P(p)\rangle$
$\displaystyle=f_{PS}(q^{2})\,\epsilon^{\ast}\cdot p\,,$ (33)

where we defined $P^{\mu}\equiv p^{\mu}+p^{\prime\mu}$. These transitions,
denoted by $P\rightarrow V$, are parametrized by eight functions of $q^{2}$:
$g$, $f$, $a_{+}$, $a_{-}$, $g_{+}$, $g_{-}$, $g_{0}$ and $f_{PS}$. Equations
of motion relate the pseudoscalar form factor with the rest through the
following expression:

$f_{PS}(q^{2})=\frac{1}{m_{q_{i}}+m_{q_{j}}}\left[f(q^{2})+a_{+}(q^{2})\left(m_{P}^{2}-m^{2}\right)+a_{-}(q^{2})q^{2}\right].$
(34)

###  A.1 Two-body decays

In the rest frame of a particle of mass $M$, the two-body decay width reads

$\Gamma=\frac{\lambda^{1/2}\left(M^{2},m_{a}^{2},m_{b}^{2}\right)}{64\pi^{2}M^{3}}\int\mathrm{d}\Omega_{\mathrm{cm}}\sum_{\mathrm{spins}}\left|\mathcal{M}\right|^{2}\,,$
(35)

where $\lambda\left(x,y,z\right)\equiv(x-y-z)^{2}-4yz$, and $m_{a}$, $m_{b}$
denote the masses of the decay products. Note that for identical particles in
the final state one has $\int\mathrm{d}\Omega_{\mathrm{cm}}=2\pi$.

In this section, we list the amplitudes and corresponding partial decay widths
of a neutral pseudoscalar meson $P\sim\overline{q_{i}}q_{j}$ decaying into (i)
a pair of HNLs, and (ii) one HNL plus one active neutrino. For a Dirac HNL
(and Dirac $\nu$), these decays occur via the LNC pair-$N_{R}$ operators in
table 1 and the LNC single-$N_{R}$ operators in table 2. In the Majorana
scenario (where the HNL and $\nu$ are Majorana), the decays receive
contributions from both the LNC and LNV operators. The amplitudes for Dirac
neutrinos read

$\displaystyle\mathcal{M}\left(P\to N\overline{N}\right)$
$\displaystyle=i\frac{f_{P}}{2}\left(c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right)\overline{u_{N}}\not{p}P_{R}v_{N^{\prime}}\,,$
(36) $\displaystyle\mathcal{M}\left(P\to\nu_{\alpha}\overline{N}\right)$
$\displaystyle=i\frac{f_{P}^{S}}{2}\left(c_{q\nu N,\,ij\alpha}^{S,RR}-c_{q\nu
N,\,ij\alpha}^{S,LR}\right)\overline{u_{\nu}}P_{R}v_{N}\,,$ (37)

where $u_{N}$ ($u_{\nu}$) and $v_{N}$ are the spinors corresponding to $N$
($\nu$) and $\overline{N}$, respectively. A prime is added when there are two
HNLs in the final state. For Majorana neutrinos, we recall that the full
Majorana HNL field is $N=N_{R}^{c}+N_{R}$, such that $N^{c}=N$. The amplitude
for $P\to NN$ has to be anti-symmetrized in $N$ and $N^{\prime}$, since there
are two identical particles in the final state.121212The convention for
associating a $u$ or $v$ spinor to a Majorana particle is arbitrary. Our is
such that there is one $\overline{u}$ spinor and one $v$ spinor in each
process. In this way, we recover the usual Dirac relations for the sums over
spins. Nonetheless, the physical result does not depend on the convention’s
choice.

$\displaystyle\mathcal{M}(P\to NN)$
$\displaystyle=i\frac{f_{P}}{2}\left(c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right)\overline{u_{N}}\not{p}\gamma_{5}v_{N^{\prime}}$
$\displaystyle+if_{P}^{S}\bigg{[}\left(c_{qN,\,ij}^{S,RR}-c_{qN,\,ij}^{S,LR}\right)\overline{u_{N}}P_{R}v_{N^{\prime}}-\left(c_{qN,\,ji}^{S,RR\,\ast}-c_{qN,\,ji}^{S,LR\,\ast}\right)\overline{u_{N}}P_{L}v_{N^{\prime}}\bigg{]}\,,$
(38) $\displaystyle\mathcal{M}\left(P\to\nu_{\alpha}N\right)$
$\displaystyle=i\frac{f_{P}^{S}}{2}\bigg{[}\left(c_{q\nu
N,\,ij\alpha}^{S,RR}-c_{q\nu
N,\,ij\alpha}^{S,LR}\right)\overline{u_{\nu}}P_{R}v_{N}-\left(c_{q\nu
N,\,ji\alpha}^{S,RR\,\ast}-c_{q\nu
N,\,ji\alpha}^{S,LR\,\ast}\right)\overline{u_{\nu}}P_{L}v_{N}\bigg{]}$
$\displaystyle+i\frac{f_{P}}{2}\bigg{[}\left(c_{q\nu
N,\,ij\alpha}^{V,RR}-c_{q\nu
N,\,ij\alpha}^{V,LR}\right)\overline{u_{\nu}}\not{p}P_{R}v_{N}-\left(c_{q\nu
N,\,ji\alpha}^{V,RR\,\ast}-c_{q\nu
N,\,ji\alpha}^{V,LR\,\ast}\right)\overline{u_{\nu}}\not{p}P_{L}v_{N}\bigg{]}\,.$
(39)

The decay rates of these two processes in the Dirac scenario are given by

$\displaystyle\Gamma\left(P\to N\overline{N}\right)$
$\displaystyle=\frac{\left|f_{P}\right|^{2}}{32\pi}\left|c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right|^{2}m_{P}\,m_{N}^{2}\,\sqrt{1-\frac{4m_{N}^{2}}{m_{P}^{2}}}\,,$
(40) $\displaystyle\Gamma\left(P\to\nu_{\alpha}\overline{N}\right)$
$\displaystyle=\frac{\left|f_{P}^{S}\right|^{2}}{64\pi}\left|c_{q\nu
N,\,ij\alpha}^{S,RR}-c_{q\nu
N,\,ij\alpha}^{S,LR}\right|^{2}m_{P}\left(1-\frac{m_{N}^{2}}{m_{P}^{2}}\right)^{2}\,,$
(41)

and in the Majorana case by

$\displaystyle\Gamma\left(P\to NN\right)$
$\displaystyle=\frac{m_{P}}{32\pi}\sqrt{1-\frac{4m_{N}^{2}}{m_{P}^{2}}}\Bigg{[}2\left|f_{P}\right|^{2}\left|c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right|^{2}m_{N}^{2}$
$\displaystyle+\left|f_{P}^{S}\right|^{2}\bigg{\\{}\left(\left|c_{qN,\,ij}^{S,RR}-c_{qN,\,ij}^{S,LR}\right|^{2}+\left|c_{qN,\,ji}^{S,RR}-c_{qN,\,ji}^{S,LR}\right|^{2}\right)\left(1-\frac{2m_{N}^{2}}{m_{P}^{2}}\right)$
$\displaystyle\hskip
42.67912pt+2\left[\left(c_{qN,\,ij}^{S,RR}-c_{qN,\,ij}^{S,LR}\right)\left(c_{qN,\,ji}^{S,RR}-c_{qN,\,ji}^{S,LR}\right)+\text{h.c.}\right]\frac{m_{N}^{2}}{m_{P}^{2}}\bigg{\\}}$
$\displaystyle+f_{P}f_{P}^{S}\bigg{\\{}\left(c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right)\left(c_{qN,\,ij}^{S,RR\,\ast}-c_{qN,\,ij}^{S,LR\,\ast}+c_{qN,\,ji}^{S,RR}-c_{qN,\,ji}^{S,LR}\right)m_{N}+\text{h.c.}\bigg{\\}}\Bigg{]}\,,$
(42) $\displaystyle\Gamma\left(P\to\nu_{\alpha}N\right)$
$\displaystyle=\frac{m_{P}}{64\pi}\left(1-\frac{m_{N}^{2}}{m_{P}^{2}}\right)^{2}\Bigg{[}\left|f_{P}^{S}\right|^{2}\left(\left|c_{q\nu
N,\,ij\alpha}^{S,RR}-c_{q\nu N,\,ij\alpha}^{S,LR}\right|^{2}+\left|c_{q\nu
N,\,ji\alpha}^{S,RR}-c_{q\nu N,\,ji\alpha}^{S,LR}\right|^{2}\right)$
$\displaystyle+\left|f_{P}\right|^{2}\left(\left|c_{q\nu
N,\,ij\alpha}^{V,RR}-c_{q\nu N,\,ij\alpha}^{V,LR}\right|^{2}+\left|c_{q\nu
N,\,ji\alpha}^{V,RR}-c_{q\nu N,\,ji\alpha}^{V,LR}\right|^{2}\right)m_{N}^{2}$
$\displaystyle+f_{P}f_{P}^{S}\bigg{\\{}\bigg{[}\left(c_{q\nu
N,\,ij\alpha}^{S,RR}-c_{q\nu N,\,ij\alpha}^{S,LR}\right)\left(c_{q\nu
N,\,ji\alpha}^{V,RR}-c_{q\nu N,\,ji\alpha}^{V,LR}\right)$
$\displaystyle\qquad\qquad+\left(c_{q\nu N,\,ji\alpha}^{S,RR\,\ast}-c_{q\nu
N,\,ji\alpha}^{S,LR\,\ast}\right)\left(c_{q\nu
N,\,ij\alpha}^{V,RR\,\ast}-c_{q\nu
N,\,ij\alpha}^{V,LR\,\ast}\right)\bigg{]}m_{N}+\text{h.c.}\bigg{\\}}\,\Bigg{]}\,.$
(43)

###  A.2 Three-body decays

For the three-body phase space integral calculations we follow the procedure
described in Ref. (DeVries:2020jbs, , Appendix B.1). However, we adopt general
labels for the leptonic products. In the rest frame of the decaying
pseudoscalar meson, the decay rate reads

$\displaystyle\Gamma$
$\displaystyle=\frac{1}{2m_{P}}\frac{1}{(2\pi)^{5}}\int\mathrm{d}^{4}p^{\prime}\int\mathrm{d}^{4}p_{a}\int\mathrm{d}^{4}p_{b}\,\delta(p^{\prime
2}-m^{2})\delta(p_{a}^{2}-m_{a}^{2})\delta(p_{b}^{2}-m_{b}^{2})$
$\displaystyle\times\sum_{\mathrm{spins}}\left|\mathcal{M}\right|^{2}\delta^{(4)}(p-p^{\prime}-p_{a}-p_{b})\,,$
(44)

where $p$ ($p^{\prime}$) is the momentum of the decaying (outgoing) meson and
$p_{a}$ and $p_{b}$ denote the momenta of the produced leptons.131313The
labels $a$ and $b$ will refer to (i) $N^{\prime}$ and $N$, and (ii) $\nu$ and
$N$. The summed over spins squared matrix element can be written in terms of
four-momentum invariant scalar products and the hadronic form factors, which
are functions of $q^{2}=(p-p^{\prime})^{2}$. After performing the integral
over the four-momentum $p_{b}$, it is convenient to introduce the variable $a$
via a factor of $1=\int\mathrm{d}a\,\delta(a-q^{2})$. Then, the spin-summed
matrix element squared can be written as

$\sum_{\mathrm{spins}}\left|\mathcal{M}\right|^{2}\bigg{\rvert}_{p_{b}=q-p_{a}}=\sum_{n=0}^{K}c_{n}(a)(p_{a}\cdot
p)^{n}\,,$ (45)

where $c_{n}(a)$ are functions of $a$, particle masses and the hadronic form
factors. For the decays under study we have $K\leq 2$. After performing the
remaining momentum integrals, one arrives at the final decay rate
expression141414When the outgoing leptons are two Majorana HNLs an extra
factor of $1/2$ appears in the final expression.

$\Gamma=\frac{1}{2m_{P}}\frac{1}{(2\pi)^{5}}\int^{(m_{P}-m)^{2}}_{(m_{a}+m_{b})^{2}}\mathrm{d}a\,I_{P}\sum^{K}_{n=0}c_{n}(a)I_{n}\,,$
(46)

where the integrals $I_{P}$ and $I_{n}$ (for $n={0,1,2}$) are given by

$\displaystyle I_{P}$
$\displaystyle=\frac{\pi}{2m_{P}^{2}}\lambda^{1/2}\left(a,m^{2},m_{P}^{2}\right)\,,$
$\displaystyle I_{0}$
$\displaystyle=\frac{\pi}{2a}\lambda^{1/2}\left(a,m_{a}^{2},m_{b}^{2}\right)\,,\qquad
I_{1}=\frac{(p_{a}\cdot q)(p\cdot q)}{a}I_{0}\,,$ $\displaystyle I_{2}$
$\displaystyle=\frac{I_{0}}{a^{2}}\left((p_{a}\cdot q)^{2}(p\cdot
q)^{2}+\frac{1}{48}\lambda(a,m^{2},m_{P}^{2})\lambda(a,m_{a}^{2},m_{b}^{2})\right)\,.$
(47)

The scalar products above have to be replaced by

$p_{a}\cdot q=\frac{1}{2}\left(a+m_{a}^{2}-m_{b}^{2}\right)\,,\qquad p\cdot
q=\frac{1}{2}\left(m_{P}^{2}+a-m^{2}\right)\,.$ (48)

In summary, after computing the spin-summed matrix element squared with
standard techniques and arranging it in the form of Eq. (45), one replaces
$(p_{a}\cdot p)^{n}$ by the corresponding integral result $I_{n}$, substitutes
the form factor expressions (see section A.3), and evaluates the integral in
Eq. (46) numerically.

In the rest of this section, we give the amplitudes of the relevant three-body
decays for producing the HNLs via the $N_{R}$LEFT operators. For Dirac HNLs,
these are

$\displaystyle\mathcal{M}\left(P\rightarrow P^{\prime}N\overline{N}\right)$
$\displaystyle=\frac{1}{2}\left(c_{qN,\,ij}^{V,RR}+c_{qN,\,ij}^{V,LR}\right)\left(f_{+}P_{\mu}+\left(f_{0}-f_{+}\right)\frac{m_{P}^{2}-m^{2}}{q^{2}}q_{\mu}\right)\left[\overline{u_{N}}\gamma^{\mu}P_{R}v_{N^{\prime}}\right]\,,$
(49) $\displaystyle\mathcal{M}\left(P\to VN\overline{N}\right)$
$\displaystyle=\bigg{\\{}\left(c_{qN,\,ij}^{V,RR}+c_{qN,\,ij}^{V,LR}\right)\left(ig\epsilon_{\mu\nu\alpha\beta}\epsilon^{\ast\nu}p^{\prime\alpha}p^{\beta}\right)$
$\displaystyle+\frac{1}{2}\left(c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right)\left(f\epsilon^{\ast}_{\mu}+a_{+}P_{\mu}\epsilon^{\ast}\cdot
p+a_{-}q_{\mu}\epsilon^{\ast}\cdot
p\right)\bigg{\\}}\left[\overline{u_{N}}\gamma^{\mu}P_{R}v_{N^{\prime}}\right]\,,$
(50) $\displaystyle\mathcal{M}\left(P\to
P^{\prime}\nu_{\alpha}\overline{N}\right)$
$\displaystyle=\frac{f_{S}}{2}\left(c_{q\nu N,\,ij\alpha}^{S,RR}+c_{q\nu
N,\,ij\alpha}^{S,LR}\right)\left[\overline{u_{\nu}}P_{R}v_{N}\right]$
$\displaystyle+c_{q\nu
N,\,ij\alpha}^{T,RR}\frac{2if_{T}}{m_{P}+m}\left(p_{\mu}p^{\prime}_{\nu}+\frac{i}{2}\epsilon_{\mu\nu\rho\sigma}p^{\rho}p^{\prime\sigma}\right)\,\left[\overline{u_{\nu}}\sigma^{\mu\nu}P_{R}v_{N}\right]\,,$
(51) $\displaystyle\mathcal{M}\left(P\to V\nu_{\alpha}\overline{N}\right)$
$\displaystyle=\frac{f_{PS}}{2}\left(c_{q\nu N,\,ij\alpha}^{S,RR}-c_{q\nu
N,\,ij\alpha}^{S,LR}\right)(\epsilon^{\ast}\cdot
p)\left[\overline{u_{\nu}}P_{R}v_{N}\right]$ $\displaystyle+c_{q\nu
N,\,ij\alpha}^{T,RR}\bigg{\\{}\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}\left[g_{+}\epsilon^{\ast\alpha}(p+p^{\prime})^{\beta}+g_{-}\epsilon^{\ast\alpha}q^{\beta}+g_{0}p^{\alpha}p^{\prime\beta}\epsilon^{\ast}\cdot
p\right]$
$\displaystyle-i\left[g_{+}\epsilon^{\ast}_{\mu}(p+p^{\prime})_{\nu}+g_{-}\epsilon^{\ast}_{\mu}q_{\nu}+g_{0}p_{\mu}p^{\prime}_{\nu}\epsilon^{\ast}\cdot
p\right]\bigg{\\}}\left[\overline{u_{\nu}}\sigma^{\mu\nu}P_{R}v_{N}\right]\,.$
(52)

The corresponding amplitudes for Majorana HNLs read

$\displaystyle\mathcal{M}(P\rightarrow P^{\prime}NN)$
$\displaystyle=\frac{1}{2}\left(c_{qN,\,ij}^{V,RR}+c_{qN,\,ij}^{V,LR}\right)\left(f_{+}P_{\mu}+\left(f_{0}-f_{+}\right)\frac{m_{P}^{2}-m^{2}}{q^{2}}q_{\mu}\right)\left[\overline{u_{N}}\gamma^{\mu}\gamma^{5}v_{N^{\prime}}\right]$
$\displaystyle+f_{S}\left(c_{qN,\,ij}^{S,RR}+c_{qN,\,ij}^{S,LR}\right)\left[\overline{u_{N}}P_{R}v_{N^{\prime}}\right]+f_{S}\left(c_{qN,\,ji}^{S,RR\ast}+c_{qN,\,ji}^{S,LR\ast}\right)\left[\overline{u_{N}}P_{L}v_{N^{\prime}}\right]\,,$
(53) $\displaystyle\mathcal{M}(P\rightarrow VNN)$
$\displaystyle=\bigg{\\{}\left(c_{qN,\,ij}^{V,RR}+c_{qN,\,ij}^{V,LR}\right)\left(ig\epsilon_{\mu\nu\alpha\beta}\epsilon^{\ast\nu}p^{\prime\alpha}p^{\beta}\right)$
$\displaystyle+\frac{1}{2}\left(c_{qN,\,ij}^{V,RR}-c_{qN,\,ij}^{V,LR}\right)\left(f\epsilon^{\ast}_{\mu}+a_{+}P_{\mu}\epsilon^{\ast}\cdot
p+a_{-}q_{\mu}\epsilon^{\ast}\cdot
p\right)\bigg{\\}}\left[\overline{u_{N}}\gamma^{\mu}\gamma^{5}v_{N^{\prime}}\right]$
$\displaystyle+f_{PS}\bigg{\\{}\left(c_{q\nu N,\,ij}^{S,RR}-c_{q\nu
N,\,ij}^{S,LR}\right)(\epsilon^{\ast}\cdot
p)\left[\overline{u_{N}}P_{R}v_{N^{\prime}}\right]$ $\displaystyle\hskip
42.67912pt-\left(c_{q\nu N,\,ji}^{S,RR\,\ast}-c_{q\nu
N,\,ji}^{S,LR\,\ast}\right)(\epsilon^{\ast}\cdot
p)\left[\overline{u_{N}}P_{L}v_{N^{\prime}}\right]\bigg{\\}}\,,$ (54)
$\displaystyle\mathcal{M}(P\rightarrow P^{\prime}\nu_{\alpha}N)$
$\displaystyle=\frac{f_{S}}{2}\bigg{\\{}\left(c_{q\nu
N,\,ij\alpha}^{S,RR}+c_{q\nu
N,\,ij\alpha}^{S,LR}\right)\left[\overline{u_{\nu}}P_{R}v_{N}\right]+\left(c_{q\nu
N,\,ji\alpha}^{S,RR\,\ast}+c_{q\nu
N,\,ji\alpha}^{S,LR\,\ast}\right)\left[\overline{u_{\nu}}P_{L}v_{N}\right]\bigg{\\}}$
$\displaystyle+\frac{if_{T}}{M+m}\bigg{\\{}c_{q\nu
N,\,ij\alpha}^{T,RR}\left(2p_{\mu}p^{\prime}_{\nu}+i\epsilon_{\mu\nu\rho\sigma}p^{\rho}p^{\prime\sigma}\right)\left[\overline{u_{\nu}}\sigma^{\mu\nu}P_{R}v_{N}\right]$
$\displaystyle\hskip 56.9055pt-c_{q\nu
N,\,ji\alpha}^{T,RR\,\ast}\left(2p_{\mu}p^{\prime}_{\nu}-i\epsilon_{\mu\nu\rho\sigma}p^{\rho}p^{\prime\sigma}\right)\left[\overline{u_{\nu}}\sigma^{\mu\nu}P_{L}v_{N}\right]\bigg{\\}}$
$\displaystyle+\bigg{\\{}\left(c_{q\nu N,\,ij\alpha}^{V,RR}+c_{q\nu
N,\,ij\alpha}^{V,LR}\right)\left[\overline{u_{\nu}}\gamma^{\mu}P_{R}v_{N}\right]-\left(c_{q\nu
N,\,ji\alpha}^{V,RR\,\ast}+c_{q\nu
N,\,ji\alpha}^{V,LR\,\ast}\right)\left[\overline{u_{\nu}}\gamma^{\mu}P_{L}v_{N}\right]\bigg{\\}}$
$\displaystyle\times\frac{1}{2}\left(f_{+}P_{\mu}+\left(f_{0}-f_{+}\right)\frac{m_{P}^{2}-m^{2}}{q^{2}}q_{\mu}\right)\,,$
(55) $\displaystyle\mathcal{M}\left(P\rightarrow V\nu_{\alpha}N\right)$
$\displaystyle=\frac{f_{PS}}{2}\bigg{\\{}\left(c_{q\nu
N,\,ij\alpha}^{S,RR}-c_{q\nu N,\,ij\alpha}^{S,LR}\right)(\epsilon^{\ast}\cdot
p)\left[\overline{u_{\nu}}P_{R}v_{N}\right]$ $\displaystyle\hskip
42.67912pt-\left(c_{q\nu N,\,ji\alpha}^{S,RR\,\ast}-c_{q\nu
N,\,ji\alpha}^{S,LR\,\ast}\right)(\epsilon^{\ast}\cdot
p)\left[\overline{u_{\nu}}P_{L}v_{N}\right]\bigg{\\}}$ $\displaystyle+c_{q\nu
N,\,ij\alpha}^{T,RR}\left(\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}T^{\alpha\beta}-iT_{\mu\nu}\right)\left[\overline{u_{\nu}}\sigma^{\mu\nu}P_{R}v_{N}\right]$
$\displaystyle-c_{q\nu
N,\,ji\alpha}^{T,RR\,\ast}\left(\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}T^{\alpha\beta}+iT_{\mu\nu}\right)\left[\overline{u_{\nu}}\sigma^{\mu\nu}P_{L}v_{N}\right]$
$\displaystyle+\frac{1}{2}\bigg{\\{}\left(c_{q\nu N,\,ij\alpha}^{V,RR}+c_{q\nu
N,\,ij\alpha}^{V,LR}\right)V_{\mu}+\left(c_{q\nu N,\,ij\alpha}^{V,RR}-c_{q\nu
N,\,ij\alpha}^{V,LR}\right)A_{\mu}\bigg{\\}}\left[\overline{u_{\nu}}\gamma^{\mu}P_{R}v_{N}\right]$
$\displaystyle-\frac{1}{2}\bigg{\\{}\left(c_{q\nu
N,\,ji\alpha}^{V,RR\,\ast}+c_{q\nu
N,\,ji\alpha}^{V,LR\,\ast}\right)V_{\mu}+\left(c_{q\nu
N,\,ji\alpha}^{V,RR\,\ast}-c_{q\nu
N,\,ji\alpha}^{V,LR\,\ast}\right)A_{\mu}\bigg{\\}}\left[\overline{u_{\nu}}\gamma^{\mu}P_{L}v_{N}\right]\,,$
(56)

where we have introduced:

$\displaystyle T^{\alpha\beta}$ $\displaystyle\equiv
g_{+}\epsilon^{\ast\,\alpha}(p+p^{\prime})^{\beta}+g_{-}\epsilon^{\ast\,\alpha}q^{\beta}+g_{0}p^{\alpha}p^{\prime\beta}(\epsilon^{\ast}\cdot
p)\,,$ $\displaystyle V_{\mu}$ $\displaystyle\equiv
2ig\epsilon_{\mu\nu\alpha\beta}\epsilon^{\ast\nu}p^{\prime\alpha}p^{\beta}\,,$
$\displaystyle A_{\mu}$ $\displaystyle\equiv
f\epsilon^{\ast}_{\mu}+a_{+}(p+p^{\prime})_{\mu}(\epsilon^{\ast}\cdot
p)+a_{-}q_{\mu}(\epsilon^{\ast}\cdot p)\,.$

Specific software, such as `FeynCalc` Shtabovenko:2016sxi ;
Shtabovenko:2020gxv , automatizes the computation of the summed over spins
squared matrix elements. Finally, to simplify notation in the previous
expressions, we have omitted the $q^{2}$ dependence of the form factors.

###  A.3 Decay constants and form factors

The masses of all mesons involved in the decays are taken from Ref.
Workman:2022ynf . In addition, we use the light quark masses at a
renormalization scale of $\mu=2$ GeV in $\overline{\text{MS}}$, whereas for
$m_{c}$ and $m_{b}$ we employ the $\overline{\text{MS}}$ masses at $\mu=m_{c}$
and $\mu=m_{b}$, respectively, as in Ref. Workman:2022ynf :

$m_{u}=2.2\text{ MeV},\quad m_{d}=4.7\text{ MeV},\quad m_{s}=93\text{
MeV},\quad m_{c}=1.27\text{ GeV},\quad m_{b}=4.18\text{ GeV}.$

We provide the decay constants of the heavy neutral pseudoscalar mesons in
table 7. We use the isospin-averaged results151515In other words, we assume
the isospin symmetry, i.e. the (approximate) symmetry between the $u$ and $d$
quarks. In this limit, the decay constants for $D^{0}$ and $D^{+}$ are the
same, $f_{D}$. Similarly, $B^{0}$ and $B^{+}$ mesons are characterised by one
constant, $f_{B}$. The difference $|f_{P^{+}}-f_{P}|$ is estimated to be
$\simeq 0.5$ MeV for both $P=D$ and $P=B$ Aoki:2021kgd ; PDG2021update . for
$N_{f}=2+1+1$ dynamical quark flavors from Ref. Aoki:2021kgd .

Meson $P$ | Decay constant $f_{P}$ [MeV]
---|---
$D^{0}$ | 212.0 (0.7)
$B^{0}$ | 190.0 (1.3)
$B^{0}_{s}$ | 230.3 (1.3)
Table 7: Decay constants of heavy neutral pseudoscalar mesons Aoki:2021kgd .

As regards hadronic form factors, several parametrizations exist in the
literature for describing their dependence on $q^{2}$. One of the most used is
the BCL method Bourrely:2008za , and we try to stick to it as long as there is
available data. In all other cases, we use the so-called double-pole
parametrization. In table 8, we show all hadronic transitions mediated by the
selected $N_{R}$LEFT operators and the references studying the associated form
factors. We indicate the equation where the explicit $q^{2}$ parametrization
appears and the tables where the best-fit parameter values are
given.161616Notice that there are several conventions in the literature for
writing hadronic matrix elements in terms of form factors. The basis presented
in Eqs. (31) and (33) must be converted to the new basis in each reference
before using their $q^{2}$ parametrization.

Meson transitions | Form factor parametrization | Fitted parameters
---|---|---
$D$-mesons | $D^{0}\rightarrow\pi^{0}$ | BCL (Lubicz:2017syv, , Eqs. (A1)–(A2)) | (Lubicz:2017syv, , Tab. VI)
$D^{+}\rightarrow\pi^{+}$
$D^{0}\rightarrow\eta,\,\eta^{\prime},\,\rho^{0},\,\omega$ | Double-pole (Ivanov:2019nqd, , Eq. (63)) | (Ivanov:2019nqd, , Tab. VI)
$D^{+}\rightarrow\rho^{+}$
$D_{s}^{+}\rightarrow K^{+},\,K^{\ast+}$
$B$-mesons | $B^{0}\rightarrow\pi^{0},\,K^{0}$ | BCL (Aoki:2021kgd, , Eqs. (529)–(530)) | (Aoki:2021kgd, , Tabs. 46, 48, 50, 51)
$B^{+}\rightarrow\pi^{+},\,K^{+}$
$B_{s}^{0}\rightarrow\overline{K^{0}}$
$B^{0}\rightarrow\eta,\,\eta^{\prime}$ | Double-pole (Wu:2006rd, , Eq. (36)) | (Wu:2006rd, , Tab. 4)
$B_{s}^{0}\rightarrow\eta,\,\eta^{\prime}$
$B^{0}\rightarrow\rho^{0},\,\omega,\,K^{\ast 0}$ | BCL (Bharucha:2015bzk, , Eq. (2.16)) | (Bharucha:2015bzk, , Tab. 14)
$B^{+}\rightarrow\rho^{+},\,K^{\ast+}$
$B_{s}^{0}\rightarrow\phi,\,\overline{K^{\ast 0}}$
Table 8: Relevant meson transitions for the decays studied in section 3. The
two last columns show the references to the form factor parametrizations and
the best-fit parameter values. Isospin symmetry is assumed so that identical
form factors are taken for $P^{0}\rightarrow M^{0}$ and $P^{+}\rightarrow
M^{+}$ transitions, where $M$ can be a pseudoscalar or vector meson.

##  Appendix B Branching ratios of meson decays triggered by single-$N_{R}$
operators

We have also computed the branching ratios for two- and three-body
pseudoscalar meson decays including in the final state a SM neutrino and an
HNL. These decays are triggered by the single-$N_{R}$ operators given in table
2. As an example, we show in figure 8 the corresponding branching ratios in
the case of $b\to d$ transitions triggered by three single-$N_{R}$ operators,
namely, $\mathcal{O}_{d\nu N,31\alpha}^{S,RR}$, $\mathcal{O}_{d\nu
N,31\alpha}^{T,RR}$, and $\mathcal{O}_{d\nu N,31\alpha}^{V,RR}$. The scalar
and tensor operators are LNC, whereas the operator of the vector type is LNV.

Figure 8: Branching ratios of $B$-meson decays triggered by the single-$N_{R}$
operators $\mathcal{O}_{d\nu N,31\alpha}^{S,RR}$, $\mathcal{O}_{d\nu
N,31\alpha}^{T,RR}$ and $\mathcal{O}_{d\nu N,32\alpha}^{V,RR}$. In each case,
the corresponding Wilson coefficient has been set to $10^{-3}v^{-2}$.

As can be seen from Eq. (43), the tensor operator does not contribute to the
two-body decay width. This is why there is no solid blue line in the
corresponding plot in figure 8. The tensor form factors for the $B\to\pi$
transitions have been taken from table 50 in Ref. Aoki:2021kgd . For
$B^{0}\to\eta,~{}\eta^{\prime}$ and $B_{s}^{0}\to\overline{K^{0}}$, we have
used Eq. (27) from Ref. Wu:2006rd , which relates $f_{T}$ to $f_{+}$ and
$f_{0}$. Finally, for the $B_{(s)}\to V$ transitions, we have used the results
of Ref. Bharucha:2015bzk . In table 8 we include these references and the ones
for all other form factors.

We note that the branching ratios in figure 8 are very similar to those for
the meson decays induced by the single-$N_{R}$ operators with a charged
lepton, which have been studied in detail in Ref. DeVries:2020jbs , cf. figure
2 therein. Thus, the constraints on the neutral current single-$N_{R}$
operators from the far detectors are expected to be very similar to those
derived in Ref. DeVries:2020jbs  for the charged current single-$N_{R}$
operators. For this reason, we do not study the neutral current single-$N_{R}$
operators in more detail.

##  Appendix C HNL two-body decays via single-$N_{R}$ operators

Single-$N_{R}$ operators contribute to both HNL production and decay. In fact,
depending on the Wilson coefficient and the value of active-heavy mixing, they
can even dominate the HNL decay width. The LNC and LNV operators in table 2
make HNLs undergo two- and three-body semileptonic decays. The dominant
channels are two-body decays, mediated by just one insertion of an effective
operator, with one light neutral meson and one active neutrino in the final
state. If the HNL mass is large enough, then multimeson final states could
also be important in the computation of the total decay width. We do not
compute this contribution, though. The single-$N_{R}$ operators with light
quarks can mediate the two-body decay processes mentioned above at tree level.
The contribution from operators with heavy quarks, which control $B$\- and
$D$-meson decays into HNLs, is loop-suppressed, and we do not compute it. The
interference between active-heavy mixing and single-$N_{R}$ operators is also
not treated here.

In this appendix, we list the HNL partial decay rates into a neutral meson and
a SM neutrino mediated by the effective operators containing two light quarks.
We leave, however, the flavor indices of the Wilson coefficients unspecified.
We consider the scenarios where both $N$ and $\nu$ are either Dirac or
Majorana particles. In each case, we distinguish between a pseudoscalar and a
vector meson in the final state. The latter requires the introduction of the
following form factors:

$\displaystyle\langle
0|\overline{q_{i}}\gamma^{\mu}q_{j}|V(p,\epsilon)\rangle$
$\displaystyle=if_{V}m_{V}\epsilon^{\mu}\,,$ $\displaystyle\langle
0|\overline{q_{i}}\sigma^{\mu\nu}q_{j}|V(p,\epsilon)\rangle$
$\displaystyle=-f_{V}^{T}(p^{\mu}\epsilon^{\nu}-p^{\nu}\epsilon^{\mu})\,,$
(57)

where $|V(p,\epsilon)\rangle$ denotes a vector meson with mass $m_{V}$,
momentum $p$, and polarization vector $\epsilon$. For Dirac $N$ and $\nu$, the
decay rate into a SM neutrino and a pseudoscalar (vector) meson
$P~{}(V)\sim\overline{q_{i}}q_{j}$ of mass $m_{P}$ ($m_{V}$) reads

$\displaystyle\Gamma(N\rightarrow\nu_{\alpha}P)$
$\displaystyle=\frac{\left|f_{P}^{S}\right|^{2}}{128\pi}\left|c_{q\nu
N,\,ji\alpha}^{S,RR}-c_{q\nu
N,\,ji\alpha}^{S,LR}\right|^{2}m_{N}\left(1-\frac{m_{P}^{2}}{m_{N}^{2}}\right)^{2}\,,$
(58) $\displaystyle\Gamma(N\rightarrow\nu_{\alpha}V)$
$\displaystyle=\frac{\left|f_{V}^{T}\right|^{2}}{4\pi}\left|c_{q\nu
N,\,ji\alpha}^{T,RR}\right|^{2}m_{N}^{3}\left(1-\frac{m_{V}^{2}}{m_{N}^{2}}\right)^{2}\left(1+\frac{m_{V}^{2}}{2m_{N}^{2}}\right)\,.$
(59)

In the Majorana scenario, the decay rates receive contributions from both the
LNC and LNV single-$N_{R}$ operators and are given by

$\displaystyle\Gamma(N\rightarrow\nu_{\alpha}P)$
$\displaystyle=\frac{m_{N}}{128\pi}\left(1-\frac{m_{P}^{2}}{m_{N}^{2}}\right)^{2}\Bigg{\\{}\left|f_{P}^{S}\right|^{2}\left(\left|c_{q\nu
N,\,ji\alpha}^{S,RR}-c_{q\nu N,\,ji\alpha}^{S,LR}\right|^{2}+\left|c_{q\nu
N,\,ij\alpha}^{S,RR}-c_{q\nu N,\,ij\alpha}^{S,LR}\right|^{2}\right)$
$\displaystyle+\left|f_{P}\right|^{2}\left(\left|c_{q\nu
N,\,ji\alpha}^{V,RR}-c_{q\nu N,\,ji\alpha}^{V,LR}\right|^{2}+\left|c_{q\nu
N,\,ij\alpha}^{V,RR}-c_{q\nu N,\,ij\alpha}^{V,LR}\right|^{2}\right)m_{N}^{2}$
$\displaystyle-f_{P}f_{P}^{S}\bigg{[}\left(c_{q\nu
N,\,ji\alpha}^{S,RR}-c_{q\nu N,\,ji\alpha}^{S,LR}\right)\left(c_{q\nu
N,\,ij\alpha}^{V,RR}-c_{q\nu N,\,ij\alpha}^{V,LR}\right)$ $\displaystyle\hskip
39.83368pt+\left(c_{q\nu N,\,ij\alpha}^{S,RR}-c_{q\nu
N,\,ij\alpha}^{S,LR}\right)\left(c_{q\nu N,\,ji\alpha}^{V,RR}-c_{q\nu
N,\,ji\alpha}^{V,LR}\right)+\text{h.c.}\bigg{]}m_{N}\Bigg{\\}}\,,$ (60)
$\displaystyle\Gamma(N\rightarrow\nu_{\alpha}V)$
$\displaystyle=\frac{m_{N}^{3}}{128\pi}\left(1-\frac{m_{V}^{2}}{m_{N}^{2}}\right)^{2}\Bigg{\\{}32\left|f_{V}^{T}\right|^{2}\left(\left|c_{q\nu
N,\,ij\alpha}^{T,RR}\right|^{2}+\left|c_{q\nu
N,\,ji\alpha}^{T,RR}\right|^{2}\right)\left(1+\frac{m_{V}^{2}}{2m_{N}^{2}}\right)$
$\displaystyle+\left|f_{V}\right|^{2}\left(\left|c_{q\nu
N,\,ji\alpha}^{V,RR}+c_{q\nu N,\,ji\alpha}^{V,LR}\right|^{2}+\left|c_{q\nu
N,\,ij\alpha}^{V,RR}+c_{q\nu
N,\,ij\alpha}^{V,LR}\right|^{2}\right)\left(1+\frac{2m_{V}^{2}}{m_{N}^{2}}\right)$
$\displaystyle-12f_{V}^{T}f_{V}\bigg{[}c_{q\nu
N,\,ji\alpha}^{T,RR}\left(c_{q\nu N,\,ij\alpha}^{V,RR}+c_{q\nu
N,\,ij\alpha}^{V,LR}\right)$ $\displaystyle\hskip 51.21504pt+c_{q\nu
N,\,ij\alpha}^{T,RR}\left(c_{q\nu N,\,ji\alpha}^{V,RR}+c_{q\nu
N,\,ji\alpha}^{V,LR}\right)+\text{h.c.}\bigg{]}\frac{m_{V}}{m_{N}}\Bigg{\\}}\,.$
(61)

## References

  * (1) J. Alimena et al., _Searching for long-lived particles beyond the Standard Model at the Large Hadron Collider_ , _J. Phys. G_ 47 (2020) 090501 [1903.04497].
  * (2) D. Curtin et al., _Long-Lived Particles at the Energy Frontier: The MATHUSLA Physics Case_ , _Rept. Prog. Phys._ 82 (2019) 116201 [1806.07396].
  * (3) M. Bauer, O. Brandt, L. Lee and C. Ohm, _ANUBIS: Proposal to search for long-lived neutral particles in CERN service shafts_ , 1909.13022.
  * (4) V.V. Gligorov, S. Knapen, B. Nachman, M. Papucci and D.J. Robinson, _Leveraging the ALICE/L3 cavern for long-lived particle searches_ , _Phys. Rev. D_ 99 (2019) 015023 [1810.03636].
  * (5) V.V. Gligorov, S. Knapen, M. Papucci and D.J. Robinson, _Searching for Long-lived Particles: A Compact Detector for Exotics at LHCb_ , _Phys. Rev. D_ 97 (2018) 015023 [1708.09395].
  * (6) S. Cerci et al., _FACET: A new long-lived particle detector in the very forward region of the CMS experiment_ , _JHEP_ 2022 (2022) 110 [2201.00019].
  * (7) J.L. Feng, I. Galon, F. Kling and S. Trojanowski, _ForwArd Search ExpeRiment at the LHC_ , _Phys. Rev. D_ 97 (2018) 035001 [1708.09389].
  * (8) FASER collaboration, _FASER’s physics reach for long-lived particles_ , _Phys. Rev. D_ 99 (2019) 095011 [1811.12522].
  * (9) J.L. Pinfold, _The MoEDAL Experiment at the LHC—A Progress Report_ , _Universe_ 5 (2019) 47.
  * (10) J.L. Pinfold, _The MoEDAL experiment: a new light on the high-energy frontier_ , _Phil. Trans. Roy. Soc. Lond. A_ 377 (2019) 20190382.
  * (11) J.P. Chou, D. Curtin and H.J. Lubatti, _New Detectors to Explore the Lifetime Frontier_ , _Phys. Lett. B_ 767 (2017) 29 [1606.06298].
  * (12) MATHUSLA collaboration, _An Update to the Letter of Intent for MATHUSLA: Search for Long-Lived Particles at the HL-LHC_ , 2009.01693.
  * (13) F. Kling and S. Trojanowski, _Heavy Neutral Leptons at FASER_ , _Phys. Rev. D_ 97 (2018) 095016 [1801.08947].
  * (14) J.C. Helo, M. Hirsch and Z.S. Wang, _Heavy neutral fermions at the high-luminosity LHC_ , _JHEP_ 07 (2018) 056 [1803.02212].
  * (15) D. Dercks, H.K. Dreiner, M. Hirsch and Z.S. Wang, _Long-Lived Fermions at AL3X_ , _Phys. Rev. D_ 99 (2019) 055020 [1811.01995].
  * (16) M. Hirsch and Z.S. Wang, _Heavy neutral leptons at ANUBIS_ , _Phys. Rev. D_ 101 (2020) 055034 [2001.04750].
  * (17) J. De Vries, H.K. Dreiner, J.Y. Günther, Z.S. Wang and G. Zhou, _Long-lived Sterile Neutrinos at the LHC in Effective Field Theory_ , _JHEP_ 03 (2021) 148 [2010.07305].
  * (18) M. Ovchynnikov, V. Kryshtal and K. Bondarenko, _Sensitivity of the FACET experiment to Heavy Neutral Leptons and Dark Scalars_ , 2209.14870.
  * (19) G. Cottin, J.C. Helo and M. Hirsch, _Searches for light sterile neutrinos with multitrack displaced vertices_ , _Phys. Rev. D_ 97 (2018) 055025 [1801.02734].
  * (20) G. Cottin, J.C. Helo and M. Hirsch, _Displaced vertices as probes of sterile neutrino mixing at the LHC_ , _Phys. Rev. D_ 98 (2018) 035012 [1806.05191].
  * (21) K. Cheung, Y.-L. Chung, H. Ishida and C.-T. Lu, _Sensitivity reach on heavy neutral leptons and $\tau$-neutrino mixing $|U_{\tau N}|^{2}$ at the HL-LHC_, _Phys. Rev. D_ 102 (2020) 075038 [2004.11537].
  * (22) J. Jones-Pérez, J. Masias and J.D. Ruiz-Álvarez, _Search for Long-Lived Heavy Neutrinos at the LHC with a VBF Trigger_ , _Eur. Phys. J. C_ 80 (2020) 642 [1912.08206].
  * (23) J. Liu, Z. Liu, L.-T. Wang and X.-P. Wang, _Seeking for sterile neutrinos with displaced leptons at the LHC_ , _JHEP_ 07 (2019) 159 [1904.01020].
  * (24) ATLAS collaboration, _Search for heavy neutral leptons in decays of $W$ bosons produced in 13 TeV $pp$ collisions using prompt and displaced signatures with the ATLAS detector_, _JHEP_ 10 (2019) 265 [1905.09787].
  * (25) CMS collaboration, _Search for long-lived heavy neutral leptons with displaced vertices in proton-proton collisions at $\sqrt{s}$ =13 TeV_, 2201.05578.
  * (26) C.-W. Chiang, G. Cottin, A. Das and S. Mandal, _Displaced heavy neutrinos from $Z^{\prime}$ decays at the LHC_, _JHEP_ 12 (2019) 070 [1908.09838].
  * (27) F. Deppisch, S. Kulkarni and W. Liu, _Heavy neutrino production via $Z^{\prime}$ at the lifetime frontier_, _Phys. Rev. D_ 100 (2019) 035005 [1905.11889].
  * (28) I. Brivio and M. Trott, _The Standard Model as an Effective Field Theory_ , _Phys. Rept._ 793 (2019) 1 [1706.08945].
  * (29) F. del Aguila, S. Bar-Shalom, A. Soni and J. Wudka, _Heavy Majorana Neutrinos in the Effective Lagrangian Description: Application to Hadron Colliders_ , _Phys. Lett. B_ 670 (2009) 399 [0806.0876].
  * (30) A. Aparici, K. Kim, A. Santamaria and J. Wudka, _Right-handed neutrino magnetic moments_ , _Phys. Rev. D_ 80 (2009) 013010 [0904.3244].
  * (31) Y. Liao and X.-D. Ma, _Operators up to Dimension Seven in Standard Model Effective Field Theory Extended with Sterile Neutrinos_ , _Phys. Rev. D_ 96 (2017) 015012 [1612.04527].
  * (32) H.-L. Li, Z. Ren, M.-L. Xiao, J.-H. Yu and Y.-H. Zheng, _Operator bases in effective field theories with sterile neutrinos: d $\leq$ 9_, _JHEP_ 11 (2021) 003 [2105.09329].
  * (33) G. Cottin, J.C. Helo, M. Hirsch, A. Titov and Z.S. Wang, _Heavy neutral leptons in effective field theory and the high-luminosity LHC_ , _JHEP_ 09 (2021) 039 [2105.13851].
  * (34) R. Beltrán, G. Cottin, J.C. Helo, M. Hirsch, A. Titov and Z.S. Wang, _Long-lived heavy neutral leptons at the LHC: four-fermion single-N R operators_, _JHEP_ 01 (2022) 044 [2110.15096].
  * (35) 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].
  * (36) I. Bischer and W. Rodejohann, _General neutrino interactions from an effective field theory perspective_ , _Nucl. Phys. B_ 947 (2019) 114746 [1905.08699].
  * (37) M. Chala and A. Titov, _One-loop matching in the SMEFT extended with a sterile neutrino_ , _JHEP_ 05 (2020) 139 [2001.07732].
  * (38) T. Li, X.-D. Ma and M.A. Schmidt, _General neutrino interactions with sterile neutrinos in light of coherent neutrino-nucleus scattering and meson invisible decays_ , _JHEP_ 07 (2020) 152 [2005.01543].
  * (39) T. Li, X.-D. Ma and M.A. Schmidt, _Constraints on the charged currents in general neutrino interactions with sterile neutrinos_ , _JHEP_ 10 (2020) 115 [2007.15408].
  * (40) G. Zhou, J.Y. Günther, Z.S. Wang, J. de Vries and H.K. Dreiner, _Long-lived sterile neutrinos at Belle II in effective field theory_ , _JHEP_ 04 (2022) 057 [2111.04403].
  * (41) T. Han, J. Liao, H. Liu and D. Marfatia, _Right-handed Dirac and Majorana neutrinos at Belle II_ , 2207.07029.
  * (42) D. Barducci, E. Bertuzzo, M. Taoso and C. Toni, _Probing right-handed neutrinos dipole operators_ , 2209.13469.
  * (43) R.M. Fonseca, _The Sym2Int program: going from symmetries to interactions_ , _J. Phys. Conf. Ser._ 873 (2017) 012045 [1703.05221].
  * (44) R.M. Fonseca, _Enumerating the operators of an effective field theory_ , _Phys. Rev. D_ 101 (2020) 035040 [1907.12584].
  * (45) W. Dekens, J. de Vries, K. Fuyuto, E. Mereghetti and G. Zhou, _Sterile neutrinos and neutrinoless double beta decay in effective field theory_ , _JHEP_ 06 (2020) 097 [2002.07182].
  * (46) W. Dekens, J. de Vries and T. Tong, _Sterile neutrinos with non-standard interactions in $\beta$\- and 0$\nu$$\beta$$\beta$-decay experiments_, _JHEP_ 08 (2021) 128 [2104.00140].
  * (47) J. de Vries, G. Li, M.J. Ramsey-Musolf and J.C. Vasquez, _Light Sterile Neutrinos, Left-Right Symmetry, and $0\nu\beta\beta$ Decay_, 2209.03031.
  * (48) V. Cirigliano et al., _Neutrinoless Double-Beta Decay: A Roadmap for Matching Theory to Experiment_ , 2203.12169.
  * (49) S. Bhattacharya and J. Wudka, _Dimension-seven operators in the standard model with right handed neutrinos_ , _Phys. Rev. D_ 94 (2016) 055022 [1505.05264].
  * (50) E.E. Jenkins, A.V. Manohar and P. Stoffer, _Low-Energy Effective Field Theory below the Electroweak Scale: Anomalous Dimensions_ , _JHEP_ 01 (2018) 084 [1711.05270].
  * (51) M. Chala and A. Titov, _One-loop running of dimension-six Higgs-neutrino operators and implications of a large neutrino dipole moment_ , _JHEP_ 09 (2020) 188 [2006.14596].
  * (52) A. Datta, J. Kumar, H. Liu and D. Marfatia, _Anomalous dimensions from gauge couplings in SMEFT with right-handed neutrinos_ , _JHEP_ 02 (2021) 015 [2010.12109].
  * (53) A. Datta, J. Kumar, H. Liu and D. Marfatia, _Anomalous dimensions from Yukawa couplings in SMNEFT: four-fermion operators_ , _JHEP_ 05 (2021) 037 [2103.04441].
  * (54) M. Hirsch, H.V. Klapdor-Kleingrothaus and S.G. Kovalenko, _New low-energy leptoquark interactions_ , _Phys. Lett. B_ 378 (1996) 17 [hep-ph/9602305].
  * (55) K. Bondarenko, A. Boyarsky, D. Gorbunov and O. Ruchayskiy, _Phenomenology of GeV-scale Heavy Neutral Leptons_ , _JHEP_ 11 (2018) 032 [1805.08567].
  * (56) P. Ballett, T. Boschi and S. Pascoli, _Heavy Neutral Leptons from low-scale seesaws at the DUNE Near Detector_ , _JHEP_ 03 (2020) 111 [1905.00284].
  * (57) P. Coloma, E. Fernández-Martínez, M. González-López, J. Hernández-García and Z. Pavlovic, _GeV-scale neutrinos: interactions with mesons and DUNE sensitivity_ , _Eur. Phys. J. C_ 81 (2021) 78 [2007.03701].
  * (58) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (59) LHCb collaboration, _Measurement of the $B_{c}^{-}$ meson production fraction and asymmetry in 7 and 13 TeV $pp$ collisions_, _Phys. Rev. D_ 100 (2019) 112006 [1910.13404].
  * (60) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (61) C. Bierlich et al., _A comprehensive guide to the physics and usage of pythia 8.3_ , 2203.11601.
  * (62) “Forward Physics in PYTHIA 8.” https://pythia.org/download/talks/SjostrandForward20.pdf.
  * (63) A. Atre, T. Han, S. Pascoli and B. Zhang, _The Search for Heavy Majorana Neutrinos_ , _JHEP_ 05 (2009) 030 [0901.3589].
  * (64) M. Cacciari, M. Greco and P. Nason, _The P(T) spectrum in heavy flavor hadroproduction_ , _JHEP_ 05 (1998) 007 [hep-ph/9803400].
  * (65) M. Cacciari, S. Frixione and P. Nason, _The p(T) spectrum in heavy flavor photoproduction_ , _JHEP_ 03 (2001) 006 [hep-ph/0102134].
  * (66) NA62 collaboration, _Search for heavy neutral lepton production in $K^{+}$ decays to positrons_, _Phys. Lett. B_ 807 (2020) 135599 [2005.09575].
  * (67) T2K collaboration, _Search for heavy neutrinos with the T2K near detector ND280_ , _Phys. Rev. D_ 100 (2019) 052006 [1902.07598].
  * (68) CHARM collaboration, _A Search for Decays of Heavy Neutrinos in the Mass Range 0.5-GeV to 2.8-GeV_ , _Phys. Lett. B_ 166 (1986) 473.
  * (69) G. Bernardi et al., _Further Limits On Heavy Neutrino Couplings_ , _Phys. Lett. B_ 203 (1988) 332.
  * (70) S.A. Baranov et al., _Search for heavy neutrinos at the IHEP-JINR neutrino detector_ , _Phys. Lett. B_ 302 (1993) 336.
  * (71) Belle collaboration, _Search for heavy neutrinos at Belle_ , _Phys. Rev. D_ 87 (2013) 071102 [1301.1105].
  * (72) DELPHI collaboration, _Search for neutral heavy leptons produced in Z decays_ , _Z. Phys. C_ 74 (1997) 57.
  * (73) V. Shtabovenko, R. Mertig and F. Orellana, _New Developments in FeynCalc 9.0_ , _Comput. Phys. Commun._ 207 (2016) 432 [1601.01167].
  * (74) V. Shtabovenko, R. Mertig and F. Orellana, _FeynCalc 9.3: New features and improvements_ , _Comput. Phys. Commun._ 256 (2020) 107478 [2001.04407].
  * (75) Y. Aoki et al., _FLAG Review 2021_ , 2111.09849.
  * (76) J.L. Rosner, S. Stone and R.S. Van de Water, _Leptonic decays of charged pseudoscalar mesons_ , _Review of Particle Physics 2021 update_ (2021) .
  * (77) C. Bourrely, I. Caprini and L. Lellouch, _Model-independent description of $B\rightarrow\pi l\nu$ decays and a determination of $|{V}_{ub}|$_, _Phys. Rev. D_ 79 (2009) 013008 [0807.2722].
  * (78) ETM collaboration, _Scalar and vector form factors of $D\to\pi(K)\ell\nu$ decays with $N_{f}=2+1+1$ twisted fermions_, _Phys. Rev. D_ 96 (2017) 054514 [1706.03017].
  * (79) M.A. Ivanov, J.G. Körner, J.N. Pandya, P. Santorelli, N.R. Soni and C.-T. Tran, _Exclusive semileptonic decays of D and D s mesons in the covariant confining quark model_, _Front. Phys. (Beijing)_ 14 (2019) 64401 [1904.07740].
  * (80) Y.-L. Wu, M. Zhong and Y.-B. Zuo, _$B_{(s)},D_{(s)}\to\pi$ , $K$, $\eta$, $\rho$, $K^{\ast}$, $\omega$, $\phi$ Transition Form Factors and Decay Rates with Extraction of the CKM parameters $|V_{ub}|$, $|V_{cs}|$, $|V_{cd}|$_, _Int. J. Mod. Phys. A_ 21 (2006) 6125 [hep-ph/0604007].
  * (81) A. Bharucha, D.M. Straub and R. Zwicky, _$B\to V\ell^{+}\ell^{-}$ in the Standard Model from light-cone sum rules_, _JHEP_ 08 (2016) 098 [1503.05534].