# Effects of heavy Majorana neutrinos on lepton flavor violating processes

G. Hernández-Tomé  CAFPE and Departamento de Física Teórica y del Cosmos,
Universidad de Granada, E–18071 Granada, Spain  J. I. Illana  CAFPE and
Departamento de Física Teórica y del Cosmos, Universidad de Granada, E–18071
Granada, Spain  G. López Castro  Departamento de Física, Centro de
Investigación y de Estudios Avanzados del Instituto Politécnico Nacional,
Apdo. Postal 14-740, 07000 México D.F., México  M. Masip  CAFPE and
Departamento de Física Teórica y del Cosmos, Universidad de Granada, E–18071
Granada, Spain  P. Roig  Departamento de Física, Centro de Investigación y de
Estudios Avanzados del Instituto Politécnico Nacional, Apdo. Postal 14-740,
07000 México D.F., México

###### Abstract

The observation of lepton flavor violating processes at colliders could be a
clear signal of a non-minimal neutrino sector. We define a 5-parameter model
with a pair of TeV fermion singlets and arbitrary mixings with the three
active neutrino flavors. Then we analyze several flavor violating transitions
($\ell\to\ell^{\prime}\gamma,\,\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$
or $\mu-e$ conversions in nuclei) and $Z\to\bar{\ell}\,\ell^{\prime}$ decays
induced by the presence of heavy neutrinos. In particular, we calculate all
the one-loop contributions to these processes and present their analytic
expressions. We focus on the genuine effects of the heavy Majorana masses,
comparing the results in that case with the ones obtained when the two heavy
neutrinos define a Dirac field. Finally, we use our results to update the
bounds on the heavy-light mixings in the neutrino sector.

##  I Introduction

In the original formulation of the standard model (SM) Glashow:1961tr ;
Weinberg:1967tq ; Salam:1968rm , the lepton flavor and the lepton number are
accidentally conserved quantities due to the assumption of massless neutrinos.
However, this framework must be extended to account for the well-established
evidence of neutrino oscillations Fukuda:1998mi ; Ahmad:2001an ; Ahmad:2002jz
, which implies non-zero masses and mixings for the active neutrinos. A
possible minimal extension is the so-called $\nu$SM Mohapatra:1998rq , which
adds right-handed components (gauge singlets) for the three neutrino families
and generates Dirac masses via Yukawa couplings with the Higgs doublet, just
like for all the other fermions. In the $\nu$SM, the mixing in the leptonic
sector is described by a $3\times 3$ unitary matrix called the PMNS matrix
Pontecorvo:1957qd ; Maki:1962mu , analogous to the CKM matrix of the quark
sector Cabibbo:1963yz ; Kobayashi:1973fv . Nevertheless, the $\nu$SM requires
extremely tiny Yukawa couplings to explain the observed masses, which suggests
that other mechanism may be at work. If, in addition to the Dirac mass terms
($m_{D}$) that combine them with the active neutrinos, the singlets have
Majorana masses ($m_{M}$) that define a new scale, then the tiny neutrino
masses appear naturally for a very large value of $m_{M}$ (i.e., $m_{M}\gg
m_{D}$). In this seesaw mechanism Minkowski:1977sc ; GellMann:1980vs ;
Mohapatra:1979ia  the new mass terms break lepton number. The physical states
after diagonalization of the mass matrix include light ($\nu$) and heavy ($N$)
sectors of Majorana neutrinos with masses

$\displaystyle m_{\nu}\approx m_{D}^{2}/m_{M},\quad m_{N}\approx m_{M}\gg
m_{\nu}.$ (1)

Like in the $\nu$SM, in this model the mixings among the three active families
may be large, as required from oscillation experiments, but the mixing with
the heavy fields is of order

$\displaystyle s_{\nu}\approx m_{D}/m_{M}\approx\sqrt{m_{\nu}/m_{N}}.$ (2)

Both in the $\nu$SM and this high scale (type I) seesaw model, the rate of
lepton flavor violating (LFV) processes at colliders is suppressed by a factor
of $(m_{\nu}/E)^{2}$, being $E$ the scale of the process. In the second
scenario LFV can also be mediated by the neutrinos in the heavy sector, but
the heavy-light mixing implies then a suppression of order of $(E/m_{N})^{2}$,
equally small. In particular, the LFV decays $\ell\to\ell^{\prime}\gamma$,
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$ and
$Z\to\bar{\ell}\ell^{\prime}$, where $\ell,$ $\ell^{\prime},\dots$ denote the
usual charged leptons ($\tau,$ $\mu$, $e$) will have a branching ratio below
$10^{-50}$ Petcov:1976ff ; Bilenky:1977du ; Cheng:1985bj ; Illana:1999ww ;
Illana:2000ic ; Hernandez-Tome:2018fbq ; Blackstone:2019njl . It is then
apparent that any experimental observation of flavor violation involving
charged leptons (cLFV) would unambiguously imply the existence of new physics
at the TeV scale in an extended neutrino sector Arganda:2004bz ; Dinh:2012bp ;
Dinh:2013vya ; Abada:2014cca ; Arganda:2014dta ; DeRomeri:2016gum ;
Lindner:2016bgg .

Well-motivated variants of the two minimal models described above include the
inverse seesaw Mohapatra:1986bd ; Bernabeu:1987gr  or the linear seesaw
Malinsky:2005bi . These scenarios allow for arbitrary masses in the heavy
neutrino sector and then unsuppressed heavy-light mixings, constrained only by
the experimental limits. They are usually known as low-scale seesaw models,
although the masses in both sectors are not necessarily correlated. They are
justified by approximate symmetries or some ansatz on the neutrino mass matrix
that relaxes the restriction in Eq. (2). This type of models may be adequate
in scenarios like little Higgs (the heavy Majorana in seesaw models would
introduce quadratic corrections to the Higgs mass delAguila:2005yi ;
delAguila:2017ugt ; delAguila:2019mvp ), supersymmetry Abada:2014kba ;
Arganda:2015naa ; Arganda:2015ija , TeV gravity models (with the cutoff right
above that scale) ArkaniHamed:1998rs ; Randall:1999ee  or composite Higgs
models Coito:2019wte . In the next section we present a simple model that
captures all the relevant effects that may appear in cLFV processes induced by
the presence of heavy neutrinos.

Another possibility in these scenarios that is interesting from the
phenomenological point of view is to test the Dirac or Majorana nature of the
neutrinos in the heavy sector through lepton number ($L$) violating processes
with $\Delta L=2$. Apart from the longly explored neutrinoless double-beta
decay Furry:1939qr ; Zeldovich:1981da , this has been undertaken in tau decays
$\tau^{-}\to\ell^{+}M_{1}^{-}M_{2}^{-}$ ($M_{1},M_{2}=\pi,\,K$ mesons)
Ilakovac:1995km ; Ilakovac:1995wc ; Gribanov:2001vv , meson decays
$M_{1}^{+}\to\ell_{1}^{+}\ell_{2}^{+}M_{2}^{-}$ and hyperon decays, like
$\Sigma^{-}\to\Sigma^{+}e^{-}e^{-}$, $\Sigma^{-}\to p\mu^{-}\mu^{-}$
Littenberg:1991rd ; Barbero:2002wm , etc. All these studies are based on
scenarios where the new sterile Majorana neutrinos have non-negligible mixings
and some of them require masses low enough to be produced on-shell (resonant-
enhancement).

Table 1: Present limits and future sensitivities for the branching ratios or capture rates of several LFV processes. We denote $Z\to\ell_{1}\ell_{2}\equiv Z\to\bar{\ell}_{1}\ell_{2}+\ell_{1}\bar{\ell}_{2}$, and similarly for $h$ decays. For a more extensive list including hadronic modes see Calibbi:2017uvl ; Amhis:2016xyh .  Reaction | Present Limit 90% C.L. | Future Sensitivity | Reaction | Present Limit 90% C.L. | Future Sensitivity
---|---|---|---|---|---
$\mu\to e\gamma$ |  $4.2\times 10^{-13}$ Adam:2013mnn  |  $6\times 10^{-14}$ Baldini:2018nnn  | $\mu\to ee\bar{e}$ |  $1.0\times 10^{-12}$ Bellgardt:1987du  |  $10^{-16}$ Blondel:2013ia
$\mu-e$ (Au) |  $7.0\times 10^{-13}$ Bertl:2006up  | — |  $\mu-e$ (Ti) |  $4.3\times 10^{-12}$ Bertl:2006up  |  $10^{-18}$ Alekou:2013eta
$\tau\to e\gamma$ |  $3.3\times 10^{-8}\,$ Aubert:2009ag  |  $3\times 10^{-9}$ Kou:2018nap  | $\tau\to\mu\gamma$ |  $4.4\times 10^{-8}\,$ Aubert:2009ag  |  $10^{-9}$ Kou:2018nap
$\tau\to ee\bar{e}$ |  $2.7\times 10^{-8}$ Tanabashi:2018oca  |  | $\tau\to\mu\mu\bar{\mu}$ |  $2.1\times 10^{-8}\,$ Tanabashi:2018oca  |
$\tau\to e\mu\bar{\mu}$ |  $2.7\times 10^{-8}\,$ Tanabashi:2018oca  |  $(2-5)\times 10^{-10}$ Kou:2018nap  | $\tau\to\mu e\bar{e}$ |  $1.8\times 10^{-8}\,$ Tanabashi:2018oca  |  $(2-5)\times 10^{-10}$ Kou:2018nap
$\tau\to ee\bar{\mu}$ |  $1.5\times 10^{-8}\,$ Tanabashi:2018oca  |  | $\tau\to\mu\mu\bar{e}$ |  $1.7\times 10^{-8}\,$ Tanabashi:2018oca  |
Reaction | Present Limit 95% C.L. | Future Sensitivity | Reaction | Present Limit 95% C.L. | Future Sensitivity
$Z\to\mu e$ |  $7.3\times 10^{-7}$ Nehrkorn:2017fyt  |  $10^{-10}$ Dam:2018rfz  | $h\to\mu e$ |  $3.4\times 10^{-4}$ Khachatryan:2016rke  | —
$Z\to\tau e$ |  $9.8\times 10^{-6}$ Akers:1995gz  |  $10^{-9}$ Dam:2018rfz  | $h\to\tau e$ |  $6.2\times 10^{-3}$ Sirunyan:2017xzt  | $5\times 10^{-4}\,$Cerri:2018ypt
$Z\to\tau\mu$ |  $1.2\times 10^{-5}$ Abreu:1996mj  |  $10^{-9}$ Dam:2018rfz  | $h\to\tau\mu$ |  $2.5\times 10^{-3}$ Sirunyan:2017xzt

Currently there is no evidence for cLFV, but intense experimental efforts have
provided strong limits on an extensive list of processes; some of them are
reported in Table 1. The sensitivity to these transitions will be considerably
improved in near-future experiments. The MEG-II and Mu3e experiments will
reach branching ratios of order $6\times 10^{-14}$ Baldini:2018nnn  and
$10^{-16}$ Blondel:2013ia  for $\mu\to e\gamma$ and $\mu\to ee\bar{e}$,
respectively, whereas the expected bounds from PRISM and COMET will be near
$10^{-18}$ Alekou:2013eta  for $\mu-e$ (Ti) conversion and $10^{-17}$
Kuno:2013mha  for $\mu-e$ (Al). For the third family, the bounds on the
$\tau\to\ell^{\prime}\gamma$ and
$\tau\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell^{\prime\prime\prime}}$
branching ratios could be improved by two orders of magnitude at Belle-II when
the experiment achieves its maximum luminosity Hays:2017ekz ; Kou:2018nap .
LHCb has already set a stringent limit (competitive with the present ones at
Belle) of $4.6\times 10^{-8}$ Aaij:2014azz  on the $\tau\to\mu\mu\bar{\mu}$
process. In its high-luminosity phase the LHC is expected to improve this
bound by one order of magnitude. Additionally, the possibility of running at
the $Z$ pole in the electron-positron version of a Future Circular Collider
(FCC-ee) Dam:2018rfz ; Blondel:2019yqr  or in the Circular Electron Positron
Collider (CEPC) CEPCStudyGroup:2018ghi  would improve the current limits on
$Z\to\bar{\ell}\ell^{\prime}$ by about four orders of magnitude. Finally, the
expected sensitivity of the HL-LHC ($3000\,{\rm fb}^{-1}$) will be around
$5\times 10^{-4}$ for both the $h\to e\tau$ and $h\to\mu\tau$ branching
fractions Cerri:2018ypt .

In this work we focus on the most phenomenologically relevant cLFV observables
in the framework of low-scale seesaw scenarios. In Section II we introduce a
model for the mixings of the active neutrinos with two singlet fermions
defining Majorana fields. The mass splitting between these heavy fields
parametrizes the breaking of lepton number; when the splitting vanishes the
heavy sector reduces to a single Dirac neutrino. In Section III we provide
detailed expressions for the amplitudes and decay rates of the processes under
consideration. In Section IV we use these observables to derive constraints on
the heavy-light mixing angles as a function of the masses of the two heavy
states. Our conclusions are given in Section V.

##  II A model for the heavy-light neutrino mixing

As mentioned before, in the usual type-I seesaw model with one Majorana
singlet per family the heavy-light mixings are correlated with the neutrino
masses: to obtain $m_{\nu}<1$ eV with $m_{D}\approx 1$ GeV we need
$m_{M}>10^{9}$ GeV, and this implies negligible heavy-light mixings,
$s_{\nu}<10^{-9}$. As it is well known by the practitioners, however, this is
no longer the case when the singlet fermions are introduced in pairs (see
Bolton:2019pcu  for a recent review). In particular, all the heavy-light
mixing effects can be captured by considering a model with just one extra
pair. Let us see how this works.

We take two bi-spinors $N$ and $N^{c}$ of left-handed chirality (undotted),
sterile and with opposite lepton number, and define the four-spinors

$\displaystyle N_{L}=\begin{pmatrix}N\\\ 0\end{pmatrix}\,;\quad
N_{R}=\begin{pmatrix}0\\\
\bar{N}^{c}\end{pmatrix}\,;\quad\nu_{Li}=\begin{pmatrix}\nu_{i}\\\
0\end{pmatrix}\,,$ (3)

where $\nu_{i=e,\mu,\tau}$ are the SM neutrinos. After the breaking of the
electroweak symmetry the SM charged leptons get masses through Yukawa
interactions with the Higgs field; the left-handed mass eigenstates are
obtained after a unitary transformation,

$\displaystyle\ell_{Li}\leadsto\sum_{j=1}^{3}U_{ij}^{\ell}\ell_{Lj},$ (4)

that we also perform in the space of the three active neutrinos $\nu_{Li}$.
Then we assume that in this basis the 5 Majorana fields
$\chi_{i}=\chi_{Li}+\chi_{Li}^{c}$ with
$\chi_{L}\equiv(\nu_{L1},\nu_{L2},\nu_{L3},N_{L},N^{c}_{R})$ have mass terms

$\displaystyle{\cal L}_{M}$
$\displaystyle=-\frac{1}{2}\overline{\chi^{c}_{Li}}\,{\cal M}\,\chi_{Lj}+{\rm
h.c.},$ (5)

with

$\displaystyle{\cal M}=\begin{pmatrix}0&0&0&0&m_{1}\\\ 0&0&0&0&m_{2}\\\
0&0&0&0&m_{3}\\\ 0&0&0&0&M\\\ m_{1}&m_{2}&m_{3}&M&\mu\end{pmatrix}.$ (6)

Notice that we have ordered the fields according to their lepton number
(positive for the first four neutrinos), and that the Majorana mass $\mu$
corresponds to the neutrino with negative lepton number (the fifth one). The
mass eigenstates are obtained diagonalizing this symmetric matrix by an
orthogonal transformation and applying a field redefinition ($\chi_{L4}\to
i\chi_{L4}$) to guarantee real and positive mass eigenvalues. Three
eigenvalues are zero ($m_{\chi_{1,2,3}}=0$) and the other two are

$\displaystyle m_{\chi_{4}}$ $\displaystyle={1\over
2}\left(\sqrt{4\left(m_{1}^{2}+m_{2}^{2}+m_{3}^{2}+M^{2}\right)+\mu^{2}}-\mu\right),$
(7) $\displaystyle m_{\chi_{5}}$ $\displaystyle={1\over
2}\left(\sqrt{4\left(m_{1}^{2}+m_{2}^{2}+m_{3}^{2}+M^{2}\right)+\mu^{2}}+\mu\right).$
(8)

Defining $m\equiv\sqrt{m_{1}^{2}+m_{2}^{2}+m_{3}^{2}}$ and
$M^{\prime}\equiv\sqrt{m^{2}+M^{2}}$, the mass eigenstates are given by the
replacement

$\displaystyle\chi_{Li}\leadsto U_{ij}^{\nu}\chi_{Lj},$ (9)

where the mixing matrix reads

$\displaystyle
U^{\nu}=\begin{pmatrix}-{m_{2}\over\sqrt{m_{1}^{2}+m_{2}^{2}}}&-{m_{1}m_{3}\over
m\sqrt{m_{1}^{2}+m_{2}^{2}}}&-{m_{1}M\over
mM^{\prime}}&-i{m_{1}m_{\chi_{5}}\over
M^{\prime}\sqrt{m_{\chi_{5}}^{2}+M^{\prime
2}}}&{m_{1}\over\sqrt{m_{\chi_{5}}^{2}+M^{\prime
2}}}\cr{m_{1}\over\sqrt{m_{1}^{2}+m_{2}^{2}}}&-{m_{2}m_{3}\over
m\sqrt{m_{1}^{2}+m_{2}^{2}}}&-{m_{2}M\over
mM^{\prime}}&-i{m_{2}m_{\chi_{5}}\over
M^{\prime}\sqrt{m_{\chi_{5}}^{2}+M^{\prime
2}}}&{m_{2}\over\sqrt{m_{\chi_{5}}^{2}+M^{\prime 2}}}\cr
0&{\sqrt{m_{1}^{2}+m_{2}^{2}}\over m}&-{m_{3}M\over
mM^{\prime}}&-i{m_{3}m_{\chi_{5}}\over
M^{\prime}\sqrt{m_{\chi_{5}}^{2}+M^{\prime
2}}}&{m_{3}\over\sqrt{m_{\chi_{5}}^{2}+M^{\prime 2}}}\cr 0&0&{m\over
M^{\prime}}&-i{Mm_{\chi_{5}}\over M^{\prime}\sqrt{m_{\chi_{5}}^{2}+M^{\prime
2}}}&{M^{\prime}\over\sqrt{m_{\chi_{5}}^{2}+M^{\prime 2}}}\cr
0&0&0&i{M^{\prime}\over\sqrt{m_{\chi_{5}}^{2}+M^{\prime
2}}}&{m_{\chi_{5}}\over\sqrt{m_{\chi_{5}}^{2}+M^{\prime 2}}}\end{pmatrix}.$
(10)

Several comments are in order:

  * •

$\mu$ is the only mass parameter breaking lepton number. This parameter
defines the mass splitting of the two heavy Majorana neutrinos: when $\mu=0$
both states form a heavy Dirac neutrino singlet of mass $M^{\prime}$.

  * •

The two heavy neutrinos ($N_{1,2}\equiv\chi_{4,5}$) have components of order
$m_{i}/M$ along the corresponding active neutrinos, where
$m_{i}=Y_{\nu_{i}}v/\sqrt{2}$ are Dirac mass terms coming from Yukawa
couplings with the SM Higgs doublet. If the heavy fields have TeV masses, the
heavy-light mixings can be as large as $\sim 0.1$ for couplings of order one.

  * •

The three (mostly) active neutrinos ($\chi_{1,2,3}$), to be identified with
the light neutrinos ($\nu_{1,2,3}$) observed so far, are exactly massless. A
deformation of the pattern in Eq. (6) with a Majorana mass $\mu^{\prime}$ for
$N_{L}$ (like in inverse seesaw models) would imply that one of these
neutrinos gets a mass $m_{\nu}\approx\mu^{\prime}(m/M)^{2}$. Since $m_{\nu}<1$
eV, however, the new term $\mu^{\prime}$ must be small and has no effect on
flavor physics (it does not change the heavy-light mixings). An analogous
argument applies to possible Dirac mass terms $m^{\prime}_{i}$ in the fourth
row/column.

  * •

The pattern that we propose must be understood as approximate: these 5 entries
are the dominant mass terms, and any deformation must respect that the third
neutrino gets a mass below 1 eV (i.e., it must be much smaller than these 5
terms). Notice also that no symmetry protects the entries assumed to be zero,
and that loop corrections will actually introduce contributions to all of them
(see Bolton:2019pcu ). The pattern must then be established at the loop level
where we work, which may require the addition of tree level terms canceling
radiative corrections. In particular, large values of $\mu$ would induce
1-loop values of $\mu^{\prime}$ that must be canceled to obtain the proposed
pattern. This fine tuning disappears for $\mu\to 0$, when the resulting
pattern is justified by the conservation of lepton number.

  * •

The generation of small masses and light-light mixings for the three active
neutrinos would require the addition of extra singlets. This could be
accommodated with the usual mechanisms (in $\nu$SM or Type I seesaw models) or
through the deformations described above (in inverse seesaw models). In any
case, it will not introduce sizeable heavy-light mixings.

We will trade the five arbitrary mass parameters in Eq. (10) for the masses of
the two heavy neutrinos and three heavy-light mixings,

$\displaystyle m_{N_{1}}\equiv m_{\chi_{4}},\quad m_{N_{2}}\equiv
m_{\chi_{5}},\quad s_{\nu_{i}}\equiv{m_{i}\over\sqrt{m_{N_{1}}m_{N_{2}}}}.$
(11)

The $5\times 5$ matrix $U_{ij}^{\nu}$ above will introduce tree-level charged
and neutral currents involving neutrinos:

$\displaystyle\mathcal{L}_{W^{\pm}}$
$\displaystyle=-\frac{g}{\sqrt{2}}W^{-}_{\mu}\sum_{i=1}^{3}\sum_{j=1}^{5}B_{ij}\,\bar{\ell}_{i}\gamma^{\mu}P_{L}\chi_{j}+{\rm
h.c.},$ (12) $\displaystyle\mathcal{L}_{Z}$
$\displaystyle=-\frac{g}{4c_{W}}Z_{\mu}\sum_{i,j=1}^{5}\bar{\chi}_{i}\gamma^{\mu}\left(C_{ij}P_{L}-C_{ij}^{*}P_{R}\right)\chi_{j},$
(13) $\displaystyle\mathcal{L}_{G^{\pm}}$
$\displaystyle=-\frac{g}{\sqrt{2}M_{W}}G^{-}\sum_{i=1}^{3}\sum_{j=1}^{5}B_{ij}$
$\displaystyle\qquad\qquad\quad\times\bar{\ell}_{i}\left(m_{\ell_{i}}P_{L}-m_{\chi_{j}}P_{R}\right)\chi_{j}+{\rm
h.c.},$ (14)

where $G^{\pm}$ is the charged would-be-Goldstone field, $g$ is the weak
coupling constant, $c_{W}=\cos\theta_{W}$ and
$P_{L,R}=\frac{1}{2}(1\mp\gamma_{5})$ are the left and right-handed
projectors, respectively. Notice that in Eq. (13), the neutral current induced
by the Majorana states involves couplings of different flavors with both left
and right-handed components.111For the case of heavy left-handed neutrinos
being sequential Dirac (active) neutrinos, replace $B_{ij}\to U^{\nu}_{ij}$,
$C_{ij}\to\delta_{ij}$, $C_{ij}^{*}\to 0$. The dimension of the rectangular
$B$ mixing matrix is $3\times 5$, whereas $C$ is a $5\times 5$ matrix,

$\displaystyle B_{ij}=\sum_{k=1}^{3}\delta_{ik}U^{\nu}_{kj}\,,\quad
C_{ij}=\sum_{k=1}^{3}{(U_{ki}^{\nu})}^{*}U_{kj}^{\nu}.$ (15)

One can see that the elements of these matrices involving heavy neutrinos can
be expressed in terms of heavy-light mixings and the squared mass ratio
$r=m^{2}_{N_{2}}/m^{2}_{N_{1}}$ as

$\displaystyle
B_{kN_{1}}=-\frac{i\,r^{\frac{1}{4}}}{\sqrt{1+r^{\frac{1}{2}}}}s_{\nu_{k}},\quad
B_{kN_{2}}=\frac{1}{\sqrt{1+r^{\frac{1}{2}}}}s_{\nu_{k}},$ (16) $\displaystyle
C_{N_{1}N_{1}}$
$\displaystyle=\frac{r^{\frac{1}{2}}}{1+r^{\frac{1}{2}}}\sum_{k=1}^{3}s_{\nu_{k}}^{2},\quad
C_{N_{2}N_{2}}=\frac{1}{1+r^{\frac{1}{2}}}\sum_{k=1}^{3}s_{\nu_{k}}^{2},$
$\displaystyle
C_{N_{1}N_{2}}=-C_{N_{2}N_{1}}=\frac{i\,r^{\frac{1}{4}}}{1+r^{\frac{1}{2}}}\sum_{k=1}^{3}s_{\nu_{k}}^{2}.$
(17)

These are the same as in Ilakovac:1994kj ; Illana:2000ic  up to an irrelevant
global phase for $B$. In addition, the matrices $B$ and $C$ satisfy some
identities that are essential to keep the renormalizability of the model:

$\displaystyle\sum_{k=1}^{5}B_{ik}B_{jk}^{*}$
$\displaystyle=\delta_{ij},\quad\sum_{k=1}^{3}B_{ki}^{*}B_{kj}=\sum_{k=1}^{5}C_{ik}C_{jk}^{*}=C_{ij},$
$\displaystyle\sum_{k=1}^{5}B_{ik}C_{kj}=B_{ij},$ (18)
$\displaystyle\sum_{k=1}^{5}m_{\chi_{k}}C_{ik}C_{jk}$
$\displaystyle=\sum_{k=1}^{5}m_{\chi_{k}}B_{ik}C_{kj}^{*}=\sum_{k=1}^{5}m_{\chi_{k}}B_{ik}B_{jk}=0.$
(19)

##  III LFV Processes

We now present the amplitudes and decays widths or transition rates for the
LFV processes $\ell\to\ell^{\prime}\gamma$, $Z\to\bar{\ell}\ell^{\prime}$,
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$ and
$\mu-e$ conversion in nuclei. All of them involve the effective interaction of
a neutral vector boson with a pair of on-shell fermions, $V\ell\ell^{\prime}$
($V=\gamma,\,Z$), through a loop with Majorana neutrinos. Since the $W$
couples only to left-handed fields, the effective $V\ell\ell^{\prime}$
vertices ($\ell\neq\ell^{\prime}$) can be written in terms of the following
form factors:

$\displaystyle i\Gamma_{\mu}^{\gamma}(q^{2})$
$\displaystyle=ie\big{[}F_{L}^{\gamma}(q^{2})\gamma_{\mu}P_{L}+i2F_{M}^{\gamma}(q^{2})P_{R}\sigma_{\mu\nu}q^{\nu}\big{]},$
(20) $\displaystyle i\Gamma_{\mu}^{Z}(q^{2})$
$\displaystyle=ie\big{[}F_{L}^{Z}(q^{2})\gamma_{\mu}P_{L}\big{]},$ (21)

where $q$ is the momentum of the $V$ boson. Actually, the most general Lorentz
structure for on-shell fermions contains two additional (anapole) form
factors, $F_{S}$ and $F_{P}$. However, they do not contribute when the $V$
boson is on-shell, due to the transversality condition
$q_{\mu}\epsilon^{\mu}=0$. The same happens for an off-shell $V$ boson when
the masses of the external fermions can be neglected Hollik:1998vz . On the
other hand, the dipole form factors (chirality flipping) are proportional to
the external lepton masses.

Figure 1: Generic penguin and box diagrams contributing to
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$
($f=\ell^{\prime\prime}=\ell^{\prime\prime\prime}$) and $\mu-e$ conversion in
nuclei ($f=u,d$).

In the limit $q^{2}\to 0$, appropriate for $\ell\to\ell^{\prime}\gamma$ and
for the penguin contributions to
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$ and
$\mu-e$ conversion (Fig. 1), we may write:

$\displaystyle F_{L}^{\gamma}(q^{2})\equiv q^{2}A_{1L},\quad
F_{M}^{\gamma}(q^{2})\simeq F_{M}^{\gamma}(0)\equiv\frac{m_{\ell}}{2}A_{2R}.$
(22)

The vector form factor $F_{L}^{\gamma}$ for an on-shell photon vanishes by
current conservation due to the electromagnetic gauge invariance, and hence
only the dipole form factor $F_{M}^{\gamma}$ contributes to
$\ell\to\ell^{\prime}\gamma$. Then the amplitude reads:

$\displaystyle\mathcal{M}(\ell\to\ell^{\prime}\gamma)=iem_{\ell}A_{2R}\bar{u}(p_{\ell^{\prime}})P_{R}\sigma^{\mu\nu}q_{\nu}u(p_{\ell})\epsilon_{\mu}^{\gamma}(q),$
(23)

where $\epsilon_{\mu}^{\gamma}$ is the photon polarization vector and we have
neglected the mass of the lighter lepton. The partial decay width is given by

$\displaystyle\Gamma(\ell\to\ell^{\prime}\gamma)=\alpha
m_{\ell}^{3}|F_{M}^{\gamma}(0)|^{2}.$ (24)

The $Z\to\bar{\ell}\ell^{\prime}$ decay proceeds through the
$Z\ell\ell^{\prime}$ vertex with $q^{2}=M_{Z}^{2}$. Here we can take both
external leptons as massless and ignore the corresponding dipole form factor
$F_{M}^{Z}$, hence omitted in (21). The amplitude is then given by

$\displaystyle\mathcal{M}(Z\to\bar{\ell}\ell^{\prime})=eF_{L}^{Z}(M_{Z}^{2})\bar{u}(p_{\ell^{\prime}})\gamma^{\mu}P_{L}v(p_{\ell})\epsilon_{\mu}^{Z}(q),$
(25)

where $\epsilon_{\mu}^{Z}$ is the $Z$ polarization vector and the partial
decay width is

$\displaystyle\Gamma(Z\to\bar{\ell}\ell^{\prime})=\frac{\alpha}{3}M_{Z}|F_{L}^{Z}(M_{Z}^{2})|^{2}.$
(26)

Table 2: Possible LFV $\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$ channels. Type | Flavors | $\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$
---|---|---
1 | $\ell\neq\ell^{\prime}=\ell^{\prime\prime}=\ell^{\prime\prime\prime}$ | $\mu\to ee\bar{e}$ | $\tau\to ee\bar{e}$ | $\tau\to\mu\mu\bar{\mu}$
2 | $\ell\neq\ell^{\prime}\neq\ell^{\prime\prime}=\ell^{\prime\prime\prime}$ |  | $\tau\to e\mu\bar{\mu}$ | $\tau\to\mu e\bar{e}$
3 | $\ell\neq\ell^{\prime}=\ell^{\prime\prime}\neq\ell^{\prime\prime\prime}$ |  | $\tau\to ee\bar{\mu}$ | $\tau\to\mu\mu\bar{e}$

Regarding
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$, we
distinguish the three types of decays in Table 2. Apart from the photon-
penguin and $Z$-penguin diagrams containing the effective $V\ell\ell^{\prime}$
vertices, these decays involve box diagrams (Fig. 1):

$\displaystyle\mathcal{M}(\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime})=\mathcal{M}_{\gamma}+\mathcal{M}_{Z}+\mathcal{M}_{B},$
(27)

where

$\displaystyle\mathcal{M}_{\gamma}$
$\displaystyle=\frac{e^{2}}{q^{2}}\bar{u}({p_{\ell^{\prime}}})\left(q^{2}A_{1L}\gamma_{\mu}P_{L}+im_{\ell}A_{2R}P_{R}\sigma_{\mu\nu}q^{\nu}\right)u({p_{\ell}})$
$\displaystyle\times\bar{u}({p_{\ell^{\prime\prime}}})\gamma^{\mu}v({p_{\ell^{\prime\prime\prime}}})-\left(\ell^{\prime}\leftrightarrow\ell^{\prime\prime}\right),$
(28) $\displaystyle\mathcal{M}_{Z}$
$\displaystyle=-\frac{e^{2}}{M_{Z}^{2}}F_{L}^{Z}(0)\bar{u}({p_{\ell^{\prime}}})\gamma_{\mu}P_{L}u({p_{\ell}})$
$\displaystyle\times\bar{u}({p_{\ell^{\prime\prime}}})\gamma^{\mu}(g^{Z}_{L}P_{L}+g^{Z}_{R}P_{R})v({p_{\ell^{\prime\prime\prime}}})-\left(\ell^{\prime}\leftrightarrow\ell^{\prime\prime}\right),$
(29) $\displaystyle\mathcal{M}_{B}$
$\displaystyle=e^{2}F_{B}\bar{u}({p_{\ell^{\prime}}})\gamma_{\mu}P_{L}u({p_{\ell}})\;\bar{u}({p_{\ell^{\prime\prime}}})\gamma^{\mu}P_{L}v({p_{\ell^{\prime\prime\prime}}}),$
(30)

with $g_{L,R}^{Z}$ the charged lepton couplings to the $Z$ boson in units of
$e$,

$\displaystyle g_{L}^{Z}=\frac{1}{2s_{W}c_{W}}(-1+2s^{2}_{W}),\quad
g_{R}^{Z}=\frac{s_{W}}{c_{W}},$ (31)

and the box diagrams are evaluated in the limit of zero external momenta.
Channels of type 3 receive only box contributions, as they require two flavor-
changing vertices. We have written the vector form factor $F_{L}^{\gamma}$ in
terms of $A_{1L}$ (22) to emphasize that the photon propagator cancels the
$q^{2}$ prefactor. However, the dipole form factor $F_{M}^{\gamma}$, written
here in terms of $A_{2R}$ (22), will introduce a logarithmic dependence with
the external lepton masses (they cannot be neglected) after phase-space
integration of the squared amplitude. Notice that crossed diagrams with
$\ell^{\prime}$ and $\ell^{\prime\prime}$ exchanged must be added, except for
channels of type 2. In the box amplitude, the form factor $F_{B}$ includes the
crossed contribution thanks to a Fierz identity (see Eq. 74 of Appendix A).
The expressions for the partial decay widths of
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$ as a
function of $A_{1L}$, $A_{2R}$, $F_{L}^{Z}(0)$ and $F_{B}$ are given in
Appendix B.

The $\mu-e$ conversion in nuclei follows from similar diagrams as
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$
replacing the last two leptons by a quark $q=u$ or $d$ (Fig. 1). It involves
the same photon-penguin and $Z$-penguin and a couple of new box form factors,
$F_{B_{u}}$, $F_{B_{d}}$:

$\displaystyle{\cal
M}_{B}^{q}=e^{2}F_{B_{q}}\bar{u}(p_{\mu})\gamma_{\mu}P_{L}u(p_{q})\;\bar{u}(p_{e})\gamma^{\mu}P_{L}u(p^{\prime}_{q}).$
(32)

The expressions for the $\mu-e$ conversion rate in nuclei as a function of
$A_{1L}$, $A_{2R}$, $F_{L}^{Z}(0)$, $F_{B_{u}}$ and $F_{B_{d}}$ are given in
Appendix C.

Figure 2: One-loop diagrams contributing to the $V\ell\ell^{\prime}$ vertex.
We have omitted here and elsewhere diagrams with would-be-Goldstone fields,
needed in the Feynman-’t Hooft gauge.

We have calculated in our model the one-loop contributions to the form factors
introduced above, in the Feynman-’t Hooft gauge and using dimensional
regularization. The effective LFV $V\ell\ell^{\prime}$ vertex is obtained from
the diagrams of Fig. 2 supplemented by similar ones with the $W^{\pm}$ fields
replaced by the would-be-Goldstone fields $G^{\pm}$. The resulting photon form
factors in the low $q^{2}$ limit are:

$\displaystyle F_{L}^{\gamma}(q^{2})$ $\displaystyle=\frac{\alpha_{W}}{8\pi
M_{W}^{2}}\sum_{i}^{5}B^{*}_{\ell
i}B_{\ell^{\prime}i}f_{L}^{\gamma}(x_{i};q^{2}),$ (33) $\displaystyle
F_{M}^{\gamma}(0)$ $\displaystyle=\frac{\alpha_{W}}{8\pi
M_{W}^{2}}\frac{m_{\ell}}{2}\sum_{i}^{5}B^{*}_{\ell
i}B_{\ell^{\prime}i}f_{M}^{\gamma}(x_{i}),$ (34)

where $\alpha_{W}=\alpha/s_{W}^{2}$, $x_{i}\equiv m_{\chi_{i}}^{2}/M_{W}^{2}$
and

$\displaystyle f_{L}^{\gamma}(x;q^{2})$
$\displaystyle=\bigg{[}\frac{x^{2}\left(x^{2}-10x+12\right)\ln x}{6(x-1)^{4}}$
$\displaystyle+\frac{\left(7x^{3}-x^{2}-12x\right)}{12(x-1)^{3}}-\frac{5}{9}\bigg{]}q^{2}+2M_{W}^{2}\Delta_{\epsilon},$
(35) $\displaystyle f_{M}^{\gamma}(x)$ $\displaystyle=\frac{3x^{3}\ln
x}{2(x-1)^{4}}-\frac{2x^{3}+5x^{2}-x}{4(x-1)^{3}}+\frac{5}{6}.$ (36)

The term $\Delta_{\epsilon}=1/\epsilon-\gamma_{E}+\ln
4\pi+\ln(\mu^{2}/M_{W}^{2})$ regulates the ultraviolet divergence in
$4-\epsilon$ dimensions and cancels in (33) for $\ell\neq\ell^{\prime}$ due to
the properties of $B$ (18). As expected, $F_{L}^{\gamma}(0)$ is zero.

From (34) one may derive the contribution $\delta a$ of heavy neutrinos to the
muon dipole moment anomaly, $(g-2)/2$. Subtracting that of light (massless)
neutrinos, it reads:

$\displaystyle\delta
a=\frac{\alpha_{W}}{4\pi}\frac{m^{2}_{\mu}}{M_{W}^{2}}\sum_{i=4}^{5}|B_{\mu
i}|^{2}\left[\frac{3x_{i}^{3}\ln
x_{i}}{2(x_{i}-1)^{4}}-\frac{2x_{i}^{3}+5x_{i}^{2}-x_{i}}{4(x_{i}-1)^{3}}\right]$
(37)

that is negative, enhancing the disagreement with the current experimental
measurement Tanabashi:2018oca , but anyway negligible because the prefactor is
$\sim 4\times 10^{-11}$, $|B_{\mu N_{i}}|^{2}\leq s_{\nu_{\mu}}^{2}\lesssim
10^{-3}$ and the absolute value of the remaining function is smaller than 0.5.

Regarding the $F_{L}^{Z}$ form factor of the effective $Z\ell\ell^{\prime}$
vertex we find

$\displaystyle F_{L}^{Z}(q^{2})$ $\displaystyle=\frac{\alpha_{W}}{8\pi
s_{W}c_{W}}\sum_{i,j}^{5}B^{*}_{\ell
i}B_{\ell^{\prime}j}\big{[}\delta_{ij}F(x_{i};q^{2})$
$\displaystyle+C^{*}_{ij}G(x_{i},x_{j};q^{2})+C_{ij}\sqrt{x_{i}x_{j}}H(x_{i},x_{j};q^{2})\big{]},$
(38)

where

$\displaystyle F(x;q^{2})$
$\displaystyle=2c_{W}^{2}\left[q^{2}\left(\overline{C}_{1}+\overline{C}_{2}+\overline{C}_{12}\right)-6\overline{C}_{00}+1\right]$
$\displaystyle-(1-2s_{W}^{2})x\overline{C}_{00}-2s_{W}^{2}xM_{W}^{2}\overline{C}_{0}$
$\displaystyle+\frac{1}{2}\left(1-2c_{W}^{2}\right)\left[(2+x)\bar{B}_{1}+1\right],$
(39) $\displaystyle G(x,y;q^{2})$
$\displaystyle=-q^{2}\left(C_{0}+C_{1}+C_{2}+C_{12}\right)+2C_{00}-1$
$\displaystyle-\frac{1}{2}xyM_{W}^{2}C_{0},$ (40) $\displaystyle H(x,y;q^{2})$
$\displaystyle=q^{2}C_{0}+\frac{1}{2}q^{2}C_{12}-C_{00}+\frac{1}{4},$ (41)

in full agreement with Illana:2000ic . Here we have used the following
shorthand notation for the standard Passarino-Veltman loop functions
Passarino:1978jh ,

$\displaystyle C(x,y)_{\\{0,1,2,12\\}}$
$\displaystyle=C_{\\{0,1,2,12\\}}(0,q^{2},0;M_{W}^{2},xM_{W}^{2},yM_{W}^{2}),$
(42) $\displaystyle\overline{C}_{\\{0,1,2,12\\}}(x)$
$\displaystyle=C_{\\{0,1,2,12\\}}(0,q^{2},0;xM_{W}^{2},M_{W}^{2},M_{W}^{2}),$
(43) $\displaystyle\bar{B}_{1}(x)$
$\displaystyle=B_{1}(0;xM_{W}^{2},M_{W}^{2}),$ (44)

defined with the same conventions as the computer packages LoopTools
Hahn:1998yk  and Collier Denner:2016kdg , that we have employed for numerical
evaluations. Analytic expressions for these functions in the low $q^{2}$
limit, appropriate for the $Z$-penguin contribution to
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$, can
be found in Ilakovac:1994kj  and have been cross-checked with the help of
Package-X Patel:2016fam . They are:

$\displaystyle F(x;0)$ $\displaystyle=\frac{5x^{2}\ln
x}{2(x-1)^{2}}\\!-\\!\frac{5x}{2(x-1)}+\frac{1}{4}\\!-\\!\left(\frac{5}{2}\\!-\\!2s_{W}^{2}\right)\Delta_{\epsilon},$
(45) $\displaystyle G(x,y;0)$
$\displaystyle=\frac{1}{2(x-y)}\left[\frac{(y-1)x^{2}\ln
x}{(x-1)}-\frac{(x-1)y^{2}\ln y}{(y-1)}\right]$
$\displaystyle+\frac{1}{2}\left(\Delta_{\epsilon}-\frac{1}{2}\right),$ (46)
$\displaystyle H(x,y;0)$ $\displaystyle=\frac{1}{4(x-y)}\left[\frac{x(x-4)\ln
x}{x-1}-\frac{y(y-4)\ln y}{y-1}\right]$
$\displaystyle-\frac{1}{4}\left(\Delta_{\epsilon}+\frac{1}{2}\right).$ (47)

The ultraviolet divergences cancel in (38) using the properties of the mixing
matrices (18) and (19).

Figure 3: Box diagrams contributing to
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$. The
diagram on the right introduces explicit LNV contributions.

The box form factors are all finite. The amplitude for
$\ell\to\ell^{\prime}\ell^{\prime\prime}\bar{\ell}^{\prime\prime\prime}$
receives the contribution of diagrams with explicit lepton number violating
(LNV) vertices (Fig. 3). To implement the LNV vertices we have followed the
algorithm in Denner:1992vza  that circumvents the explicit introduction of the
charge conjugation matrix in the Feynman rules and allows to use Dirac
propagators also for Majorana particles. In particular, the diagrams on the
right of Fig. 3 contain genuine LNV contributions from Majorana particles that
should vanish if lepton number is conserved. We have verified that this is
indeed the case when the two fermion singlets form a Dirac field, i.e. when
$\mu=0$ ($r=1$). The Lorentz structure of all box diagrams can be reduced to
the form in (30) after some algebra (see Appendix A). In agreement with
Ilakovac:1994kj , we find

$\displaystyle F_{B}$ $\displaystyle=\frac{\alpha_{W}}{16\pi
M_{W}^{2}s_{W}^{2}}$
$\displaystyle\times\Big{\\{}\sum_{i,j}^{5}\left[B^{*}_{\ell
i}B_{\ell^{\prime}i}B^{*}_{\ell^{\prime\prime\prime}j}B_{\ell^{\prime\prime}j}+(\ell^{\prime}\leftrightarrow\ell^{\prime\prime})\right]f_{B_{d}}(x_{i},x_{j})$
$\displaystyle\quad+\sum_{i,j}^{5}B^{*}_{\ell
i}B_{\ell^{\prime}j}B^{*}_{\ell^{\prime\prime\prime}i}B_{\ell^{\prime\prime}j}f_{B}^{\rm
LNV}(x_{i},x_{j})\Big{\\}},$ (48)

where

$\displaystyle f_{B_{d}}(x,y)$
$\displaystyle=\left(1+\frac{xy}{4}\right)\tilde{d}(x,y)-2xy\,d(x,y),$ (49)
$\displaystyle f^{\rm LNV}_{B}(x,y)$
$\displaystyle=\sqrt{xy}\left[2\tilde{d}(x,y)-(4+xy)d(x,y)\right],$ (50)

and

$\displaystyle\tilde{d}(x,y)$ $\displaystyle=\frac{x^{2}\ln
x}{(1-x)^{2}(y-x)}+\frac{y^{2}\ln y}{(1-y)^{2}(x-y)}$
$\displaystyle-\frac{1}{(1-x)(1-y)},$ (51) $\displaystyle d(x,y)$
$\displaystyle=\frac{x\ln x}{(1-x)^{2}(y-x)}+\frac{y\ln y}{(1-y)^{2}(x-y)}$
$\displaystyle-\frac{1}{(1-x)(1-y)}.$ (52)

Figure 4: Box diagrams contributing to $\mu-e$ conversion in nuclei.

The one-loop contributions to the box form factors of $\mu-e$ conversion come
from the diagrams in Fig. 4. We obtain:

$\displaystyle F_{B_{d}}$ $\displaystyle=\frac{\alpha_{W}}{16\pi
m_{W}^{2}s_{W}^{2}}\sum_{i}^{5}\sum_{j}^{3}B^{*}_{\mu
i}B_{ei}\left|V_{jd}\right|^{2}f_{B_{d}}(x_{i},x_{j}^{u}),$ (53)
$\displaystyle F_{B_{u}}$ $\displaystyle=\frac{\alpha_{W}}{16\pi
m_{W}^{2}s_{W}^{2}}\sum_{i}^{5}\sum_{j}^{3}B^{*}_{\mu
i}B_{ei}\left|V_{uj}\right|^{2}f_{B_{u}}(x_{i},x_{j}^{d}),$ (54)

where $x^{q}_{i}=m^{2}_{q_{i}}/M_{W}^{2}$, $V_{ij}$ is the CKM matrix and

$\displaystyle f_{B_{u}}(x,y)$
$\displaystyle=-\left(4+\frac{xy}{4}\right)\tilde{d}(x,y)+2xy\,d(x,y).$ (55)

Neglecting all quark masses, except that of the top quark, and defining
$x_{t}=m_{t}^{2}/m_{W}^{2}$, we may write:

$\displaystyle\sum_{j}^{3}$
$\displaystyle\left|V_{jd}\right|^{2}f_{B_{d}}(x_{i},x_{j}^{u})$
$\displaystyle=\left|V_{td}\right|^{2}\left[f_{B_{d}}(x_{i},x_{t})-f_{B_{d}}(x_{i},0)\right]-f_{B_{d}}(x_{i},0),$
(56) $\displaystyle\sum_{j}^{3}$
$\displaystyle\left|V_{uj}\right|^{2}f_{B_{u}}(x_{i},x_{j}^{d})=f_{B_{u}}(x_{i},0).$
(57)

In Appendix D we show how to express all these form factors in terms of the
contributions of heavy neutrinos only.

##  IV Numerical results

Next we analyze the predictions of our model for different values of its free
parameters: the three heavy-light mixings $s_{\nu_{e}}$, $s_{\nu_{\mu}}$ and
$s_{\nu_{\tau}}$, the mass $m_{N_{1}}$ of the lightest heavy neutrino and the
mass ratio $r=m_{N_{2}}^{2}/m_{N_{1}}^{2}$. If $r=1$ the two heavy Majorana
neutrinos become a single Dirac field.

In order to be consistent with perturbative unitarity, the Yukawa couplings
cannot exceed an upper limit. We will take

$\displaystyle
Y_{\nu_{i}}=\frac{\sqrt{2m_{N_{1}}m_{N_{2}}}}{v}s_{\nu_{i}}=\frac{\sqrt{2}}{v}m_{N_{1}}r^{1/4}s_{\nu_{i}}<\sqrt{4\pi}.$
(58)

This means that, given $m_{N_{1}}$ and $r$,

$\displaystyle s_{\nu_{i}}<\frac{\sqrt{2\pi}v}{m_{N_{1}}r^{1/4}},$ (59)

that constrains the mixings if $m_{N_{1}}r^{1/4}\gtrsim 620$ GeV. In
particular, $s_{\nu_{i}}<0.12$ for $m_{N_{1}}r^{1/4}=5$ TeV. Given the mixings
$s_{\nu_{i}}$ this condition also implies

$\displaystyle m_{N_{1}}r^{1/4}<\frac{\sqrt{2\pi}v}{\max\\{s_{\nu_{i}}\\}}.$
(60)

On the other hand, the heavy-light mixings must respect indirect constraints.
We take $2\sigma$ limits from the global fit to electroweak precision
observables and lepton flavor conserving processes in Fernandez-
Martinez:2016lgt , where the effects of extra neutrinos are encoded in
effective operators:222In a recent work Coutinho:2019aiy  a global fit to
modified neutrino couplings has been performed that alleviates the Cabibbo-
angle anomaly and is compatible with the bounds we use.

$\displaystyle s_{\nu_{e}}<0.050,\quad s_{\nu_{\mu}}<0.021,\quad
s_{\nu_{\tau}}<0.075.$ (61)

Then (60) implies that $m_{N_{1}}r^{1/4}<8.2$ TeV if all mixings are fixed to
the upper limits, but it could be larger otherwise.

###  IV.1 $\mu-e$ transitions

Figure 5: Contour plots in the $s_{\nu_{e}}-s_{\nu_{\mu}}$ plane assuming
degenerate heavy neutrino masses $m_{N_{1}}=m_{N_{2}}=m_{N}=0.1,\,5$ TeV that
saturate present limits (solid lines) and future sensitivities (dashed lines)
of several $\mu-e$ transitions. A conservative value $s_{\nu_{\tau}}=0$ is
assumed when needed. The black dot-dashed lines show current indirect limits.
The gray-shadowed region is beyond the pertubative limit
$Y_{\nu_{i}}^{2}<4\pi$.

LFV processes involving only $\mu-e$ transitions further constrain the masses
and mixings of the heavy neutrinos in our model. Let us first consider the
case of one singlet Dirac neutrino ($m_{N_{1}}=m_{N_{2}}=m_{N}$). Figure 5
shows the contours in the $s_{\nu_{e}}-s_{\nu_{\mu}}$ plane that saturate
present experimental bounds (solid lines) and the future sensitivities in
Table 1 (dashed lines) for $m_{N}=0.1,5$ TeV. In general the amplitudes for
$Z\to\bar{\mu}e$, $\mu\to ee\bar{e}$ and $\mu-e$ conversion in nuclei depend
on all three heavy-light mixing angles through the $Z$-penguin contribution,
that involves the $C_{ij}$ matrix elements (17). Here we have assumed
$s_{\nu_{\tau}}=0$, so the regions below the curves in Fig. 5 enclose the most
conservative values for $s_{\nu_{e}}$ and $s_{\nu_{\mu}}$ (i.e., a non-zero
$s_{\nu_{\tau}}$ would imply stronger bounds). Actually, $\mu\to e\gamma$ sets
stringent constrains only on the product $s_{\nu_{e}}s_{\nu_{\mu}}$; for
$m_{N}\gtrsim 1$ TeV we find that this process does not depend on the heavy
neutrino masses and that its branching fraction can be approximated by

$\displaystyle{\rm BR}(\mu\to
e\gamma)\approx\frac{3\alpha}{8\pi}s_{\nu_{e}}^{2}s_{\nu_{\mu}}^{2},$ (62)

which yields the conservative direct limit:

$\displaystyle s_{\nu_{e}}^{2}s_{\nu_{\mu}}^{2}<5.1\times 10^{-10}.$ (63)

Fig. 5 also reveals that in forthcoming experiments $\mu-e$ (Ti) and $\mu\to
ee\bar{e}$ will be more constraining than $\mu\to e\gamma$.

Another point that we would like to emphasize is that the amplitudes for
$\mu\to ee\bar{e}$, $\mu-e$ conversion in nuclei and $Z\to\mu e$ introduce
terms of order $s_{\nu_{i}}^{4}$ that cannot be ignored, since they imply a
strong quadratic dependence on the heavy neutrino masses. Indeed, these terms
dominate the amplitude when the splitting between the two heavy masses is
large. Our results differ then from those in Fernandez-Martinez:2016lgt ,
where as a first approximation the terms proportional to $s_{\nu_{i}}^{4}$ are
neglected.

From a phenomenological point of view it is also interesting to investigate
whether the model can accommodate values of the different observables
involving $\mu-e$ transitions near the current experimental bounds. What are
the maximum values of $\mu\to e\gamma$, $\mu\to ee\bar{e}$, $\mu-e$ (Ti) and
$Z\to\mu e$ consistent with all the bounds in heavy neutrino models? To answer
this question we have considered the cases with low ($r=1$) and high ($r=25$)
neutrino mass ratio and different values of $s_{\nu_{e}}s_{\nu_{\mu}}$,
$m_{N}$, and $s_{\nu_{\tau}}$ that respect the indirect and perturbative
limits.

Figure 6: Predictions normalized to current limits assuming fixed maximum
mixings $s_{\nu_{e}}^{\rm max}$ and $s_{\nu_{\mu}}^{\rm max}$ compatible with
$\mu\to e\gamma$ (63) and the indirect limits (61), for two values of the
neutrino mass-ratio $r$. The predictions for $\mu\to ee\bar{e}$ and $\mu-e$
(Ti) get constrained to the corresponding bands whose upper (lower) boundaries
are determined by $s_{\nu_{\tau}}=0\;(s_{\nu_{\tau}}^{\rm max})$. Higher
masses are forbidden by perturbative unitarity (60) for these mixings. The
ratio of BR$(Z\to\mu e)$ to current limits is always much smaller than the
others, below $10^{-6}$.

In Fig. 6 we plot the ratio of the different observables to their current
bound. In general, it is $\mu\to e\gamma$ the most constraining process, so in
the plot we set the maximum value of $s_{\nu_{e}}s_{\nu_{\mu}}$ compatible
with that process and vary the rest of parameters. The lower (upper) curves of
each band correspond to $s_{\nu_{\tau}}=0\;(s_{\nu_{\tau}^{\rm max}}=0.075)$,
whereas the drop in the $\mu-e$ conversion amplitude at neutrino masses
$\lesssim 2$ TeV is due to the opposite sign in the form factors $F_{B_{d}}$
and $F_{B_{u}}$.

We find that $\mu-e$ (Ti) may also saturate its present bounds if the neutrino
masses are large enough: $m_{N}>5.5$ TeV for a Dirac neutrino ($r=1$) and
$m_{N}>4.7$ TeV if $r=25$. In contrast, in these models the processes $\mu\to
ee\bar{e}$ and $Z\to\mu e$ can not reach their current experimental limits
consistently with $\mu\to e\gamma$ and $\mu-e$ (Ti) for any values of the free
parameters. We find

$\displaystyle{\rm BR}(\mu\to ee\bar{e})<5.2\times 10^{-13},$ (64)

and

$\displaystyle{\rm BR}(Z\to\mu e)<6.5\times 10^{-13}.$ (65)

These limits do not change for lower values of $s_{\nu_{e}}s_{\nu_{\mu}}$, so
our result implies that the observation at future experiments of any of these
processes at a rate between the current bounds and these upper limits would
exclude heavy neutrinos (both Dirac or Majorana) as a possible explanation.
Larger values of the Majorana mass ratio $r$ than the one shown in the second
plot of Fig. 6 would result in similar allowed regions just cutting off the
higher masses to meet the perturbative unitarity limit.

Looking at the improvement factor of the sensitivities in future experiments
(Table 1) we conclude that $\mu-e$ (Ti) will take the lead in constraining the
parameter space of our model, rather than $\mu\to e\gamma$, except for a tiny
region of masses between 1.5 and 2 TeV that would be probed better by $\mu\to
ee\bar{e}$.

###  IV.2 $\tau-e$ transitions

Figure 7: Maximum values of $\tau-e$ transition rates compatible with current
constraints on $\mu\to e\gamma$, $\mu-e$ (Ti) and the indirect limits (61).
These maximum rates get constrained to the corresponding bands whose lower
(upper) boundaries are determined by $r=1\;(r\gg 1)$.

The constraints on our model from current limits on LFV $\tau-e$ transitions
(Table 1) turn out to be less restrictive than those involving the first two
lepton families. In Fig. 7 we show our predictions for the maximum possible
rates for several LFV $\tau$ decays consistent with current bounds from
$\mu\to e\gamma$ and $\mu-e$ (Ti), the indirect limits (61) and perturbative
unitarity (60).

The maximum branching ratio for $\tau\to e\gamma$, independent of heavy
neutrino masses as for $\mu\to e\gamma$, is

$\displaystyle{\rm BR}(\tau\to e\gamma)<2.0\times 10^{-9},$ (66)

close but still below the future sensitivity of Belle-II.

The predictions for $\tau\to ee\bar{e}$ and $\tau\to e\mu\bar{\mu}$ are very
similar because the dominant contribution comes in both cases from the $Z$
penguin diagram. They can reach:

$\displaystyle{\rm BR}(\tau\to ee\bar{e})$ $\displaystyle<7.3\times 10^{-9},$
(67) $\displaystyle{\rm BR}(\tau\to e\mu\bar{\mu})$ $\displaystyle<6.0\times
10^{-9},$ (68)

which, unlike $\tau\to e\gamma$, are well within the expected sensitivity of
Belle-II. These maximum values correspond to a Dirac neutino singlet ($r=1$)
with a mass just below the perturbative limit for the largest posible mixings,
$m_{N}=8.2$ TeV. For smaller masses there is room for $r>1$ that enhance the
decay rates up to the upper part of the shaded band in Fig 7.

The decay $\tau\to ee\bar{\mu}$ (double flavor change) is generated through
box diagrams only, so its amplitude is always proportional to
$s_{\nu_{e}}^{2}s_{\nu_{\tau}}s_{\nu_{\mu}}$ and hence more suppressed than
the other channels (see lower plot of Fig 7), with a maximum at

$\displaystyle{\rm BR}(\tau\to ee\bar{\mu})<2.3\times 10^{-14}.$ (69)

Nevertheless, it is important to remark that this decay is sensitive to the
genuine effects of Majorana neutrinos encoded in the LNV vertices of one of
its box contributions. Unlike the other processes, the rate of this for two
non-degenerate Majorana neutrinos can be enhanced by more than two orders of
magnitude when compared to the case of a Dirac singlet. In fact, the maximum
branching ratio above is obtained for $r\approx 16.9$ and $m_{N_{1}}\approx
4.1$ TeV.

For $Z\to\tau e$ our model predicts

$\displaystyle{\rm BR}(Z\to\tau e)<6.0\times 10^{-8},$ (70)

which is at the reach of future circular colliders.

Concerning the maximum values for the rates of $\tau-\mu$ transitions, we get
similar results as above by exchanging $e$ and $\mu$ and applying some
approximate correction factors. In particular, taking the maximal possible
mixings from the indirect limits (61), the processes $\tau\to\mu\gamma$,
$\tau\to\mu\mu\bar{\mu}$, $\tau\to\mu e\bar{e}$ and $Z\to\tau\mu$ are
suppressed by a factor of $(s_{\nu_{\mu}}^{\rm max}/s_{\nu_{e}}^{\rm
max})^{2}\approx 0.18$, whereas $\tau\to\mu\mu\bar{e}$ is enhanced by
$s_{\nu_{e}}^{\rm max}/s_{\nu_{\mu}}^{\rm max}\approx 2.4$.

##  V Conclusions

The lepton sector of the SM is still poorly known. In particular, we do not
know whether the observed neutrinos are Dirac or Majorana particles or if the
sector includes additional fermion singlets (sterile neutrinos). These extra
neutrinos, if any, would enhance cLFV processes that are otherwise very
suppressed by the tiny masses of the observed neutrinos.

In this work we have introduced the simplest neutrino model that captures all
the effects that are relevant to these processes: a minimal number of Majorana
neutrino fields (three active and two sterile), two of them heavy and the rest
massless, allowing unsuppressed heavy-light mixings and the possibility of LNV
encoded in the mass splitting of the heavy states
($\mu=m_{N_{2}}-m_{N_{1}}=(\sqrt{r}-1)m_{N_{1}}$). Larger splittings, however,
imply an increasing amount of fine tuning required to cancel loop corrections
that deform the proposed pattern, which is only stable in the lepton number
conserving case with $\mu=0$. In any case, only five parameters describe the
model, that are expressed in terms of the two heavy masses and the three
heavy-light mixings ($s_{\nu_{k}}^{2}$, $k=e,\mu,\tau$). The model could be
perturbed to account for the extremely light neutrino masses ($m_{\nu_{i}}<1$
eV) and the observed PMNS mixings, but this would have no impact on cLFV.

We have explored the predictions of our model for the most relevant reactions
involving one or two flavor changes. We have presented analytical expressions
for all of them and calculated their expected rates compatible with present
direct and indirect limits. Our computation is exact at one loop, including
all orders in the heavy-light mixings, and the genuine Majorana effects have
been singled out. This work completes and updates previous results and is in
agreement with an effective field theory analysis of the most general seesaw
extension of the SM Coy:2018bxr . We conclude that forthcoming LFV experiments
will probe a significant fraction of the parameter space of models with heavy
Majorana neutrinos.

## Acknowledgments

We would like to thank F. del Águila and José Santiago for helpful
discussions. This work was supported in part by the Spanish Ministry of
Science, Innovation and Universities, under grants FPA2016-78220-C3-1,2,3-P
(fondos FEDER), and Junta de Andalucía, grants FQM 101 and SOMM17/6104/UGR.
G.H.T. wants to acknowledge financial support from Conacyt through the program
“Estancia Postdoctoral en el Extranjero”. The work of P.R. has been partially
funded by Conacyt through the project 250628 (Ciencia Básica) and Fondo SEP-
Cinvestav 2018 (project number 142).

##  Appendix A Useful identities

In the limit of zero external momenta the Lorentz structure of all box
diagrams can be reduced to the same form using several identities based on the
decomposition in the chiral basis of Dirac matrices, transpositions and Fierz
rearrangements:

$\displaystyle\bar{u}(p_{1})\gamma^{\mu}\gamma^{\alpha}$
$\displaystyle\gamma^{\nu}P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\nu}\gamma_{\alpha}\gamma_{\mu}P_{L}v(p_{3})$
$\displaystyle=4\,\bar{u}(p_{1})\gamma^{\mu}P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\mu}P_{L}v(p_{3}),$
(71) $\displaystyle\bar{u}(p_{1})\gamma^{\mu}\gamma^{\alpha}$
$\displaystyle\gamma^{\nu}P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\mu}\gamma_{\alpha}\gamma_{\nu}P_{L}v(p_{3})$
$\displaystyle=16\,\bar{u}(p_{1})\gamma^{\mu}P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\mu}P_{L}v(p_{3}),$
(72) $\displaystyle\bar{u}(p_{3})\gamma^{\mu}\gamma^{\nu}$ $\displaystyle
P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\mu}\gamma_{\nu}P_{L}v(p_{1})$
$\displaystyle=\;\;\;4\,\bar{u}(p_{3})P_{L}u(p)\;\bar{u}(p_{2})P_{R}v(p_{1})$
$\displaystyle=-4\,\bar{v}(p)P_{L}v(p_{3})\;\bar{u}(p_{2})P_{R}v(p_{1})$
$\displaystyle=-\tfrac{1}{2}\,\bar{v}(p)\gamma^{\mu}P_{R}v(p_{1})\;\bar{u}(p_{2})\gamma_{\mu}P_{L}v(p_{3})$
$\displaystyle=-\tfrac{1}{2}\,\bar{u}(p_{1})\gamma^{\mu}P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\mu}P_{L}v(p_{3}),$
(73) $\displaystyle\bar{u}(p_{2})\gamma^{\mu}P_{L}$ $\displaystyle
u(p)\;\bar{u}(p_{1})\gamma_{\mu}P_{L}v(p_{3})$
$\displaystyle=-\bar{u}(p_{1})\gamma^{\mu}P_{L}u(p)\;\bar{u}(p_{2})\gamma_{\mu}P_{L}v(p_{3}).$
(74)

##  Appendix B Partial decay widths for LFV three-body decays

Given the generic form factors in Eqs. (27), (28), (29) and (30) the
expressions for the partial widths of the three types of decays in Table 2 are
delAguila:2019htj :

$\displaystyle\Gamma_{1}$
$\displaystyle=\frac{\alpha^{2}m_{\ell}^{5}}{96\pi}\Big{\\{}3\left|A_{1L}\right|^{2}+2\left|A_{2R}\right|^{2}\left(8\ln\frac{m_{\ell}}{m_{\ell^{\prime\prime}}}-13\right)$
$\displaystyle+2\left|F_{LL}\right|^{2}+\left|F_{LR}\right|^{2}+\tfrac{1}{2}\left|F_{B}\right|^{2}-\big{[}6A_{1L}A_{2R}^{*}-F_{LL}F_{B}^{*}$
$\displaystyle-(A_{1L}-2A_{2R})(2F_{LL}^{*}+F_{LR}^{*}+F_{B}^{*})+{\rm
h.c.}\big{]}\Big{\\}},$ (75) $\displaystyle\Gamma_{2}$
$\displaystyle=\frac{\alpha^{2}m_{\ell}^{5}}{96\pi}\Big{\\{}2\left|A_{1L}\right|^{2}+4\left|A_{2R}\right|^{2}\left(4\ln\frac{m_{\ell}}{m_{\ell^{\prime\prime}}}-7\right)$
$\displaystyle+\left|F_{LL}\right|^{2}+\left|F_{LR}\right|^{2}+\left|F_{B}\right|^{2}-\big{[}4A_{1L}A_{2R}^{*}-\tfrac{1}{2}F_{LL}F_{B}^{*}$
$\displaystyle-(A_{1L}-2A_{2R})(F_{LL}^{*}+F_{LR}^{*}+\tfrac{1}{2}F_{B}^{*})+{\rm
h.c.}\big{]}\Big{\\}},$ (76) $\displaystyle\Gamma_{3}$
$\displaystyle=\frac{\alpha^{2}m_{\ell}^{5}}{192\pi}\left|F_{B}\right|{}^{2},$
(77)

where

$\displaystyle F_{LL}=-\frac{F_{L}^{Z}(0)g_{L}^{Z}}{M_{Z}^{2}},\quad
F_{LR}=-\frac{F_{L}^{Z}(0)g_{R}^{Z}}{M_{Z}^{2}}.$ (78)

and the $Z$ couplings $g_{L,R}^{Z}$ are given in Eq. (31).

##  Appendix C The $\mu-e$ conversion rate

In terms of the form factors in Eqs. (27), (28), (29) and (32), the $\mu-e$
conversion rate in a nucleus with $Z$ protons and $N=A-Z$ neutrons is given by

$\displaystyle{\cal R}$ $\displaystyle=\frac{\alpha^{5}Z_{\rm
eff}^{4}}{\Gamma_{\rm
Capt}Z}F_{P}^{2}m_{\mu}^{5}\big{|}2Z\left(A_{1L}+A_{2R}\right)$
$\displaystyle-\left(2Z+N\right)(F_{LL}^{u}+F_{LR}^{u}+B_{L}^{u})$
$\displaystyle-\left(Z+2N\right)(F_{LL}^{d}+F_{LR}^{d}+B_{L}^{d})\big{|}^{2},$
(79)

where

$\displaystyle F_{LL}^{q}$
$\displaystyle=-\frac{F_{L}^{Z}(0)g_{Lq}^{Z}}{M_{Z}^{2}},\quad
F_{LR}^{q}=-\frac{F_{L}^{Z}(0)g_{Rq}^{Z}}{M_{Z}^{2}},$ (80) $\displaystyle
g_{Lu}^{Z}$ $\displaystyle=\frac{1-\frac{4}{3}s_{W}^{2}}{2s_{W}c_{W}},\quad
g_{Ru}^{Z}=-\frac{2s_{W}}{3c_{W}},$ (81) $\displaystyle g_{Ld}^{Z}$
$\displaystyle=\frac{-1+\frac{2}{3}s_{W}^{2}}{2s_{W}c_{W}},\quad
g_{Rd}^{Z}=\frac{s_{W}}{3c_{W}},$ (82)

and the rest of parameters are in Table 3.

Table 3: Input parameters for different nuclei Kitano:2002mt ; Suzuki:1987jf . Nucleus | $N$ | $Z$ | $Z_{\rm eff}$ | $F_{P}$ |  $\Gamma_{\rm capt}$ [GeV]
---|---|---|---|---|---
${27\atop 13}$Al | 14 | 13 | 11.5 | 0.64 | $4.6\times 10^{-19}$
${48\atop 22}$Ti | 26 | 22 | 17.6 | 0.54 | $1.7\times 10^{-18}$
${197\atop 79}$Au | 118 | 79 | 33.5 | 0.16 | $8.6\times 10^{-18}$

##  Appendix D Expressions in terms of massive neutrinos only

Using the relations between mixings matrices (18) and (19) with
$\ell\neq\ell^{\prime}$, one may write the generic contributions to all form
factors in terms of massive neutrinos only, $m_{N_{1}}=m_{\chi_{4}}$,
$m_{N_{2}}=m_{\chi_{5}}$:

$\displaystyle\sum_{i}^{5}B_{\ell i}^{*}B_{\ell^{\prime}i}f(x_{i})$
$\displaystyle=\sum_{i}^{2}B_{\ell
N_{i}}^{*}B_{\ell^{\prime}N_{i}}\left[f(x_{N_{i}})-f(0)\right],$ (83)
$\displaystyle\sum_{i,j}^{5}$ $\displaystyle B_{\ell
i}^{*}B_{\ell^{\prime}j}C_{ij}^{*}f(x_{i},x_{j})=\sum_{i,j}^{2}B_{\ell
N_{i}}^{*}B_{\ell^{\prime}N_{j}}$
$\displaystyle\times\big{\\{}\delta_{N_{i}N_{j}}\left[f(x_{N_{i}},0)+f(0,x_{N_{j}})-2f(0,0)\right]$
$\displaystyle+C_{N_{i}N_{j}}^{*}\big{[}f(x_{N_{i}},x_{N_{j}})-f(x_{N_{i}},0)-f(0,x_{N_{j}})$
$\displaystyle\qquad\qquad+f(0,0)\big{]}\big{\\}},$ (84)
$\displaystyle\sum_{i,j}^{5}$ $\displaystyle B_{\ell
i}^{*}B_{\ell^{\prime}i}B_{\ell^{\prime\prime\prime}j}^{*}B_{\ell^{\prime\prime}j}f(x_{i},x_{j})$
$\displaystyle=\sum_{i,j}^{2}B_{\ell
N_{i}}^{*}B_{\ell^{\prime}N_{j}}\delta_{\ell^{\prime\prime}\ell^{\prime\prime\prime}}\delta_{N_{i}N_{j}}\left[f(x_{N_{i}},0)-f(0,0)\right]$
$\displaystyle+\sum_{i,j}^{2}B_{\ell
N_{i}}^{*}B_{\ell^{\prime}N_{i}}B_{\ell^{\prime\prime\prime}N_{j}}^{*}B_{\ell^{\prime\prime}N_{j}}\big{[}f(x_{N_{i}},x_{N_{j}})$
$\displaystyle\qquad\qquad-f(x_{N_{i}},0)-f(0,x_{N_{i}})+f(0,0)\big{]}.$ (85)

## References

  * (1) S. L. Glashow, Nucl. Phys. 22, 579 (1961).
  * (2) S. Weinberg, Phys. Rev. Lett. 19, 1264 (1967).
  * (3) A. Salam, Conf. Proc. C 680519, 367 (1968).
  * (4) Y. Fukuda et al. [Super-Kamiokande Collaboration], Phys. Rev. Lett. 81, 1562 (1998) [hep-ex/9807003].
  * (5) Q. R. Ahmad et al. [SNO Collaboration], Phys. Rev. Lett. 87, 071301 (2001) [nucl-ex/0106015].
  * (6) Q. R. Ahmad et al. [SNO Collaboration], Phys. Rev. Lett. 89, 011301 (2002) [nucl-ex/0204008].
  * (7) R. N. Mohapatra and P. B. Pal, World Sci. Lect. Notes Phys. 60, 1 (1998) [World Sci. Lect. Notes Phys. 72, 1 (2004)].
  * (8) B. Pontecorvo, Sov. Phys. JETP 7, 172 (1958) [Zh. Eksp. Teor. Fiz. 34, 247 (1957)].
  * (9) Z. Maki, M. Nakagawa and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
  * (10) N. Cabibbo, Phys. Rev. Lett. 10, 531 (1963).
  * (11) M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 49, 652 (1973).
  * (12) P. Minkowski, Phys. Lett. 67B, 421 (1977).
  * (13) M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927 315 (1979) [arXiv:1306.4669 [hep-th]].
  * (14) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * (15) S. T. Petcov, Sov. J. Nucl. Phys. 25, 340 (1977) [Yad. Fiz. 25, 641 (1977)] Erratum: [Sov. J. Nucl. Phys. 25, 698 (1977)] Erratum: [Yad. Fiz. 25, 1336 (1977)].
  * (16) S. M. Bilenky, S. T. Petcov and B. Pontecorvo, Phys. Lett. 67B, 309 (1977).
  * (17) T. P. Cheng and L. F. Li, Oxford, Uk: Clarendon (1984) 536 P. (Oxford Science Publications).
  * (18) J. I. Illana, M. Jack and T. Riemann, hep-ph/0001273.
  * (19) J. I. Illana and T. Riemann, Phys. Rev. D 63, 053004 (2001) [hep-ph/0010193].
  * (20) G. Hernández-Tomé, G. López Castro and P. Roig, Eur. Phys. J. C 79, no. 1, 84 (2019) [arXiv:1807.06050 [hep-ph]].
  * (21) P. Blackstone, M. Fael and E. Passemar, arXiv:1912.09862 [hep-ph].
  * (22) E. Arganda, A. M. Curiel, M. J. Herrero and D. Temes, Phys. Rev. D 71, 035011 (2005) [hep-ph/0407302].
  * (23) D. N. Dinh, A. Ibarra, E. Molinaro and S. T. Petcov, JHEP 1208 125 (2012) Erratum: [JHEP 1309 023 (2013)] [arXiv:1205.4671 [hep-ph]].
  * (24) D. N. Dinh and S. T. Petcov, JHEP 1309, 086 (2013) [arXiv:1308.4311 [hep-ph]].
  * (25) A. Abada, V. De Romeri, S. Monteil, J. Orloff and A. M. Teixeira, JHEP 1504, 051 (2015) [arXiv:1412.6322 [hep-ph]].
  * (26) E. Arganda, M. J. Herrero, X. Marcano and C. Weiland, Phys. Rev. D 91, no. 1, 015001 (2015) [arXiv:1405.4300 [hep-ph]].
  * (27) V. De Romeri, M. J. Herrero, X. Marcano and F. Scarcella, Phys. Rev. D 95, no. 7, 075028 (2017) [arXiv:1607.05257 [hep-ph]].
  * (28) M. Lindner, M. Platscher and F. S. Queiroz, Phys. Rept. 731, 1 (2018) [arXiv:1610.06587 [hep-ph]].
  * (29) R. N. Mohapatra and J. W. F. Valle, Phys. Rev. D 34, 1642 (1986).
  * (30) J. Bernabéu, A. Santamaría, J. Vidal, A. Méndez and J. W. F. Valle, Phys. Lett. B 187, 303 (1987).
  * (31) M. Malinsky, J. C. Romao and J. W. F. Valle, Phys. Rev. Lett. 95, 161801 (2005) [hep-ph/0506296].
  * (32) F. del Aguila, M. Masip and J. L. Padilla, Phys. Lett. B 627, 131 (2005) [hep-ph/0506063].
  * (33) F. del Aguila, L. Ametller, J. I. Illana, J. Santiago, P. Talavera and R. Vega-Morales, JHEP 1708, 028 (2017) Erratum: [JHEP 1902, 047 (2019)] [arXiv:1705.08827 [hep-ph]].
  * (34) F. del Aguila, J. I. Illana, J. M. Perez-Poyatos and J. Santiago, JHEP 1912, 154 (2019) [arXiv:1910.09569 [hep-ph]].
  * (35) A. Abada, M. E. Krauss, W. Porod, F. Staub, A. Vicente and C. Weiland, JHEP 1411, 048 (2014) [arXiv:1408.0138 [hep-ph]].
  * (36) E. Arganda, M. J. Herrero, X. Marcano and C. Weiland, Phys. Rev. D 93, no. 5, 055010 (2016) [arXiv:1508.04623 [hep-ph]].
  * (37) E. Arganda, M. J. Herrero, X. Marcano and C. Weiland, Phys. Lett. B 752, 46 (2016) [arXiv:1508.05074 [hep-ph]].
  * (38) N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B 429, 263 (1998) [hep-ph/9803315].
  * (39) L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 3370 (1999) [hep-ph/9905221].
  * (40) W. H. Furry, Phys. Rev. 56, 1184 (1939).
  * (41) Y. B. Zeldovich and M. Y. Khlopov, Pisma Zh. Eksp. Teor. Fiz. 34, 148 (1981).
  * (42) L. Coito, C. Faubel and A. Santamaria, arXiv:1912.10001 [hep-ph].
  * (43) A. Ilakovac, B. A. Kniehl and A. Pilaftsis, Phys. Rev. D 52, 3993 (1995) [hep-ph/9503456].
  * (44) A. Ilakovac, Phys. Rev. D 54, 5653 (1996) [hep-ph/9608218].
  * (45) V. Gribanov, S. Kovalenko and I. Schmidt, Nucl. Phys. B 607, 355 (2001) [hep-ph/0102155].
  * (46) L. S. Littenberg and R. E. Shrock, Phys. Rev. D 46, R892 (1992).
  * (47) C. Barbero, G. López Castro and A. Mariano, Phys. Lett. B 566, 98 (2003) [nucl-th/0212083].
  * (48) Y. Amhis et al. [HFLAV Collaboration], Eur. Phys. J. C 77, no. 12, 895 (2017) [arXiv:1612.07233 [hep-ex]].
  * (49) L. Calibbi and G. Signorelli, Riv. Nuovo Cim. 41, 1 (2018).
  * (50) J. Adam et al. [MEG Collaboration], Phys. Rev. Lett. 110, 201801 (2013) [arXiv:1303.0754 [hep-ex]].
  * (51) U. Bellgardt et al. [SINDRUM Collaboration], Nucl. Phys. B 299, 1 (1988).
  * (52) W. H. Bertl et al. [SINDRUM II Collaboration], Eur. Phys. J. C 47, 337 (2006).
  * (53) B. Aubert et al. [BaBar Collaboration], Phys. Rev. Lett. 104, 021802 (2010) [arXiv:0908.2381 [hep-ex]].
  * (54) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, no. 3, 030001 (2018).
  * (55) A. Nehrkorn [CMS Collaboration], Nucl. Part. Phys. Proc. 287-288, 160 (2017).
  * (56) R. Akers et al. [OPAL Collaboration], Z. Phys. C 67, 555 (1995).
  * (57) P. Abreu et al. [DELPHI Collaboration], Z. Phys. C 73, 243 (1997).
  * (58) V. Khachatryan et al. [CMS Collaboration], Phys. Lett. B 763, 472 (2016) [arXiv:1607.03561 [hep-ex]].
  * (59) A. M. Sirunyan et al. [CMS Collaboration], JHEP 1806, 001 (2018) [arXiv:1712.07173 [hep-ex]].
  * (60) A. M. Baldini et al. [MEG II Collaboration], Eur. Phys. J. C 78, no. 5, 380 (2018) [arXiv:1801.04688 [physics.ins-det]].
  * (61) A. Blondel et al., arXiv:1301.6113 [physics.ins-det].
  * (62) A. Alekou et al., arXiv:1310.0804 [physics.acc-ph].
  * (63) E. Kou et al. [Belle-II Collaboration], arXiv:1808.10567 [hep-ex]. To be published in PTEP.
  * (64) M. Dam, SciPost Phys. Proc. 1, 041 (2019) [arXiv:1811.09408 [hep-ex]].
  * (65) A. Cerri et al., arXiv:1812.07638 [hep-ph].
  * (66) Y. Kuno [COMET Collaboration], PTEP 2013, 022C01 (2013).
  * (67) C. Hays, M. Mitra, M. Spannowsky and P. Waite, JHEP 1705, 014 (2017) [arXiv:1701.00870 [hep-ph]].
  * (68) R. Aaij et al. [LHCb Collaboration], JHEP 1502, 121 (2015) [arXiv:1409.8548 [hep-ex]].
  * (69) A. Blondel et al., arXiv:1906.02693 [hep-ph].
  * (70) J. B. Guimaraes da Costa et al. [CEPC Study Group], arXiv:1811.10545 [hep-ex].
  * (71) P. D. Bolton, F. F. Deppisch and P. S. B. Dev, arXiv:1912.03058 [hep-ph].
  * (72) A. Ilakovac and A. Pilaftsis, Nucl. Phys. B 437, 491 (1995) [hep-ph/9403398].
  * (73) W. Hollik, J. I. Illana, S. Rigolin, C. Schappacher and D. Stockinger, Nucl. Phys. B 551, 3 (1999) Erratum: [Nucl. Phys. B 557, 407 (1999)] [hep-ph/9812298].
  * (74) G. Passarino and M. J. G. Veltman, Nucl. Phys. B 160, 151 (1979).
  * (75) T. Hahn and M. Pérez-Victoria, Comput. Phys. Commun. 118, 153 (1999) [hep-ph/9807565].
  * (76) A. Denner, S. Dittmaier and L. Hofer, Comput. Phys. Commun. 212, 220 (2017) [arXiv:1604.06792 [hep-ph]].
  * (77) H. H. Patel, Comput. Phys. Commun. 218, 66 (2017) [arXiv:1612.00009 [hep-ph]].
  * (78) A. Denner, H. Eck, O. Hahn and J. Kublbeck, Nucl. Phys. B 387, 467 (1992).
  * (79) E. Fernández-Martínez, J. Hernández-García and J. López-Pavon, JHEP 1608, 033 (2016) [arXiv:1605.08774 [hep-ph]].
  * (80) A. M. Coutinho, A. Crivellin and C. A. Manzari, arXiv:1912.08823 [hep-ph].
  * (81) R. Coy and M. Frigerio, Phys. Rev. D 99 (2019) no.9, 095040 [arXiv:1812.03165 [hep-ph]].
  * (82) F. del Aguila, L. Ametller, J. I. Illana, J. Santiago, P. Talavera and R. Vega-Morales, JHEP 1907, 154 (2019) [arXiv:1901.07058 [hep-ph]].
  * (83) R. Kitano, M. Koike and Y. Okada, Phys. Rev. D 66, 096002 (2002) Erratum: [Phys. Rev. D 76, 059902 (2007)] [hep-ph/0203110].
  * (84) T. Suzuki, D. F. Measday and J. P. Roalsvig, Phys. Rev. C 35, 2212 (1987).