#  Connecting $(g-2)_{\mu}$ to neutrino mass in the extended neutrinophilic
2HDM

A. L. Cherchiglia<EMAIL_ADDRESS>Instituto de Física Gleb Wataghin,
Universidade Estadual de Campinas, Campinas-SP, Brazil  Departamento de Física
Teórica y del Cosmos, Universidad de Granada, Campus de Fuentenueva, E–18071
Granada, Spain  G. De Conto<EMAIL_ADDRESS>Centro de Ciências
Naturais e Humanas, Universidade Federal do ABC, Santo André-SP, Brasil  C. C.
Nishi<EMAIL_ADDRESS>Centro de Matemática, Computação e Cognição,
Universidade Federal do ABC, Santo André-SP, Brasil

###### Abstract

One simple way to lower the scale of the seesaw mechanism that generates
neutrino masses is to attribute part of their smallness to a suppressed vacuum
expectation value of a second Higgs doublet as in the neutrinophilic 2HDM or
in the type IB seesaw model. On that structure we add one charged singlet
scalar to induce a chirally enhanced contribution to $(g-2)_{\mu}$ with the
same righthanded neutrinos of the seesaw. We discuss the interplay of
generating the necessary contribution to the latter with lepton flavor
violation which is also necessarily brought to low scale. We show that it is
possible to explain $(g-2)_{\mu}$ even for heavy neutrino masses of order of a
few TeV.

##  I Introduction

The observation of neutrino oscillations firmly established that neutrinos
have tiny masses and mix in the weak charged current ohlsson . These
properties clearly demonstrate that family lepton numbers associated to the
lepton flavors $e,\mu,\tau$ are not conserved in nature. The origin of both
family lepton number violation and of neutrino masses can be naturally
attributed to the dimension five Weinberg operator weinberg.op  pointing to a
natural scale of $10^{12}\,\mathrm{GeV}$. If the new physics violating family
lepton number appears only at this scale, then its effects, visible primarily
on dimension six operators, are expected to be unobservably small. These
effects include those of charged lepton flavor violation (CLFV) that are
expected to be probed with much greater precision in the coming years.

In contrast, the discrepancy between the experimental and SM prediction for
the muon anomalous magnetic moment, $(g-2)_{\mu}$ or $a_{\mu}$, is a
persistent anomaly that requires TeV scale (or lower) new physics coupling to
the muon that may or may not violate lepton flavor. A dedicated program to
decrease the experimental uncertainty by a factor of four is currently ongoing
in Fermilab, with the first data analysis released in 2021, and the second in
2023. The value for $a_{\mu}$ combined from the results obtained in Fermilab
and Brookhaven is Muong-2:2021ojo  $a_{\mu}^{\text{Exp}}=(11659205.9\pm
2.2)\times 10^{-10},$ while the SM prediction from the Muon $g-2$ Theory
Initiative White Paper Aoyama:2020ynm ; g-2:theory  is
$a_{\mu}^{\text{SM}}=(11659181.0\pm 4.3)\times 10^{-10},$ resulting in the
$5.1\sigma$ discrepancy:

$\delta a_{\mu}^{\text{BSM}}=(24.9\pm 4.8)\times 10^{-10}.$ (1)

The SM prediction makes use of dispersion methods to calculate the hadronic
contributions that are the major source of uncertainty.111There is currently
no consensus between the dispersion method and lattice results reported by the
BMW collaboration BMW . On the one hand, if the latter is used, the $a_{\mu}$
discrepancy is reduced to $1.5\sigma$ level. On the other hand, it disagrees
with the R-ratio based prediction used in Aoyama:2020ynm  by more than
$2.1\sigma$. Comparison between different lattice groups is under way, with
consistent results among them lattice . Nevertheless, most of these results
stand for specific euclidean windows, not the full determination of hadronic
contributions to $a_{\mu}$. So an investigation of the source of
incompatibility among the two approaches is urgent, together with more high
precision lattice results. See Colangelo:2022vok  for an recent evaluation of
the different euclidean windows using the dispersion method. In this work, we
will adopt the result in eq. (1).

By disregarding neutrino mass generation and minimally introducing new physics
that couples to the muon in a flavor conserving fashion or obeying minimal
flavor violation, one can explain $a_{\mu}$ strumia ; freitas ; crivellin:g-2
; calibbi.ziegler  and still easily avoid large CLFV effects. A critical
review of the minimal extensions with one or two new fields can be seen in
Ref. stockinger ; see also Ref. queiroz . If the dominant new physics
contribution occurs at 1-loop, a key ingredient is the presence of a _chiral
enhancement_ that occurs by replacing one chiral flipping coupling involving
the muon Yukawa $y_{\mu}\sim 0.0006$ by an order one coupling crivellin:g-2 ;
calibbi.ziegler ; curtin ; strumia ; yin.yamaguchi .

With new physics at the TeV scale, the dipole operator responsible for
$a_{\mu}$ should be given by the dimension 6 operator
$\bar{\ell}_{\mu}\sigma^{\alpha\beta}H\mu_{R}F_{\alpha\beta}$ in the
electroweak symmetric phase where the necessary chirality flip
$\mu_{R}\to\mu_{L}$ involves a Higgs insertion. Then the one-loop contribution
to $a_{\mu}$ can be separated into two classes calibbi.ziegler ; curtin
depending on the location of the Higgs insertion: (i) if it occurs in the
external muon legs; (ii) if it occurs in one of the internal lines. For
contributions of type (i), there is no chiral enhancement because the Higgs
insertion is accompanied by the usual muon Yukawa $y_{\mu}$ which also
supplies a chiral flipping.222It is possible to explain $a_{\mu}$ with a type
(i) contribution but either the new states in the loop need to lie at the
electroweak scale or lower; see e.g. babu.lindner  or the couplings are at the
edge of perturbativity; see e.g. renato . If all new fields are electroweak
singlets, this is the only possibility. In contrast, chiral enhancement is
possible for contributions of type (ii) because the Higgs is attached to a
vertex other than the muon Yukawa. Enhancement occurs if this vertex is
replaced by the top Yukawa coupling (e.g., with the mediation of leptoquarks
lepto.quark ) or the tau Yukawa coupling (e.g., with lepton flavor changing
coupling with an additional Higgs doublet barman ). The large coupling may
also come from non-SM interactions such as with vectorlike leptons
crivellin:g-2 ; dermisek . In all these examples, the new Higgs insertion also
provides a chiral flipping.

With the presence of new scalars mixing with the SM Higgs, enhancement is also
possible through the couplings of these new scalars with the muon. In the Two-
Higgs-doublet model (2HDM), this feature is available for the types II/X or
aligned version with contributions enhanced by $\tan\beta$ (or, similarly,
$\zeta_{l}$ for the aligned version). However, in order to accommodate
$a_{\mu}$, some of the scalars must be at the weak scale, which renders only
type X 2hdm.typeX  or the aligned version Ilisie ; Han ; Cherchiglia1 ;
Cherchiglia2  as possible candidates. For the muon specific 2HDM mu-spec , it
is possible to push the scalars to higher masses. In the MSSM, a similar
reasoning apply, where the main contributions are also enhanced by
$\tan\beta$, allowing to accommodate $a_{\mu}$ (see e.g. g-2:susy  for a
review). However, differently from the 2HDM versions mentioned, in the MSSM
there are innumerous sources of flavor violation that need to be properly
suppressed endo.yin ; hisano ; see e.g. stockinger  for an updated analysis.

In this work, we aim to connect the physics responsible for $a_{\mu}$ with
neutrino masses generation so that the same mediators participate in both
processes. We focus on TeV scale mediators and we also seek minimality in
their number. Two general features immediately follow. Firstly, the seesaw
scale has to be lowered to the TeV scale or lower and, secondly, one certainly
cannot impose family lepton number to be conserved so that the interplay with
CLFV is crucial. Simply lowering the mediator masses for the type I typeI ,
type II typeII  or type III typeIII  seesaw models do not work because their
contribution to $a_{\mu}$ is negative freitas ; typeII:g-2 . 333Ref. ma.raidal
had considered the type II seesaw and its interplay with $a_{\mu}$ but they
did not consider the correct sign of the contribution.  In particular, adding
any number of singlet RHNs to implement the type I seesaw mechanism leads to a
negative contribution to $a_{\mu}$ coy.frigerio  even in low-scale seesaw
scenarios such as the inverse seesaw inverse.ss  or the symmetry protected
case kersten.smirnov . This situation does not change if the heavy particles
do not participate in the generation of neutrino masses.

Focusing on the seesaw, there are many ways that the seesaw scale can be
lowered. The aforementioned inverse seesaw is the case with extended tree
level mediators that brings additional suppression due to approximate lepton
number conservation. Another option is to consider radiative generation; see
Ref. cai  for a review and Refs. radiative:g-2 ; barman ; dicus.he  for models
connecting $a_{\mu}$ with radiative neutrino mass generation. Yet another way
is to attribute part of the smallness of neutrino masses to a vacuum
expectation value (VEV) of a second Higgs doublet that only couples to
neutrinos. This is the neutrinophilic 2HDM ($\nu$-2HDM) for which the RHNs can
be Majorana ma:nu2HDM ; ma.raidal  or Dirac nandi ; nu2hdm:U(1) . Our focus
will be on the former as the Majorana mass will induce the chirality flip
enhancing the $a_{\mu}$ contribution. An earlier attempt to connect neutrino
masses with $a_{\mu}$ considered this neutrinophilic 2HDM with Majorana RHNs
ma.raidal , where the interplay with CLFV was studied. However, in its
simplest form, the new contribution to $(g-2)_{\mu}$ were negative.
Phenomenological studies unrelated to $a_{\mu}$ were performed in Refs.
nu2HDM:pheno ; haba.yamada . A second Higgs doublet with small vev also helps
explaining the small neutrino mass in the recently proposed type IB seesaw
fu.king  where two Higgs doublets couple to the two RHNs that form a pseudo-
Dirac pair of an approximate $U(1)$ symmetry. Other recent works connecting
$a_{\mu}$ with neutrino mass generation can be seen in Ref. renato ;
nu.g-2:others .

Here we show that the shortcomings of the $\nu$-2HDM (and also the type IB
seesaw) to solve $a_{\mu}$ can be overcome by introducing another field —a
charged scalar singlet— to the model so that a chirally enhanced contribution
can be generated. This contribution has no definite sign and allows the new
fields to be at the TeV scale. Both the Majorana RHNs and the charged scalar
contributes to $a_{\mu}$ while the RHNs and the neutrinophilic Higgs doublets
are responsible for neutrino masses. In particular, the RHN Majorana mass in
the loop induces the necessarily chirality flip while the Higgs insertion
occurs with scalars and it is decoupled from the chiral flipping. Since lepton
flavor violation is built in with the TeV scale new physics, it is important
to study the induced CLFV processes.

The outline of the paper is as follows. In Sec. II we present the models and
in Sec. III we show how they generate the neutrino masses with TeV scale
mediators. Section IV shows the calculation for $a_{\mu}$ and CLFV
observables. Section V discusses the interplay between obtaining the correct
$a_{\mu}$ and avoiding the constraints of CLFV processes. The summary is given
in Sec. VI.

##  II Models

We seek scenarios where $a_{\mu}$ is explained with the participation of the
righthand neutrinos (RHNs) at one-loop at the same time that these heavy
neutrinos generate the necessary active neutrino masses through a low scale
seesaw mechanism. More specifically, we require a chiral enhanced contribution
to $a_{\mu}$ connecting $\mu_{L}$ with $\mu_{R}$ and then two more charged
scalars —one electroweak singlet and another residing in a doublet— are
necessary to close the loop as shown in Fig. 1. From the diagram it is clear
that the chiral enhancement will be proportional to the RHN Majorana mass
$M_{N}$ which we require to be at the TeV scale.444This is one of the minimal
possibilities considered in Ref. calibbi.ziegler  to explain $(g-2)_{\mu}$ and
DM stabilized by an unbroken $\mathbb{Z}_{2}$ leading to an inert Higgs
doublet. Here the additional Higgs doublet has a small vev that helps
explaining small neutrino masses.

Figure 1: Chirally enhanced contribution to muon $g-2$ involving RHN $N_{R}$
and charged scalars.

So we extend the SM by considering two Higgs doublets $\Phi_{1},\Phi_{2}$, two
555These two neutrinos will generate the minimal number of two nonzero light
neutrino masses but three RHNs can be equally considered. righthanded neutrino
(RHN) fields $N_{iR}$, $i=1,2$, and one singlet charged scalar $\varphi^{+}$.
The relevant part of the Lagrangian involving the Higgs doublets and the
righthanded neutrino fields is

$-\mathscr{L}\supset\bar{\ell}_{\alpha}h_{\alpha}\Phi_{e}e_{\alpha
R}+\bar{N}_{iR}\lambda^{(1)}_{i\alpha}\tilde{\Phi}_{1}^{\dagger}\ell_{\alpha}+\bar{N}_{iR}\lambda^{(2)}_{i\alpha}\tilde{\Phi}_{2}^{\dagger}\ell_{\alpha}+\mbox{\large$\tfrac{1}{2}$}\bar{N}_{iR}(M_{R})_{ij}N_{jR}^{c}+h.c.,$
(2)

where $\ell_{\alpha}$, $\alpha=e,\mu,\tau$ are the lepton doublets and
$\Phi_{e}$ is $\Phi_{1}$ or $\Phi_{2}$ depending on the model. The relevant
interaction terms involving the charged singlet scalar is

$-\mathscr{L}\supset\mu_{\varphi}\Phi_{2}^{\mbox{\scriptsize$\mathsf{T}$}}\epsilon\Phi_{1}\varphi^{-}+f_{i\alpha}\bar{N}_{iR}e^{c}_{\alpha
R}\varphi^{-}+h.c.,$ (3)

where we choose $\mu_{\varphi}$ real without loss of generality. The
interaction term $\overline{\ell^{c}_{\alpha}}\ell_{\beta}\varphi^{+}$ is
absent owing to the symmetries of the specific models introduced below and
then there is no radiative neutrino mass generation as in the Zee model
zee.model .666A recent paper considered the Zee model for $(g-2)_{e,\mu}$
barman , where the couplings to charged singlet are not relevant and $a_{\mu}$
is explained by the contributions of the additional Higgs doublet with a
chiral enhancement at one-loop of at most the $\tau$ Yukawa. Ref. dicus.he
considered the complementary case of two charged scalar singlets and one RHN
where neutrino masses are radiative but there is no chirally enhanced
contribution to $a_{\mu}$.  This absence also forbids the generation of four-
fermion operators $\ell^{4}$ from tree level $\varphi^{+}$ exchange.

One of the simplest models within the type I seesaw is to suppress the Dirac
mass term by attributing its origin to a Higgs doublet different from the
rest, with tiny vacuum expectation value. This is the neutrinophilic 2HDM
model ($\nu$-2HDM) ma:nu2HDM  which is obtained from (2) by imposing a
$\mathbb{Z}_{2}$ where $N_{iR},\Phi_{1}$ are odd. As a consequence,
$\lambda^{(2)}=0$ and $\Phi_{e}=\Phi_{2}$. The singlet scalar $\varphi^{+}$ is
additional and we assume it is odd so that the terms in (3) are allowed. The
resulting Lagrangian is

$\displaystyle-\mathscr{L}_{\nu\text{-2HDM}}$
$\displaystyle\supset\bar{\ell}_{\alpha}h_{\alpha}\Phi_{2}e_{\alpha
R}+\bar{N}_{iR}\lambda^{(1)}_{i\alpha}\tilde{\Phi}_{1}^{\dagger}\ell_{\alpha}+\mbox{\large$\tfrac{1}{2}$}\bar{N}_{iR}M_{N_{i}}N_{iR}^{c}$
(4) $\displaystyle\quad+\
\mu_{\varphi}\Phi_{2}^{\mbox{\scriptsize$\mathsf{T}$}}\epsilon\Phi_{1}\varphi^{-}+f_{i\alpha}\bar{N}_{iR}e^{c}_{\alpha
R}\varphi^{-}+h.c.$

The doublet $\Phi_{2}\approx H_{\rm SM}$ will be mostly the SM Higgs doublet
while $\Phi_{1}\approx H_{\nu}$ will be mostly composed of non-SM higgses.
Neutrino masses will depend solely on $\langle\Phi_{1}^{0}\rangle=v_{1}$ which
will be suppressed.

Another minimal possibility is to combine the righthanded neutrinos into a
(pseudo)Dirac pair through a $U(1)_{N}$ symmetry with charges fu.king

$N_{1R}\sim\Phi_{1}\sim+1\,,\quad N_{2R}\sim\Phi_{2}\sim-1\,.$ (5)

It leads to the Lagrangian

$-\mathscr{L}_{\rm
Ib}\supset\bar{N}_{R1}\lambda_{1\alpha}\tilde{\Phi}_{1}^{\dagger}\ell_{\alpha}+\bar{N}_{R2}\lambda_{2\alpha}\tilde{\Phi}_{2}^{\dagger}\ell_{\alpha}+M\overline{N_{R1}^{c}}N_{R2}+h.c.$
(6)

and the seesaw mechanism is dubbed seesaw type IB. Compared to (2), we are
already simplifying the notation in that
$\lambda^{(1)}_{1\alpha}=\lambda_{1\alpha}$ and
$\lambda^{(2)}_{2\alpha}=\lambda_{2\alpha}$ whereas
$\lambda^{(1)}_{2\alpha}=\lambda^{(2)}_{1\alpha}=0$. The other couplings
depend on the choice of $U(1)_{N}$ charges for $e_{R\alpha}$ and $\varphi$. We
adopt the charges $e_{R}\sim-1$ and $\varphi^{-}\sim-2$ so that

$-\mathscr{L}_{\rm Ib}\supset
h_{\alpha}\bar{\ell}_{\alpha}\Phi_{1}e_{R\alpha}+f_{2\alpha}\bar{N}_{2R}e^{c}_{\alpha
R}\varphi^{-}+\
\mu_{\varphi}\Phi_{2}^{\mbox{\scriptsize$\mathsf{T}$}}\epsilon\Phi_{1}\varphi^{-}+h.c.$
(7)

To keep $h_{\tau}$ within perturbative values, we need $v_{1}\gtrsim 10^{-2}v$
and most of the suppression for light neutrino masses should still come from
the suppression of the Yukawas $\lambda_{j\alpha}$. We could flip the charges
of $e_{R}$ and $\varphi^{-}$ and then replace $\Phi_{1}$ by $\Phi_{2}$ in the
first term and $N_{2R}$ by $N_{1R}$ in the second. But even in this case we
cannot suppress $v_{1}$ too much because the dominant $g-2$ contribution is
proportional to $v_{1}\lambda_{2\alpha}$. The third term softly breaks the
$U(1)_{N}$ symmetry 777 The quadratic term $\Phi_{2}^{\dagger}\Phi_{1}$ should
be added as well. The term $\ell\ell\varphi^{+}$ remains forbidden.  which
cannot be exact as a global symmetry to avoid unwanted massless scalars. The
model in Ref. (fu.king, , a), for example, imposes a $\mathbb{Z}_{3}$ instead
of $U(1)_{N}$.

For both models, the charged scalar component that is dominant in the non-SM
Higgs doublet behaves as the charged scalar in the type I 2HDM and they are
only constrained by LEP ALEPH:2013htx :

$M_{H^{+}}>75\,\mathrm{GeV}\,.$ (8)

On the other hand, the singlet component couples solely with charged leptons
and righthand neutrino, implying that the model independent constraint from
LEP does not apply. Since we allow for a small mixing, we will consider that
both charged scalars will have masses around the electroweak scale or above.

##  III Light neutrino masses

Light neutrino masses are generated from the type I seesaw, i.e., through the
exchange of RHNs in (2) at tree level. The generated effective Weinberg
operator is

$\mathscr{L}=\frac{1}{2}\bar{\ell}^{c}_{\alpha}\Gamma^{\mbox{\scriptsize$\mathsf{T}$}}_{\alpha
i}(M_{R}^{-1})_{ij}\Gamma_{j\beta}\ell_{\beta}+h.c.\,,$ (9)

depending on the two Higgs doublets through

$\Gamma_{j\beta}=\lambda^{(1)}_{j\beta}\tilde{\Phi}_{1}^{\dagger}+\lambda^{(2)}_{j\beta}\tilde{\Phi}_{2}^{\dagger}\,.$
(10)

###  III.1 $\nu$-2HDM model

For the $\nu$-2HDM model, only $\Phi_{1}$ participates in the seesaw and
$\lambda^{(2)}=0$. Then the neutrino mass matrix is given by

$M_{\nu}=-v_{1}^{2}{\lambda^{(1)}}^{\mbox{\scriptsize$\mathsf{T}$}}M_{R}^{-1}\lambda^{(1)}\,,$
(11)

where its lightness is partly explained by
$\langle\Phi_{1}^{0}\rangle=v_{1}\ll
v=\sqrt{v_{1}^{2}+v_{2}^{2}}=174\,\mathrm{GeV}$. If $M_{R}\sim\mathrm{TeV}$,
we need

$v_{1}\lambda^{(1)}/v\sim\frac{\sqrt{M_{R}M_{\nu}}}{v}\lesssim\frac{\sqrt{1\,\mathrm{TeV}\times\,0.05\,\mathrm{eV}}}{v}\sim
10^{-6}\,.$ (12)

We can choose $v_{1}/v\sim 10^{-6}$ in order to have $\lambda^{(1)}$ of order
one to address the deviation in $a_{\mu}$.

Since we are focusing on the minimal case of two RHNs, the heavy neutrino
Yukawa coupling $\lambda^{(1)}$ is mostly fixed by the masses and mixing of
light neutrinos. Using the Casas-Ibarra parametrization, we can write for
normal ordering (NO),

$\displaystyle v_{1}\lambda^{(1)}_{1\alpha}$
$\displaystyle=i\sqrt{M_{1}}\big{(}\sqrt{m_{2}}c_{z}V^{\dagger}_{2\alpha}-\sqrt{m_{3}}s_{z}V^{\dagger}_{3\alpha}\big{)}\,,$
(13) $\displaystyle v_{1}\lambda^{(1)}_{2\alpha}$
$\displaystyle=i\sqrt{M_{2}}\big{(}\sqrt{m_{2}}s_{z}V^{\dagger}_{2\alpha}+\sqrt{m_{3}}c_{z}V^{\dagger}_{3\alpha}\big{)}\,,$

where $c_{z}=\cos z$ and $s_{z}=\sin z$ depend on the free complex angle $z$.
We are taking $M_{R}=\operatorname{\mathrm{diag}}(M_{1},M_{2})$ and the
neutrino mass matrix in (11) is
$M_{\nu}=V^{*}\operatorname{\mathrm{diag}}(0,m_{2},m_{3})V^{\dagger}$, with
$V$ being the PMNS matrix. Therefore, besides the two CP phases in $V$, we
have five free parameters in the neutrino sector:
$M_{1},M_{2},\operatorname{\mathrm{Re}}(z),\operatorname{\mathrm{Im}}(z),\tan\beta=v_{2}/v_{1}$.
Note that for $M_{1}=M_{2}$, the real part of $z$ is not physical.

For inverted ordering (IO), we have instead

$\displaystyle v_{1}\lambda^{(1)}_{1\alpha}$
$\displaystyle=i\sqrt{M_{1}}\big{(}\sqrt{m_{1}}c_{z}V^{\dagger}_{1\alpha}-\sqrt{m_{2}}s_{z}V^{\dagger}_{2\alpha}\big{)}\,,$
(14) $\displaystyle v_{1}\lambda^{(1)}_{2\alpha}$
$\displaystyle=i\sqrt{M_{2}}\big{(}\sqrt{m_{1}}s_{z}V^{\dagger}_{1\alpha}+\sqrt{m_{2}}c_{z}V^{\dagger}_{2\alpha}\big{)}\,.$

The only difference in this case is that
$M_{\nu}=V^{*}\operatorname{\mathrm{diag}}(m_{1},m_{2},0)V^{\dagger}$.

###  III.2 Type Ib seesaw model

For the type IB seesaw model, the neutrino mass matrix coming from (9) is

$(M_{\nu})_{\alpha\beta}=-\frac{v_{1}v_{2}}{M}(\lambda_{1\alpha}\lambda_{2\beta}+\lambda_{2\alpha}\lambda_{1\beta})\,,$
(15)

where $v_{a}=\langle\Phi_{a}^{0}\rangle$, $a=1,2$, and the Yukawa couplings
are defined in (6). We choose $v_{1}<v_{2}$ because $\Phi_{2}$ is the doublet
coupling to the charged leptons in (7). Adapting (12) to this case with
$M\sim\mathrm{TeV}$, we obtain

$\frac{1}{v}\sqrt{v_{1}\lambda_{1\alpha}v_{2}\lambda_{2\beta}}\sim\frac{\sqrt{M_{R}M_{\nu}}}{v}\lesssim\frac{\sqrt{1\,\mathrm{TeV}\times
0.05\,\mathrm{eV}}}{v}\sim 10^{-6}\,.$ (16)

As in the $\nu$-2HDM model, the Yukawa couplings $\lambda_{1\alpha}$ and
$\lambda_{2\alpha}$ are mostly fixed by light neutrino masses and mixing.

For NO, the Yukawa couplings can be parametrized as

$\displaystyle\lambda_{1\alpha}$
$\displaystyle=i\kappa\left(\frac{M}{2v_{1}v_{2}}\right)^{1/2}\bigg{[}+i\sqrt{m_{2}}\,V^{\dagger}_{2\alpha}+\sqrt{m_{3}}\,V^{\dagger}_{3\alpha}\bigg{]}\,,$
(17) $\displaystyle\lambda_{2\alpha}$
$\displaystyle=i\kappa^{-1}\left(\frac{M}{2v_{1}v_{2}}\right)^{1/2}\bigg{[}-i\sqrt{m_{2}}\,V^{\dagger}_{2\alpha}+\sqrt{m_{3}}\,V^{\dagger}_{3\alpha}\bigg{]}\,,$

where we use a different but equivalent parametrization with respect to Ref.
fu.king . Flipping the sign in front of $\sqrt{m_{2}}$ is equivalent and this
sign can be absorbed by flipping the sign of $\nu_{2L}$. The parameter
$\kappa$ is free and can be chosen real by rephasing $N_{1R},N_{2R}$ with
opposite phases. For IO, we can analogously parametrize

$\displaystyle\lambda_{1\alpha}$
$\displaystyle=i\kappa\left(\frac{M}{2v_{1}v_{2}}\right)^{1/2}\bigg{[}+i\sqrt{m_{1}}\,V^{\dagger}_{1\alpha}+\sqrt{m_{2}}\,V^{\dagger}_{2\alpha}\bigg{]}\,,$
(18) $\displaystyle\lambda_{2\alpha}$
$\displaystyle=i\kappa^{-1}\left(\frac{M}{2v_{1}v_{2}}\right)^{1/2}\bigg{[}-i\sqrt{m_{1}}\,V^{\dagger}_{1\alpha}+\sqrt{m_{2}}\,V^{\dagger}_{2\alpha}\bigg{]}\,.$

Note that in the neutrino sector, besides the unknown CP phases, there are
only three free parameters: $M,\tan\beta,\kappa$. This is the same number as
the neutrinophilic case with equal masses for the RHNs.

##  IV Dipole moments and CLFV

Using an effective theory approach, the operators relevant to lepton dipole
moments and charged lepton flavor violation (CLFV) processes are the photonic
operators

$\displaystyle\mathscr{L}_{\gamma\text{-eff}}$ $\displaystyle=-\left(C^{\sigma
R}_{\alpha\beta}\bar{e}_{\alpha L}\sigma_{\mu\nu}e_{\beta
R}F^{\mu\nu}+h.c.\right)$ (19) $\displaystyle\quad-\
\left(C_{\alpha\beta}^{\text{ND-}L}\bar{e}_{\alpha L}\gamma_{\nu}e_{\beta
L}+C_{\alpha\beta}^{\text{ND-}R}\bar{e}_{\alpha R}\gamma_{\nu}e_{\beta
R}\right)\partial_{\mu}F^{\mu\nu}\,.$

The operator in the first line is the dipole contribution whereas the ones in
the second line are the non-dipole (ND) part. We can see that only the dipole
part involves chirarity flipping and chiral enhancement will be possible only
for this term. The Wilson coefficients at 1-loop can be obtained by matching
the full theory with the effective theory through appropriate 1-loop
amplitudes. The relevant amplitudes come from the dipole (left) and self-
energy (right) diagrams in Fig. 2.

Figure 2: Dipole and self-energy contribution to flavor changing processes.

###  IV.1 Dipole moments and $\ell_{\alpha}\to\ell_{\beta}\gamma$

The Wilson coefficient of the dipole operator contributes to the dipole
moments and $\ell_{\alpha}\to\ell_{\beta}\gamma$ as

$\displaystyle a_{\alpha}$
$\displaystyle=-\frac{4m_{\alpha}}{e}\operatorname{\mathrm{Re}}(C^{\sigma
R}_{\alpha\alpha})\,,$ (20) $\displaystyle d_{\alpha}$
$\displaystyle=-2\operatorname{\mathrm{Im}}(C^{\sigma R}_{\alpha\alpha})\,,$
$\displaystyle\operatorname{Br}[\ell_{\alpha}\to\ell_{\beta}\gamma]$
$\displaystyle=\frac{m^{3}_{\alpha}}{4\pi\Gamma_{\alpha}}(|C^{\sigma
R}_{\alpha\beta}|^{2}+|C^{\sigma R}_{\beta\alpha}|^{2})\,,$

where $a_{\alpha},d_{\alpha}$ are the contributions to the magnetic and
eletric dipole moments, respectively. These formulas assume that our covariant
derivative in QED is $D_{\mu}=\partial_{\mu}+ieQA_{\mu}$ and similarly for the
SM.

For the $\nu$-2HDM described in (4), $\Phi_{2}\approx H_{\rm SM}$ and
$\Phi_{1}\approx H_{\nu}$ is almost inert with a tiny vev. The chirally
enhanced contribution to the dipole operator is given in Fig. 1, with
$\varphi$ being the charged singlet, $H^{+}$ being essentially the charged
Higgs residing in $\Phi_{1}$ and $h^{0}$ being the SM Higgs boson. The chiral
enhancement comes from the $N_{iR}$ Majorana masses, $M_{N_{i}}$, and these
heavy neutrinos only couple to $\Phi_{1}$. The expected chirally enhanced
contribution is then of the form

$C^{\sigma
R}_{\beta\alpha}\sim\frac{ev_{2}\mu_{\varphi}}{(4\pi)^{2}\Lambda^{4}}{\lambda^{(1)}_{\beta
i}}^{\dagger}M_{N_{i}}f^{*}_{i\alpha}\,,$ (21)

where the transferred momentum is $q^{2}=-m_{\mu}^{2}$.

For the exact contribution we need to include the mixing between the charged
higgs $H^{+}$ in the doublet and the scalar singlet $\varphi^{+}$ coming from
the $\mu_{\varphi}$ term in (4) or (3). The contribution in the scalar
potential is

$\mu_{\varphi}\varphi^{-}(\Phi_{2}^{+}\Phi_{1}^{0}-\Phi_{2}^{0}\Phi_{1}^{+})\to-\mu_{\varphi}v\,\varphi^{-}H^{+}\,,$
(22)

where $vH^{+}\equiv v_{2}\Phi_{1}^{+}-v_{1}\Phi_{2}^{+}$,
$v_{i}=\langle\Phi_{i}^{0}\rangle$, is the physical charged Higgs field within
the doublets when $\mu_{\varphi}=0$. The orthogonal direction is the charged
Goldstone absorbed by the $W$. Considering arbitrary mass terms for
$H^{+}H^{-}$ and $\varphi^{+}\varphi^{-}$, the mixing term (22) induces the
mixing

$\begin{pmatrix}H^{+}\cr\varphi^{+}\end{pmatrix}=\begin{pmatrix}c_{\gamma}&-s_{\gamma}\cr
s_{\gamma}&c_{\gamma}\end{pmatrix}\begin{pmatrix}S_{1}^{+}\cr
S_{2}^{+}\end{pmatrix}\,,$ (23)

with angle

$\sin 2\gamma=-\frac{2\mu_{\varphi}v}{M^{2}_{S_{2}}-M^{2}_{S_{1}}}\,,$ (24)

and $S_{i}^{+}$ are the charged scalars with masses $M_{S_{i}}$, $i=1,2$. We
have chosen $M_{S_{2}}>M_{S_{1}}$.

Then the complete contribution arising from the diagrams in Fig. 2 is

$\displaystyle\frac{16\pi^{2}}{e}C^{\sigma R}_{\beta\alpha}$
$\displaystyle=c_{\gamma}s_{\gamma}\frac{v_{2}}{v}\sum_{j}\frac{{\lambda^{(1)}_{\beta
j}}^{\dagger}f^{*}_{j\alpha}}{M_{N_{j}}}\left[x_{2j}f_{S}(x_{2j})-x_{1j}f_{S}(x_{1j})\right]\,,$
(25a) $\displaystyle\quad+\
m_{\alpha}\frac{v_{2}^{2}}{v^{2}}\sum_{j}\lambda^{(1){\dagger}}_{\beta
j}\left[\frac{c_{\gamma}^{2}}{M^{2}_{S_{1}}}\tilde{f}_{S}(x_{1j})+\frac{s_{\gamma}^{2}}{M^{2}_{S_{2}}}\tilde{f}_{S}(x_{2j})\right]\lambda^{(1)}_{j\alpha}$
(25b) $\displaystyle~{}~{}\
+m_{\beta}\sum_{j}f^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
j}\left[\frac{s_{\gamma}^{2}}{M^{2}_{S_{1}}}\tilde{f}_{S}(x_{1j})+\frac{c_{\gamma}^{2}}{M^{2}_{S_{2}}}\tilde{f}_{S}(x_{2j})\right]f_{j\alpha}^{*}\,,$
(25c)

where $x_{kj}\equiv M^{2}_{N_{j}}/M^{2}_{S_{k}}$ and the loop functions are
crivellin:g-2 ; calibbi.ziegler

$\displaystyle f_{S}(x)$ $\displaystyle\equiv\frac{x^{2}-1-2x\log
x}{4(x-1)^{3}}\,,$ (26) $\displaystyle\tilde{f}_{S}(x)$
$\displaystyle\equiv\frac{2x^{3}+3x^{2}-6x+1-6x^{2}\log x}{24(x-1)^{4}}\,.$

See appendix A for the details. We have checked that the contribution
proportional to $(\lambda^{(1)})^{2}$ matches Ref. haba.yamada  for
$\gamma=0$.

Let us discuss the various contributions. The first contribution (25a) is the
chirally enhanced (left-right), reducing to (21) when $M_{S_{i}}\gg M_{N_{j}}$
as $f_{S}(0)=1/4$ and we identify
$\Lambda^{4}=4M_{S_{2}}^{2}M_{S_{1}}^{2}(M_{S_{2}}^{2}-M_{S_{1}}^{2})$. The
contributions (25b) and (25c) are not chirally enhanced (left-left and right-
right respectively) as the chiral flip comes from the external lines. So the
chirally enhanced contribution is larger than the non-chirally enhanced
contributions by a factor

$c_{\gamma}s_{\gamma}\frac{M^{2}_{S}}{m_{\mu}M_{N}}\sim
c_{\gamma}s_{\gamma}\times 10^{4}\,,$ (27)

for order one couplings and the numbers assume $M_{N}\sim
M_{S_{i}}\sim\mathrm{TeV}$. Given the definition of the angle $\gamma$ in (24)
the trilinear coupling $\mu_{\varphi}$ cannot be arbitrarily small if we
require the chirally enhanced contribution to dominate:

$|\mu_{\varphi}|\gtrsim
1\,\mathrm{GeV}\text{~{}~{}or~{}~{}}|s_{\gamma}|\gtrsim 10^{-4}\,.$ (28)

The dominance of the chirally enhanced contribution is important because the
non-chirally enhanced contributions are positive definite and leads to a
contribution to $a_{\mu}$ which is negative definite, contrary to the
experimental observation. They are similar to the contribution of $N_{R}$
exchange in the usual seesaw models freitas ; coy.frigerio .

For the type IB model in (6) and (7), the estimate of the chirally enhanced
contribution analogous to (21) will be

$C^{\sigma
R}_{\beta\alpha}\sim\frac{v_{2}\mu_{\varphi}}{(4\pi)^{2}\Lambda^{4}}{\lambda}^{\dagger}_{\beta
1}Mf^{*}_{2\alpha}\,.$ (29)

Then we need $\lambda_{1\beta}$ to be of order one while $\lambda_{2\beta}$
needs to be suppressed due to neutrino masses, cf. (16).

The complete contribution reads

$\displaystyle\frac{16\pi^{2}}{e}C^{\sigma R}_{\beta\alpha}$
$\displaystyle=c_{\gamma}s_{\gamma}\frac{v_{2}}{v}\frac{\lambda^{\dagger}_{\beta
1}f^{*}_{2\alpha}}{M_{N}}\left[x_{2}f_{S}(x_{2})-x_{1}f_{S}(x_{1})\right]\,,$
(30a) $\displaystyle\quad+\
m_{\alpha}\frac{v_{2}^{2}}{v^{2}}\lambda^{{\dagger}}_{\beta
1}\left[\frac{c_{\gamma}^{2}}{M^{2}_{S_{1}}}\tilde{f}_{S}(x_{1})+\frac{s_{\gamma}^{2}}{M^{2}_{S_{2}}}\tilde{f}_{S}(x_{2})\right]\lambda_{1\alpha}$
(30b) $\displaystyle~{}~{}\
+m_{\beta}f^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
2}\left[\frac{s_{\gamma}^{2}}{M^{2}_{S_{1}}}\tilde{f}_{S}(x_{1})+\frac{c_{\gamma}^{2}}{M^{2}_{S_{2}}}\tilde{f}_{S}(x_{2})\right]f_{2\alpha}^{*}\,,$
(30c)

where $x_{i}=M^{2}_{N}/M_{S_{i}^{2}}$. We have neglected additional
contributions proportional to $\lambda_{2\beta}^{*}\lambda_{2\alpha}v_{1}^{2}$
which will be highly suppressed. Here the requirement (28) is equally
necessary for order one couplings.

In Table 1 we show the current and future limits for different
$\ell_{\alpha}\to\ell_{\beta}\gamma$. Limits for other CLFV processes are also
shown.

$\begin{array}[]{|c|c|c|c|}\hline\cr\text{Observable}&\text{Current
limit}&\text{Future limit}\cr\hline\cr\text{Br}(\mu\to eee)&<1.0\times
10^{-12}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{PDG2020}{\@@citephrase{(}}{\@@citephrase{)}}}}&10^{-16}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{mu3e}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\text{Br}(\tau\to\mu\mu\mu)&<2.1\times
10^{-8}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{HFLAV:2019otj}{\@@citephrase{(}}{\@@citephrase{)}}}}&3.4\times
10^{-10}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{belle2:book}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\text{Br}(\tau\to\mu
ee)&<8.4\times 10^{-9}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{HFLAV:2019otj}{\@@citephrase{(}}{\@@citephrase{)}}}}&2.9\times
10^{-10}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{belle2:book}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\text{Br}(\tau\to
eee)&<1.4\times 10^{-8}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{HFLAV:2019otj}{\@@citephrase{(}}{\@@citephrase{)}}}}&4.3\times
10^{-10}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{belle2:book}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\text{Br}(\tau\to
e\mu\mu)&<1.6\times 10^{-8}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{HFLAV:2019otj}{\@@citephrase{(}}{\@@citephrase{)}}}}&4.3\times
10^{-10}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{belle2:book}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\hline\cr\text{Br}(\mu\to
e\gamma)&<4.2\times 10^{-13}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{PDG2020}{\@@citephrase{(}}{\@@citephrase{)}}}}&6\times
10^{-14}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{meg2}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\text{Br}(\tau\to\mu\gamma)&<4.4\times
10^{-8}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{PDG2020}{\@@citephrase{(}}{\@@citephrase{)}}}}&10^{-9}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{belle2:book}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\text{Br}(\tau\to
e\gamma)&<3.3\times 10^{-8}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{PDG2020}{\@@citephrase{(}}{\@@citephrase{)}}}}&3\times
10^{-9}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{belle2:book}{\@@citephrase{(}}{\@@citephrase{)}}}}\cr\hline\cr\Gamma^{\rm
conv}_{\mu\rightarrow e}/\Gamma^{\rm capt}_{N}&<7.0\times
10^{-13}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{SINDRUMII:2006dvw}{\@@citephrase{(}}{\@@citephrase{)}}}}^{*}&3\times
10^{-17}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{mu2e,comet}{\@@citephrase{(}}{\@@citephrase{)}}}}^{**}/10^{-18}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{jparc}{\@@citephrase{(}}{\@@citephrase{)}}}}^{\dagger}\cr\hline\cr\end{array}$
Table 1:  Current and future limits for charged lepton flavor violating
processes at 90% CL. For $\mu e$ conversion, the nucleous is Au (*) and Al
(**) or Ti ($\dagger$).

###  IV.2 $\mu$-$e$ conversion in nuclei

A very stringent test for CFLV is the coherent $\mu^{-}\text{\,--\,}e^{-}$
conversion in a muonic atom of nucleus $(A,Z)$ by neutrinoless muon capture

$\mu^{-}+(A,Z)\to e^{-}+(A,Z)\,,$

mediated in our case by the effective flavor-changing photon interactions in
(19). We neglect similar $Z$ mediated processes which are suppressed by the
$Z$ mass. It should be noticed that, contrary to $\mu\rightarrow e\gamma$
decay, one needs to consider off-shell photon emission. The general photonic
$\mu{-}e$ transition amplitude is given by Kuno

$\mathcal{M}=-eA_{\mu}^{*}(q)\bar{u}_{e}(p_{e})\left[(f_{E0}+\gamma_{5}f_{M0})\gamma_{\nu}\left(g^{\mu\nu}-\frac{q^{\mu}q^{\nu}}{q^{2}}\right)+(f_{M1}+\gamma_{5}f_{E1})\frac{i\sigma_{\mu\nu}q^{\nu}}{m_{\mu}}\right]u_{\mu}(p_{\mu})\,,$
(31)

where $p_{\mu}$ and $p_{e}$ are the momenta of the muon and electron
respectively, $q=p_{\mu}-p_{e}$ is the transfered momentum, and $f_{X}\equiv
f_{X}(q^{2})$ are form factors. In the case of $\mu\rightarrow e\gamma$ decay,
only the dipole part ($f_{E1}$, $f_{M1}$) contributes. The form factors can be
written in terms of the Wilson coefficients introduced in (19):

$\displaystyle f_{E0}$
$\displaystyle=(C_{e\mu}^{\text{ND-}L}+C_{e\mu}^{\text{ND-}R})\frac{q^{2}}{2e},\quad
f_{M0}=(C_{e\mu}^{\text{ND-}L}-C_{e\mu}^{\text{ND-}R})\frac{q^{2}}{2e},$ (32)
$\displaystyle f_{E1}$ $\displaystyle=-(C_{e\mu}^{\sigma R}-C_{\mu e}^{\sigma
R*})\frac{m_{\mu}}{e},\quad f_{M1}=-(C_{e\mu}^{\sigma R}+C_{\mu e}^{\sigma
R*})\frac{m_{\mu}}{e}\,.$ (33)

In the approximation pioneered by Weinberg and Feinberg Weinberg:1959zz , only
the photonic contribution is included in the calculation of the branching
ratio of the coherent $\mu-e$ conversion which is given by

$\operatorname{Br}[\mu N\to eN]\equiv\frac{\Gamma[\mu N\to eN]}{\Gamma_{\rm
capt}}=\frac{8m_{\mu}\alpha^{5}Z_{\rm eff}^{4}Z|F_{p}|^{2}\xi^{2}}{\Gamma_{\rm
capt}}\,.$ (34)

Here $Z_{\rm eff}$ is an effective atomic charge due to averaging the muon
wave function over the nuclear density, $\Gamma_{\rm capt}$ is the total muon
capture rate, and $\xi^{2}$ is obtained with knowledge of the form factors

$\xi^{2}=|f_{E0}(-m_{\mu}^{2})+f_{M1}(-m_{\mu}^{2})|^{2}+|f_{E1}(-m_{\mu}^{2})+f_{M0}(-m_{\mu}^{2})|^{2}\,,$
(35)

where we use the transferred momentum $q^{2}=-m_{\mu}^{2}$.

In general, the non-photonic part should be included as well, which is
parametrized by an effective four-fermion interaction containing quarks. In
our particular case, these contributions would come from box diagrams
containing both the righthanded neutrinos and charged scalars. Since the
coupling between quarks and the charged scalars is suppressed (it will be
proportional to the ratio $v_{1}/v_{2}$), we can safely neglect these
diagrams.

Finally, we rewrite the branching ratio in terms of the Wilson coefficients,
to allow a straightforward application to our model

$\operatorname{Br}[\mu N\to eN]=\frac{4m^{3}_{\mu}\alpha^{4}Z_{\rm
eff}^{4}Z|F_{p}|^{2}}{\pi\Gamma_{\rm capt}}\left(\left|C^{\sigma
R}_{e\mu}+\frac{m_{\mu}}{2}C^{\rm ND-L}_{e\mu}\right|^{2}+\left|C^{\sigma
R}_{\mu e}+\frac{m_{\mu}}{2}(C^{\rm ND-R}_{e\mu})^{*}\right|^{2}\right)\,.$
(36)

In our case, for order one couplings and loop functions, the chiral enhanced
part in the dipole contribution will dominate the non-dipole part by a factor
$M_{S_{i}}/m_{\mu}\sim M_{N_{k}}/m_{\mu}\sim 10^{4}$ for the new fields in the
TeV scale and $s_{\gamma}$ not too small as (28). In this case of dipole
dominance, the ratio $\operatorname{Br}[\mu N\to eN]/\operatorname{Br}(\mu\to
e\gamma)$ depends roughly only on the atomic number of the nucleous cirigliano
and we have checked it for our case. As for the non-dipole part, we have found
that the relative sign between $C^{\sigma R}$ and $C^{\rm ND}$ in (36) is
effectively opposite to Refs. haba.yamada ; bertuzzo.perez ; arganda .

At present, the best limit for this branching ratio comes from conversion in
gold nuclei SINDRUMII:2006dvw : $\Gamma_{\rm conv}^{\rm Au}/\Gamma_{\rm
capt}^{\rm Au}<7.0\times 10^{-13}$, where $\Gamma_{\rm capt}^{\rm
Au}=8.7\times 10^{-18}$ GeV. There are, however, future experiments that aim
to reduce the bounds on $\mu-e$ conversion by several orders of magnitude mu2e
; comet . Their aim is to achieve $\Gamma_{\rm conv}^{\rm Al}/\Gamma_{\rm
capt}^{\rm Al}<3\times 10^{-17}$ using aluminium nuclei. The different
effective parameters for different nucleous are shown in Table 2.

$\begin{array}[]{|c|c|c|c|}\hline\cr\ce{{}^{$A$}_{$Z$}Nucleus}&Z_{\rm
eff}&F_{p}&\Gamma_{\rm
capt}(\mathrm{GeV})\cr\hline\cr\ce{{}^{27}_{13}Al}&11.5&0.64&4.64\times
10^{-19}\cr\ce{{}^{48}_{22}Ti}&17.6&0.54&1.70\times
10^{-18}\cr\ce{{}^{80}_{38}Sr}&25.0&0.39&4.62\times
10^{-18}\cr\ce{{}^{121}_{51}Sb}&29.0&0.32&6.72\times
10^{-18}\cr\ce{{}^{197}_{79}Au}&33.5&0.16&8.60\times
10^{-18}\cr\ce{{}^{207}_{82}Pb}&\;34.0&\;0.15&8.85\times
10^{-18}\cr\hline\cr\end{array}$ Table 2: Values for $Z_{\rm eff}$, $F_{p}$,
and $\Gamma_{\rm capt}(\mathrm{GeV})$ for different nuclei Kitano ; arganda .

The Wilson coefficient of the non-dipole part can be obtained for the
$\nu$-2HDM model as

$\displaystyle\frac{16\pi^{2}}{e}C_{\beta\alpha}^{\text{ND-}L}$
$\displaystyle=\frac{v_{2}^{2}}{v^{2}}\sum_{j}{\lambda^{(1)}_{\beta
j}}^{\dagger}\lambda^{(1)}_{j\alpha}\left[c_{\gamma}^{2}\frac{G_{S}(x_{1j})}{6M^{2}_{S_{1}}}+s_{\gamma}^{2}\frac{G_{S}(x_{2j})}{6M^{2}_{S_{2}}}\right]\,,$
(37) $\displaystyle\frac{16\pi^{2}}{e}C_{\beta\alpha}^{\text{ND-}R}$
$\displaystyle=\sum_{j}f^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
j}f^{*}_{j\alpha}\left[s_{\gamma}^{2}\frac{G_{S}(x_{1j})}{6M^{2}_{S_{1}}}+c_{\gamma}^{2}\frac{G_{S}(x_{2j})}{6M^{2}_{S_{2}}}\right]\,,$

where the loop function is

$G_{S}(x)=\frac{2-9x+18x^{2}-11x^{3}+6x^{3}\log(x)}{6(1-x)^{4}}\,.$ (38)

The latter function is the same as $G_{2}(x)$ in Ref. haba.yamada . See
details in appendix A. As anticipated, there is no chiral enhancement for
these coefficients.

For the type IB seesaw model, we analogously obtain

$\displaystyle\frac{16\pi^{2}}{e}C_{\beta\alpha}^{\text{ND-}L}$
$\displaystyle=\frac{v_{2}^{2}}{v^{2}}\lambda^{\dagger}_{\beta
1}\lambda_{1\alpha}\left[c_{\gamma}^{2}\frac{G_{S}(x_{1})}{6M^{2}_{S_{1}}}+s_{\gamma}^{2}\frac{G_{S}(x_{2})}{6M^{2}_{S_{2}}}\right]\,,$
(39) $\displaystyle\frac{16\pi^{2}}{e}C_{\beta\alpha}^{\text{ND-}R}$
$\displaystyle=f^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
2}f^{*}_{2\alpha}\left[s_{\gamma}^{2}\frac{G_{S}(x_{1j})}{6M^{2}_{S_{1}}}+c_{\gamma}^{2}\frac{G_{S}(x_{2j})}{6M^{2}_{S_{2}}}\right]\,.$

We could also consider the decay $\mu\rightarrow eee$. In addition to the
Wilson coefficients that we already have, we will need to consider box
diagrams. These diagrams have righthanded neutrinos and charged scalars
inside, whose contributions will be of the form

$f_{i\mu}f_{ie}f_{je}f_{je},\quad\lambda_{i\mu}\lambda_{ie}\lambda_{je}\lambda_{je},\quad\lambda_{i\mu}f_{ie}\lambda_{je}f_{je},\quad
f_{i\mu}\lambda_{ie}f_{je}\lambda_{je}$

In Sec. V we will see that we need suppressed $f_{ie}$ to evade CLFV
constraints and still be able to account for $(g-2)_{\mu}$. In this case, only
the contribution with $\lambda$ survives which is subdominant haba.yamada .

##  V Solving $(g-2)_{\mu}$ avoiding CLFV

By considering only the dominant chirally enhanced contribution to the dipole
term, we can estimate the interplay between the necessary contribution to
$a_{\mu}$ and the necessary suppression to avoid significant CLFV processes.
The contribution to $a_{\mu}$ in (20) necessary to explain the experimental
deviation (1) requires

$|C^{\sigma R}_{\mu\mu}|\sim|\operatorname{\mathrm{Re}}C^{\sigma
R}_{\mu\mu}|\sim 2\times 10^{-9}\,\mathrm{GeV}^{-1}\,.$ (40)

On the other hand, the current limit on $\mu\to e\gamma$ requires

$\sqrt{|C^{\sigma R}_{\mu e}|^{2}+|C^{\sigma R}_{e\mu}|^{2}}<4\times
10^{-14}\,\mathrm{GeV}^{-1}\,.$ (41)

So it is necessary that

$\frac{\sqrt{|C^{\sigma R}_{\mu e}|^{2}+|C^{\sigma R}_{e\mu}|^{2}}}{|C^{\sigma
R}_{\mu\mu}|}\lesssim 2\times 10^{-5}\,.$ (42)

The expected future limit will decrease this number by a factor $2.6$. A
similar analysis on the current limit on $\mu e$ conversion in gold nuclei
requires

$\sqrt{|C^{\sigma R}_{\mu e}|^{2}+|C^{\sigma
R}_{e\mu}|^{2}}<10^{-12}\,\mathrm{GeV}^{-1}\,,$ (43)

which implies

$\frac{\sqrt{|C^{\sigma R}_{\mu e}|^{2}+|C^{\sigma R}_{e\mu}|^{2}}}{|C^{\sigma
R}_{\mu\mu}|}\lesssim 5\times 10^{-4}\,.$ (44)

The expected future limit in aluminium nuclei will greatly reduce the limit in
(43) to $6\times 10^{-15}\,\mathrm{GeV}^{-1}$ and then the bound in (44)
becomes $3\times 10^{-6}$.

Now, let us see how our models can satisfy the hierarchy (42) between flavor
changing and flavor conserving couplings to the muon.

We start with the type IB seesaw model. We assume the dominance of the
chirally enhanced contribution and require $|s_{\gamma}|\gg 10^{-4}$ due to
(28). In this model, to have appreciable $a_{\mu}$, we need

$\lambda^{*}_{1\mu}f^{*}_{2\mu}\sim O(1)\,.$ (45)

On the other hand, according to (42), we need
$|\lambda^{*}_{1\mu}f^{*}_{2e}|,|\lambda^{*}_{1e}f^{*}_{2\mu}|$ to be less
than $10^{-5}$. Since $\lambda_{1\mu}$ needs to be order one, we can adopt
$f_{2e}=0$ for simplicity. In this case,

$\frac{\sqrt{|C^{\sigma R}_{\mu e}|^{2}+|C^{\sigma R}_{e\mu}|^{2}}}{|C^{\sigma
R}_{\mu\mu}|^{2}}\approx\frac{|\lambda_{1e}|}{|\lambda_{1\mu}|}\gtrsim 0.1\,,$
(46)

from the dependence of $\lambda_{1\alpha}$ on neutrino masses and mixing for
NO or IO, cf. eqs. (17) or (18), and (42) is never satisfied. Turning on the
coupling $f_{2e}$ only worsens the situation and the type IB seesaw model
augmented with a singly charged singlet cannot account for $(g-2)_{\mu}$
without violating current bounds on $\mu\to e\gamma$.

We now turn to the $\nu$-2HDM model. For appreciable $a_{\mu}$, assuming the
dominance of chirally enhanced contribution, we need large $\mu\mu$ couplings:

${\lambda^{(1)}_{\mu 1}}^{\dagger}f^{*}_{1\mu},\,{\lambda^{(1)}_{\mu
2}}^{\dagger}f^{*}_{2\mu}\sim O(1)\,.$ (47)

For equal masses $M_{1}=M_{2}$, to suppress CLFV processes, we need suppressed
$\mu e$ couplings:

$|{\lambda^{(1)}_{\mu 1}}^{\dagger}f^{*}_{1e}+{\lambda^{(1)}_{\mu
2}}^{\dagger}f^{*}_{2e}|\lesssim
10^{-5}\,,\quad|{\lambda^{(1)}_{e1}}^{\dagger}f^{*}_{1\mu}+{\lambda^{(1)}_{e2}}^{\dagger}f^{*}_{2\mu}|\lesssim
10^{-5}\,.$ (48)

Considering (47), the first combination vanishes if

$f_{1e}=f_{2e}=0\,,$ (49)

while the vanishing of the second requires the orthogonality between
$f_{i\mu}^{*}$ and $\lambda^{(1)}_{ie}$:

$(f_{1\mu},f_{2\mu})=\zeta(\lambda^{(1)}_{2e},-\lambda^{(1)}_{1e})\,.$ (50)

The contribution to $a_{\mu}$ then is proportional to

${\lambda^{(1)}_{\mu 1}}^{\dagger}f^{*}_{1\mu}+{\lambda^{(1)}_{\mu
2}}^{\dagger}f^{*}_{2\mu}=\zeta({\lambda^{(1)}_{2e}}^{*}{\lambda^{(1)}_{1\mu}}^{*}-{\lambda^{(1)}_{1e}}^{*}{\lambda^{(1)}_{2\mu}}^{*})\,,$
(51)

which is mostly fixed from neutrino parameters. The mass degeneracy
$M_{1}=M_{2}$ and the orthogonality condition (50) could in principle be
justified by flavor symmetries feruglio  which needs much more structure and
will not be treated here. This links the coupling $f$ with $\lambda$ so that
all the terms in (25) scale as $1/v_{1}^{2}$. The vanishing (49) also makes
the BSM contribution to electron EDM negligible, making it easily compatible
with the current precise measurement d.e . For the muon EDM, the combination
(51) contributes but it is equally negligible compared to the current limit
d.mu .

Considering now the complete contribution to the dipole coefficient $C^{\sigma
R}_{\beta\alpha}$ in (25), the turning off (or suppression) of the electron
coupling (49) eliminates the $f^{2}$ contribution (25c) for $\alpha=e$ or
$\beta=e$. Analogously, the orthogonality condition (50) eliminates the mixed
(left-right) contribution (25a), i.e., the chiral enhanced contribution, to
the flavor transition $\mu\to e$ if $N_{1},N_{2}$ have equal masses. Then,
only the $\lambda^{2}$ term (left-left) in (25c), which is not enhanced,
contributes to $\mu\to e\gamma$ while the chiral enhanced part still
contributes dominantly to $a_{\mu}$. With the increasing of the mass
difference $\Delta M_{N}=M_{2}-M_{1}$, the chiral enhanced contribution also
increases rapidly for $\mu\to e\gamma$. In principle, we could change the
orthogonality condition (50) to keep the chiral enhanced part (25a) vanishing,
at the cost of devising a mass dependent condition. This vanishing could even
be extended to the whole dipole coefficient (25). We will keep the simple
orthogonality condition (50), and use the mass difference $\Delta M_{N}$ as a
quantifier of the degree of tuning.

With the orthogonality condition, the couplings $f_{i\alpha}$ are completely
determined by free variables related to neutrinos masses
($\operatorname{\mathrm{Im}}(z)$, $M_{1}=M_{2}=M_{N}$), the vev $v_{1}$, and
the scaling parameter $\zeta$. Note that the real part of $z$ is not physical
for $M_{1}=M_{2}$. Thus, for fixed values of $M_{N}$ and $v_{1}$, we can
impose perturbativity bounds $|f_{i\alpha}|<4\pi$ as a function of
$\operatorname{\mathrm{Im}}z$ and $\zeta$. If the masses of the charged
scalars and their mixing is also fixed, we can obtain further bounds by
requiring compatibility with $(g-2)_{\mu}$ and $\mu\rightarrow e\gamma$. For
both NO and IO, we choose

$M_{S_{1}}=350\;\mathrm{GeV}\,,\quad M_{S_{2}}=450\;\mathrm{GeV}\,,\quad
s_{\gamma}=0.1\,.$ (52)

Figure 3 shows allowed regions in the plane
$\operatorname{\mathrm{Im}}(z){-}\zeta$ for $(g-2)_{\mu}$ (blue), $\mu\to
e\gamma$ (pink) and $\mu e$ conversion in Au (yellow). Perturbativity for
$f_{i\alpha}$ (dashed curves) is also shown. The left (right) figure is for NO
(IO) for which $M_{N}=1\,\mathrm{TeV}$ and $v_{1}=10^{-3}\;\mathrm{GeV}$
($v_{1}=2\times 10^{-3}\;\mathrm{GeV}$). The mixing angles of the PMNS matrix
are fixed in the best-fit of valle.fit  while we choose $\delta=218^{\circ}$.
Using a different value of $\delta$ leads to different shapes and regions for
this plot and others that follow but the overall possibility of explaining
$(g-2)_{\mu}$ and avoiding CLFV do not change significantly. The change is
larger for NO. For IO the variation is not significant.

Figure 3: Allowed regions in the plane
$\operatorname{\mathrm{Im}}(z)\times\zeta$ for NO (left) and IO (right).
Scalar masses and $s_{\gamma}$ are fixed as (52) while $M_{N}=1\,\mathrm{TeV}$
and $v_{1}=x\times 10^{-3}\;\mathrm{GeV}$ with $x=1$ ($x=2$) for NO (IO). The
region above the gray (green) dashed lines are excluded by perturbativity on
the couplings $f_{1\mu}$ ($f_{2\mu})$. The orange region is allowed by the
present constraint on $\mu N\to eN$, while the light red region is related to
$\mu\rightarrow e\gamma$. The dark (light) blue region show the 1(2)-$\sigma$
region allowed by $(g-2)_{\mu}$. The crosses denote the benchmark points in
(53). The mixing angles of the PMNS are fixed in best-fit of valle.fit
whereas $\delta=218^{\circ}$.

As illustrated in the plot, the strongest constraint on
$\operatorname{\mathrm{Im}}(z)$ comes from $\mu\rightarrow e\gamma$. There is
no visible constraint on $\zeta$ from $\mu\rightarrow e\gamma$ because the
chiral enhanced term was chosen to vanish. The constraint on $\zeta$ will
mainly come from $(g-2)_{\mu}$, whose dominant chiral enhanced contribution
depends linearly on this variable. In turn, perturbativity bounds on $f_{ij}$
impose upper bounds on $\zeta$. For each NO and IO, we see that the following
benchmark points account for $(g-2)_{\mu}$ still evading the current CLFV
constraints:

(BM-NO) $\displaystyle\quad M_{N}=1\,\mathrm{TeV}\,,\quad
v_{1}=10^{-3}\,\mathrm{GeV}\,,\quad z=0.2i\,,\quad\zeta=60\,;$ (53) (BM-IO)
$\displaystyle\quad M_{N}=1\,\mathrm{TeV}\,,\quad v_{1}=2\times
10^{-3}\,\mathrm{GeV}\,,\quad z=-0.1i\,,\quad\zeta=60\,.$

Only $v_{1}$ and $z$ differ in the two points. The rest of the model
parameters are fixed as (52).

For comparison and to assess the degree of tuning, we show in Fig. 4 the
allowed regions analogous to Fig. 3, with the same parameters, except for the
lifting of the mass degeneracy of RHNs to $M_{2}=1.0001\,M_{1}$. We can see
that we lose compatibiliy between $(g-2)_{\mu}$ at 1$\sigma$ and $\mu\to
e\gamma$ as the chiral enhanced term of the latter will not be completely
canceled. The region of compatibility for $(g-2)_{\mu}$ is practically
unchanged. We fix $\operatorname{\mathrm{Re}}(z)=0$ but we have checked that
variation of $\operatorname{\mathrm{Re}}(z)$ may lead at most to a difference
of a factor two.

Figure 4: The same as Fig. 3, but for non-degenerate neutrino masses:
$M_{2}=1.0001\,M_{1}$. Left: NO. Right: IO.

We illustrate in figures 5-7 the interplay of the distinct contributions to
$g-2$, $\mu\rightarrow e\gamma$, and $\mu e$ conversion in Au respectively,
using the benchmarks defined in (52) and (53). In all the plots, we show the
influence of the $\zeta$ parameter in the various contributions for each
observable. Starting with the contributions to $g-2$, Fig. 5, we show that the
chiral enhanced term given by (25a) is dominant, easily explaining the present
anomaly. We also show in the figure that allowing for non-degenerate neutrino
masses affects very little. Considering $\mu\rightarrow e\gamma$, for non-
degenerate masses the chiral contribution can easily surpass the bound on
$\mu\rightarrow e\gamma$. For the case of degenerate masses, the present bound
can be avoided for this benchmark, even though for the future limit this would
not be the case. The important point to notice is that all the contributions
can be suppressed by increasing $v_{1}$ since all of them are proportional to
$v_{1}^{-2}$ under the orthogonality condition. However, since only the chiral
contribution is proportional to $s_{2\gamma}$, we can still compensate by
increasing the value of $s_{\gamma}$. This reasoning was applied when choosing
the benchmark for IO, where $v_{1}$ is higher than the case for NO. At last,
we comment on the present bounds on $\mu e$ conversion in Au, Fig. 7. We show
that the non-dipole contribution given by (37) is negligible compared to the
dipole contribution given by (25). Moreover, as in the case of $\mu\rightarrow
e\gamma$, if the neutrino masses are non-degenerate we can easily surpass the
present bound.

Figure 5: The absolute value of the contributions for $(g{-}2)_{\mu}$ as
defined in eq. (25) as a function of $\zeta$, for the benchmarks defined in
(52) and (53) (left is NO and right is IO). The red line stands for the
dominant chiral enhanced term (25a). The blue line refers to the $\lambda^{2}$
term, (25b), while the green line corresponds to the $f^{2}$ term, (25c). Both
the green and blue contributions are negative, while the red one is positive.
The dark (light) gray band corresponds to the 1(2)-$\sigma$ region allowed at
present for $(g{-}2)_{\mu}$. The curves assumes $M_{1}=M_{2}$ while the bands
around them show the variation for $|M_{2}/M_{1}-1|\leq 0.4$.

Figure 6: The contributions for $\mu\to e\gamma$ as defined in eq. (25) as a
function of $\zeta$, for the benchmarks defined in (52) and (53) (left is NO
and right is IO). The red line stands for the chiral enhanced term (25a) while
the blue line refers to the $\lambda^{2}$ term, cf. (25b). The dashed
(continuous) gray line corresponds to the present (future) limits. Here,
$\Delta=M_{2}/M_{1}-1$. The band around the blue line shows how much this
contribution varies if $\Delta$ varies within $|\Delta|\leq 0.4$.

Figure 7: The contributions for $\operatorname{Br}(\mu\rm{Au}{\to}e\rm{Au})$
as defined in eqs. (25) and (37) as a function of $\zeta$, for the benchmark
defined in (52) and (53) (left is NO and right is IO). The red line stands for
the dipole contribution (25), which is dominant. Moreover, if the neutrino
masses are chosen non-degenerate, it can violate the present bound, given by
the dashed gray line. The blue line is related to the non-dipole contribution,
(37). Here, $\Delta=M_{2}/M_{1}-1$. The band around the blue line shows how
much this contribution varies if $\Delta$ varies within $|\Delta|\leq 0.4$.

To illustrate how the CLFV constraints are sensitive to the mass difference
$M_{2}-M_{1}$ of the RHNs, we can choose for NO the scaling
$v_{1}=10^{-3}\,\mathrm{GeV}\sqrt{M_{1}/\mathrm{TeV}}$ such that the Yukawas
$\lambda^{(1)}$ in (13) has the overall scale fixed as

$(\lambda^{(1)}{\lambda^{(1)}}^{\dagger})_{11}=\frac{M_{1}}{v_{1}^{2}}\Big{(}m_{2}|c_{z}|^{2}+m_{3}|s_{z}|^{2}\Big{)}\sim
0.05\times|s_{z}|^{2}\,,$ (54)

and the benchmark (53) is attained for $M_{1}=1\,\mathrm{TeV}$. For the
benchmark value we can see that $\lambda^{(1)}\sim O(0.1)$ and $f$ is larger
by $\zeta$ if respecting the orthogonality condition. For IO we choose twice
the value for $v_{1}$ so that the benchmark (53) is also attained for
$M_{1}=1\,\mathrm{TeV}$.

We start with countour curves, shown in Fig. 8, for the current and future
limit for $\mu\to e\gamma$ as a function of $\operatorname{\mathrm{Im}}(z)$
and $M_{1}$ for different values of $\Delta=M_{2}/M_{1}-1$. We clearly see
that the curves move to the right as the mass difference increases showing
that the constraints get stronger with the mass difference. The benchmark
point, marked with a cross, is already excluded for $\Delta=10^{-4}$ for both
NO and IO. In the future, this benchmark will be excluded even for degenerate
masses. A curve similar to the blue curve can be also found in Ref.
haba.yamada  where they consider the pure $\nu$-2HDM without the charged
singlet and hence equivalent to our case turning off $f_{j\alpha}$. The shape
of the curve is not exactly the same because their treatment of the Casas-
Ibarra parametrization of the minimal case is not appropriate. For certain
curves, there are islands near $M_{1}\sim 10^{2.3}\,\mathrm{GeV}$ indicating a
destructive interference between the chiral enhanced contribution and the
$\lambda^{2}$ contribution in (25).

Figure 8: Contours for $\mu\to e\gamma$ for fixed $\zeta=60$, and varying
$M_{2}/M_{1}-1=0,10^{-4},10^{-3},10^{-2}$ (respectively in blue, green, red
and purple). The vev $v_{1}=x\times
10^{-3}\,\mathrm{GeV}\sqrt{M_{1}/\mathrm{TeV}}$ scales with $\sqrt{M_{1}}$ so
that the global scale for the Yukawa $\lambda^{(1)}$ is fixed. The continuous
curves denote $\operatorname{Br}(\mu\to e\gamma)=4.2\times 10^{-13}$ (current)
and the dashed ones denote $\operatorname{Br}(\mu\to e\gamma)=6\times
10^{-14}$ (future). The mass is $M_{N}=M_{1}$ and $z$ is the parameter in (13)
for NO (similarly for IO) with $\operatorname{\mathrm{Re}}(z)=0$. Left: NO
with $x=1$. Right: IO with $x=2$. The cross denotes the benchmark points
defined in (52) and (53) for NO and IO.

In Fig. 9, we show similar contours for $\mu e$ conversion in nuclei. We can
see that the future constraints will be much stronger than the current ones.
Currently, a small mass difference is still allowed but in the future the
benchmark will be easily excluded.

Figure 9: Contours for $\operatorname{Br}(\mu N{\to}eN)$ for fixed $\zeta=60$,
and varying $M_{2}/M_{1}-1=0,10^{-4},10^{-3},10^{-2}$ (respectively in blue,
green, red and purple). The vev $v_{1}=x\times
10^{-3}\,\mathrm{GeV}\sqrt{M_{1}/\mathrm{TeV}}$ scales with $\sqrt{M_{1}}$ so
that the global scale for the Yukawa $\lambda^{(1)}$ is fixed. The continuous
curves denote $N={\rm Au}$ (current) and the dashed denote $N={\rm Ti}$
(future). The mass is $M_{N}=M_{1}$ and $z$ is the parameter in (13) for NO
(similarly for IO) with $\operatorname{\mathrm{Re}}(z)=0$. Left: NO with
$x=1$. Right: IO with $x=2$. The cross denotes the benchmark points defined in
(52) and (53) for NO and IO.

Finally we briefly comment on the dependence of our results on the masses of
the charged scalars. The expression of the chiral enhanced contribution (25a)
involves a cancellation between the contributions of the two charged scalars
and it vanishes for degenerate masses. So up to a certain point, increasing
the mass difference leads to an increase in the contribution to $g-2$. This
information can be seen in Fig. 10 where we show $1\sigma$ regions to satisfy
the $(g-2)_{\mu}$ constraint in the $M_{S_{1}}$ and $M_{N}=M_{1}=M_{2}$ plane
keeping the ratio fixed as $M_{S_{2}}/M_{S_{1}}=450/350\approx 1.29$ (blue)
and $M_{S_{2}}/M_{S_{1}}=2$ (gray). We can see that the larger mass ratio
allows a larger compatibility region with larger masses for the scalars and
the RHNs. The benchmark points defined in (52) and (53) correspond to the
origin of the plot in the corner of the blue region. Following this point to
the right inside the blue region would still allow the compatibility to
$(g-2)_{\mu}$ but with decreasing contribution to CLFV. We should remark that
the mass difference of the charged scalars cannot be arbitrarily large as this
would require an increasingly large $\mu_{\varphi}$. The benchmark we have
chosen, given by (52), is conservative in this sense, since it assumes
$\mu_{\varphi}\sim-30\;\mathrm{GeV}$. For higher values, for instance,
$-\mu_{\varphi}\sim\mathrm{TeV}$, it will still be possible to explain
$(g-2)_{\mu}$ with masses of the charged scalars and the RHNs of
$\mathcal{O}(\mathrm{TeV})$ given all other parameters fixed as before. For
even higher values of $-\mu_{\varphi}\sim 20\;\mathrm{TeV}$, we can avoid
future constraints on $\mu e$ conversion in nuclei and still satisfy
$(g-2)_{\mu}$.

Figure 10: $1\sigma$ regions for $(g-2)_{\mu}$ in the $M_{S_{1}}$ and
$M_{N}=M_{1}=M_{2}$ plane for fixed $\zeta=60$, $s_{\gamma}=0.1$ and two
values of $M_{S_{2}}/M_{S_{1}}=450/350\approx 1.29$ (blue) and
$M_{S_{2}}/M_{S_{1}}=2$ (gray). The vev $v_{1}=x\times
10^{-3}\,\mathrm{GeV}\sqrt{M_{1}/\mathrm{TeV}}$ scales with $\sqrt{M_{1}}$ so
that the global scale for the Yukawa $\lambda^{(1)}$ is fixed. Left: NO with
$x=1$ and $z$ as in (53). Right: IO with $x=2$ and $z$ as in (53).

Before ending this section, let us briefly discuss the collider phenomenology.
Given new particles with masses at the TeV scale, many BSM states can be
searched for directly. Since the type IB seesaw model cannot account for
$(g-2)_{\mu}$, we consider only the neutrinophilic version in eq. (4) in which
the neutrinophilic Higgs doublet $\Phi_{1}\sim H_{\nu}$ essentially contains
only BSM states and their interactions with quarks are suppressed by the large
$\tan\beta$, just as in the type I 2HDM. The interaction with quarks of the
charged scalars within the singlet $\varphi^{+}$ are analogously suppressed
and occurs only through mixing. Therefore, the dominant production mechanism
at LHC for the charged scalars present in the extended neutrinophilic model
will be through pair production $q\bar{q}\to\gamma^{*}(Z^{*})\to S^{-}S^{+}$
charged.higgs , while the production of the gauge singlets $N_{R}$ will be
highly suppressed deppisch  by the light-heavy neutrino mixing $|V_{\ell
N}|\sim v_{1}\lambda^{(1)}/M_{N}\sim 10^{-6}$ for $M_{N}\sim\mathrm{TeV}$ and
$v_{1}\sim\mathrm{MeV}$. Therefore, we will not consider the decay of $N_{R}$.

The decay mode of the charged scalars $S_{i}^{+}$ depends on their masses
relative to the masses of the RHNs. If $M_{S_{i}}<M_{N_{j}}$, assumed in our
benchmarks, and not so large $\tan\beta=v_{2}/v_{1}$, the decay of $S_{i}$
occurs dominantly through the Yukawa coupling to $\Phi_{2}$, which contains a
tiny admixture of $S_{i}^{+}$. The dominant width to $tb$ quarks is roughly

$\Gamma(S_{i}^{+}\to
tb)\sim\frac{y_{t}^{2}}{\tan^{2}\beta}\frac{M_{H^{\pm}}}{16\pi}\,,$ (55)

where we ignore the suppression by $c_{\gamma}^{2}$ or $s_{\gamma}^{2}$. For
large $\tan\beta$, the dominant decay is through Yukawa couplings
$\lambda^{(1)}$ and $f$ suppressed by the small mixing between the light and
heavy neutrinos with width

$\Gamma(S_{i}^{+}\to\ell_{\alpha}^{+}\nu_{j})\approx\frac{M_{S_{i}}}{16\pi}\left[|\Gamma^{H}_{j\alpha}a_{i}|^{2}+|\Gamma^{\varphi}_{j\alpha}b_{i}|^{2}\right]\,,$
(56)

where
$\Gamma^{H}=v_{1}V^{\mbox{\scriptsize$\mathsf{T}$}}{\lambda^{(1)}}^{\mbox{\scriptsize$\mathsf{T}$}}M_{R}^{-1}\lambda^{(1)}=\hat{M}_{\nu}V^{\dagger}/v_{1}$,
$\Gamma^{\varphi}=v_{1}V^{\dagger}{\lambda^{(1)}}^{\dagger}M_{R}^{-1}f^{*}$,
$(a_{1},a_{2})=(c_{\gamma},-s_{\gamma})$,
$(b_{1},b_{2})=(s_{\gamma},c_{\gamma})$. The term with $\Gamma^{H}$ is already
present in the $\nu$-2HDM without the charged singlet nu2HDM:pheno . As we
need $f_{1\mu},f_{2\mu}\sim O(1)$ to solve $(g-2)_{\mu}$, the contribution
from $\Gamma^{\varphi}$ is larger. From a rough estimate, the decay (56) start
to dominate for $\tan\beta\gtrsim 10^{6}$. In contrast, if
$M_{S_{i}}>M_{N_{j}}$, the charged scalars decay rapidly through the Yukawa
couplings $\lambda^{(1)}$ and $f$ into charged leptons and $N_{R}$. The latter
will be long-lived and decay to $W^{+}\ell_{\alpha}$ or $Z\nu_{j}$ with rates
suppressed by the light-heavy neutrino mixing.

To conclude this section, we have seen that the solution for $(g-2)_{\mu}$
within our model requires a sensitive balance of parameters to evade the
strong bounds coming from CLFV processes. In general, increasing $s_{\gamma}$,
increasing $f_{j\alpha}$ or decreasing $v_{1}$ will increase the chiral
enhanced contribution which dominates in $(g-2)_{\mu}$ but also in the CLFV
processes if not properly suppressed.

##  VI Summary

In order to connect the mechanism of neutrino mass generation with the
$(g-2)_{\mu}$ anomaly, we proposed to add a single charged singlet to the
neutrinophilic 2HDM and a variant which implements a low scale seesaw by
attributing part of the smallness of neutrino masses to a small vev. We
studied two models: (a) the $\nu$-2HDM version and (b) the type IB seesaw
version, both for the minimal case of two righthanded neutrinos. A chiral
enhanced contribution to $(g-2)_{\mu}$ is generated through exchange of
charged scalars and righthanded neutrinos, the latter also participating in
the neutrino mass generation. As family lepton number breaking is also brought
to low scale, the chiral enhanced contribution generically leads to large
rates for CLFV processes. We find that the type IB seesaw implementation does
not have enough freedom to circumvent the constraints from CLFV while solving
the $(g-2)_{\mu}$ anomaly. The $\nu$-2HDM version, on the other hand, has
enough freedom to allow some special cases where $(g-2)_{\mu}$ anomaly can be
solved and yet avoiding the stringent CLFV processes. Even in these special
cases, the region of compatibility between $(g-2)_{\mu}$ and current CLFV is
very restricted. One region for some choices of parameters is given in Fig. 3.
It is clear the our solutions are not restricted to the minimal case of two
RHNs as this limit can be mimicked in the presence of three righthanded
neutrinos and more regions may open up. In the future, experiments of $\mu e$
conversion in different nuclei are expected to drastically improve the limits
and this kind of solution to the $(g-2)_{\mu}$ anomaly will be put to test.

###### Acknowledgements.

A.C. acknowledges support from National Council for Scientific and
Technological Development – CNPq through projects 166523/2020-8 and
201013/2022-3. G.D.C. acknowledges financial support by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code
001. C.C.N. acknowledges partial support by Brazilian Fapesp, grant
2014/19164-6, and CNPq, grant 312866/2022-4.

##  Appendix A Wilson coefficients

Here we describe briefly how we obtain the Wilson coefficients of the
effective photonic operators (19). The calculations are based on appendices B
and C.

###  A.1 $\nu$-2HDM model

For the dipole part, the contribution (25b) comes from the $LL$ part of the
penguin diagram leading to (64a) with the following couplings multiplying
(67):

$\displaystyle N_{k}-S_{1}^{+}:$ $\displaystyle\quad-e{\lambda^{(1)}_{\beta
k}}^{\dagger}\frac{v_{2}c_{\gamma}}{v}\lambda^{(1)}_{k\alpha}\frac{v_{2}c_{\gamma}}{v}\,,$
(57) $\displaystyle N_{k}-S_{2}^{+}:$
$\displaystyle\quad-e{\lambda^{(1)}_{\beta
k}}^{\dagger}\frac{v_{2}(-s_{\gamma})}{v}\lambda^{(1)}_{k\alpha}\frac{v_{2}(-s_{\gamma})}{v}\,.$

The masses should be attributed accordingly. The contribution (25c) comes from
the $RR$ part of the penguin diagram leading to (64a), $L$ – $R$ exchanged,
with the following couplings multiplying (68):

$\displaystyle N_{k}-S_{1}^{+}:$ $\displaystyle\quad-
ef^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
k}s_{\gamma}f^{*}_{k\alpha}s_{\gamma}\,,$ (58) $\displaystyle
N_{k}-S_{2}^{+}:$ $\displaystyle\quad-
ef^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
k}c_{\gamma}f^{*}_{k\alpha}c_{\gamma}\,.$

The chirally enhanced contribution (25a) comes from $LR$ part of the penguin
diagram leading to (64d), with the following couplings multiplying (70):

$\displaystyle N_{k}-S_{1}^{+}:$ $\displaystyle\quad-e{\lambda^{(1)}_{\beta
k}}^{\dagger}\frac{v_{2}c_{\gamma}}{v}f^{*}_{k\alpha}s_{\gamma}\,,$ (59)
$\displaystyle N_{k}-S_{2}^{+}:$ $\displaystyle\quad-e{\lambda^{(1)}_{\beta
k}}^{\dagger}\frac{v_{2}(-s_{\gamma})}{v}f^{*}_{k\alpha}c_{\gamma}\,.$

For the non-dipole part, there are only chirality preserving contributions.
The integral $I^{\mu}_{LL}$ (64a) and the coefficient (69) should be
multiplied by (57) while the analogous $I^{\mu}_{RR}$ should be multiplied by
(58), resulting in the coefficients (37).

###  A.2 Seesaw type IB

The contribution (30b) comes from the $LL$ part (64a) with the following
couplings multiplying (67):

$\displaystyle N-S_{1}^{+}:$ $\displaystyle\quad-e\lambda^{\dagger}_{\beta
1}\frac{v_{2}c_{\gamma}}{v}\lambda_{1\alpha}\frac{v_{2}c_{\gamma}}{v}\,,$ (60)
$\displaystyle N-S_{2}^{+}:$ $\displaystyle\quad-e\lambda^{\dagger}_{\beta
1}\frac{v_{2}(-s_{\gamma})}{v}\lambda_{1\alpha}\frac{v_{2}(-s_{\gamma})}{v}\,.$

The contribution (30c) comes from the $RR$ part of (64a), $L$ – $R$ exchanged,
with the following couplings multiplying (68):

$\displaystyle N-S_{1}^{+}:$ $\displaystyle\quad-
ef^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
2}s_{\gamma}f^{*}_{2\alpha}s_{\gamma}\,,$ (61) $\displaystyle N-S_{2}^{+}:$
$\displaystyle\quad-ef^{\mbox{\scriptsize$\mathsf{T}$}}_{\beta
2}c_{\gamma}f^{*}_{2\alpha}c_{\gamma}\,.$

The chirally enhanced contribution (30a) comes from $LR$ part of the penguin
diagram leading to (64d), with the following couplings multiplying (70):

$\displaystyle N-S_{1}^{+}:$ $\displaystyle\quad-e\lambda^{\dagger}_{\beta
1}\frac{v_{2}}{v}c_{\gamma}f^{*}_{2\alpha}s_{\gamma}\,,$ (62) $\displaystyle
N-S_{2}^{+}:$ $\displaystyle\quad-e\lambda^{\dagger}_{\beta
1}\frac{v_{2}}{v}(-s_{\gamma})f^{*}_{2\alpha}c_{\gamma}\,.$

Other combinations are forbidden by the pseudo-Dirac nature of $N$.

For the non-dipole part, there are only chirality preserving contributions.
The integral $I^{\mu}_{LL}$ (64a) and the coefficient (69) should be
multiplied by (57) while the analogous $I^{\mu}_{RR}$ should be multiplied by
(58), resulting in the coefficients (39).

##  Appendix B Loop integrals

Our calculations for loop integrals and operators are similar to Refs. bilenky
; raidal.santamaria . The following loop integrals come from self-energy
diagrams and penguin diagrams in Fig. 2:

$\displaystyle iI_{LL}$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{R\not{k}}{(k^{2}-M_{N}^{2})[(k-p)^{2}-M_{\varphi}^{2}]}$
(63a)
$\displaystyle=\frac{iR}{(4\pi)^{2}}\,\frac{\not{p}}{2}\left\\{\frac{1}{\epsilon}+\log\frac{\bar{\mu}^{2}}{M_{\varphi}^{2}}+h_{S}(x)+\frac{p^{2}}{M_{\varphi}^{2}}8\tilde{f}_{S}(x)\right\\}\,,$
(63b) $\displaystyle iI_{LR}$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{RM_{N}}{(k^{2}-M_{N}^{2})[(k-p)^{2}-M_{\varphi}^{2}]}$
(63c)
$\displaystyle=\frac{iR}{(4\pi)^{2}}M_{N}\left\\{\frac{1}{\epsilon}+\log\frac{\bar{\mu}^{2}}{M_{\varphi}^{2}}+\frac{1-x+x\log
x}{1-x}+\frac{p^{2}}{M_{\varphi}^{2}}2f_{S}(x)\right\\}\,,$ (63d)
$\displaystyle iI^{\mu}_{LL}$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{R\not{k}(p_{1}+p_{2}-2k)^{\mu}}{(k^{2}-M_{N}^{2})[(k-p_{1})^{2}-M_{\varphi}^{2}][(k-p_{2})^{2}-M_{\varphi}^{2}]}$
(64a)
$\displaystyle=-\frac{iR}{(4\pi)^{2}}\bigg{\\{}\mbox{\large$\tfrac{1}{2}$}\gamma^{\mu}\left[\frac{1}{\epsilon}+\log\frac{\bar{\mu}^{2}}{M_{\varphi}^{2}}+h_{S}(x)\right]+\frac{1}{M_{\varphi}^{2}}(q^{2}\gamma^{\mu}-\not{q}q^{\mu})\mbox{\large$\tfrac{1}{6}$}G_{S}(x)$
(64b) $\displaystyle\qquad\qquad+\
\frac{1}{M_{\varphi}^{2}}[(\not{p_{1}}+\not{p_{2}})(p_{1}+p_{2})^{\mu}+(p_{1}^{2}+p_{2}^{2})\gamma^{\mu}]2\tilde{f}_{S}(x)\bigg{\\}}\,,$
(64c) $\displaystyle iI^{\mu}_{LR}$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{RM_{N}(p_{1}+p_{2}-2k)^{\mu}}{(k^{2}-M_{N}^{2})[(k-p_{1})^{2}-M_{\varphi}^{2}][(k-p_{2})^{2}-M_{\varphi}^{2}]}$
(64d)
$\displaystyle=-\frac{iR}{(4\pi)^{2}}\frac{M_{N}}{M_{\varphi}^{2}}(p_{1}+p_{2})^{\mu}2f_{S}(x)\,,$
(64e)

where $x=M_{N}^{2}/M_{\varphi}^{2}$ and $q=p_{2}-p_{1}$. These expression
should be supplied with couplings and enclosed by spinors $\bar{u}(p_{2})$ and
$u(p_{1})$ to give the amplitudes. Simple chirality exchange $L\leftrightarrow
R$ leads to identical expressions with the projector exchanged. We use
dimensional regularization with $d=4-2\epsilon$ and retain only terms up to
$M_{\varphi}^{-2}$ or $M_{N}^{-2}$. The additional loop function that appear
is

$h_{S}(x)=\frac{1-4x+3x^{2}-2x^{2}\log x}{2(1-x)^{2}}\,.$ (65)

The result for $I_{LL}$ and $I^{\mu}_{LL}$ match bilenky  for $x\to 0$ as
$8\tilde{f}_{S}(0)=1/3$ and $G_{S}(0)=1/3$.

For the last integral in (64d), we can use the Gordon-type identity where we
can replace

$(p_{1}+p_{2})^{\mu}\to-i\sigma^{\mu\nu}q_{\nu}+\not{p_{2}}\gamma^{\mu}+\gamma^{\mu}\not{p_{1}}\,.$
(66)

##  Appendix C Operators

The operators relevant to CLFV are the photon interactions written in (19).
The operator in the first line is the dipole contribution whereas the ones in
the second line are the non-dipole (ND) part. These Wilson coefficients at
1-loop can be obtained by matching the full theory with the effective theory
through appropriate 1-loop amplitudes. The relevant amplitudes lead to the
expressions in (63) and (64) which comes from the self-energy and dipole
diagrams in Fig. 2. To obtain the full expressions, one only needs to add the
coupling constants and adapt the masses.

Let us focus first on the chirality preserving contributions. The coefficient
of $\not{p}$ in (63a) and of the $\gamma^{\mu}$ part in (64a) give rise to the
operator $\bar{\psi}_{L}\not{D}\psi_{L}$, $\psi$ being the collection of
lepton fields, and it should be removed by wave function renormalization. The
coefficient of $\not{p}p^{2}$ in (63a) and the last term in the square
brackets in (64a) lead to the operator
$\bar{\psi}_{L}(\not{D}D^{2}+D^{2}\not{D})\psi_{L}$ bilenky . This operator
can be replaced by $\bar{\psi}_{L}\not{D}^{3}\psi_{L}$ with additional dipole
contributions. The former operator does not lead to relevant physical
phenomena. The result is that the last term in $I^{\mu}_{LL}$ (64a) generates

$\delta C^{\sigma
R}_{\beta\alpha}=\frac{(-1)}{(4\pi)^{2}}\frac{1}{2M^{2}_{\varphi}}2\tilde{f}_{S}(x)m_{\alpha}\,,$
(67)

while an analogous $I^{\mu}_{RR}$ leads to

$\delta C^{\sigma
R}_{\beta\alpha}=m_{\beta}\frac{(-1)}{(4\pi)^{2}}\frac{1}{2M^{2}_{\varphi}}2\tilde{f}_{S}(x)\,.$
(68)

These contributions are chirality flipping after the use of the equations of
motion. Finally, the contribution proportional to $G_{S}(x)$ in (64a) leads to
the non-dipole term with

$C_{\beta\alpha}^{\text{ND-}L}=\frac{(-1)}{(4\pi)^{2}}\frac{1}{6M^{2}_{\varphi}}G_{S}(x)\,.$
(69)

A similar term comes from $I^{\mu}_{RR}$.

Let us now turn to the chirality flipping contributions. The coefficient of
$M_{N}$ in (63c) should be removed by lepton mass renormalization. The
coefficient of $p^{2}M_{N}$ in (63c) and the coefficient of
$\not{p_{2}}\gamma^{\mu}+\gamma^{\mu}\not{p_{1}}$ in (64d), after replacement
(66), leads to the operator $\bar{\psi}_{L}\not{D}^{2}\psi_{R}$ which is not
phenomenologically relevant. The dipole contribution (66) in (64d) leads to

$\delta C^{\sigma
R}_{\alpha\beta}=\frac{1}{(4\pi)^{2}}\frac{M_{N}}{M^{2}_{\varphi}}2f_{S}(x)\,.$
(70)

## References

  * (1) T. Ohlsson, Special Issue on “Neutrino Oscillations: Celebrating the Nobel Prize in Physics 2015” in Nucl. Phys. B 908 (2016), 1.
  * (2) S. Weinberg, Phys. Rev. Lett. 43 (1979), 1566-1570.
  * (3) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126, no.14, 141801 (2021) [arXiv:2104.03281 [hep-ex]]; D. P. Aguillard et al. [Muon g-2], “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.20 ppm,” https://muon-g-2.fnal.gov/result2023.pdf
  * (4) T. Aoyama, N. Asmussen, M. Benayoun, J. Bijnens, T. Blum, M. Bruno, I. Caprini, C. M. Carloni Calame, M. Cè and G. Colangelo, et al. Phys. Rept. 887, 1-166 (2020) [arXiv:2006.04822 [hep-ph]].
  * (5) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109 (2012), 111808 [arXiv:1205.5370 [hep-ph]], T. Aoyama, T. Kinoshita and M. Nio, Atoms 7 (2019) no.1, 28, A. Czarnecki, W. J. Marciano and A. Vainshtein, Phys. Rev. D 67 (2003), 073006 [erratum: Phys. Rev. D 73 (2006), 119901] [arXiv:hep-ph/0212229 [hep-ph]], C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, Phys. Rev. D 88 (2013), 053005 [arXiv:1306.5546 [hep-ph]], M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 77 (2017) no.12, 827 [arXiv:1706.09436 [hep-ph]], A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 97 (2018) no.11, 114025 [arXiv:1802.02995 [hep-ph]], G. Colangelo, M. Hoferichter and P. Stoffer, JHEP 02 (2019), 006 [arXiv:1810.00007 [hep-ph]], M. Hoferichter, B. L. Hoid and B. Kubis, JHEP 08 (2019), 137 [arXiv:1907.01556 [hep-ph]], M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 80 (2020) no.3, 241 [erratum: Eur. Phys. J. C 80 (2020) no.5, 410] [arXiv:1908.00921 [hep-ph]], A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 101 (2020) no.1, 014029 [arXiv:1911.00367 [hep-ph]], A. Kurz, T. Liu, P. Marquard and M. Steinhauser, Phys. Lett. B 734 (2014), 144-147 [arXiv:1403.6400 [hep-ph]], K. Melnikov and A. Vainshtein, Phys. Rev. D 70 (2004), 113006 [arXiv:hep-ph/0312226 [hep-ph]], P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95 (2017) no.5, 054026 [arXiv:1701.05829 [hep-ph]], G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, JHEP 04 (2017), 161 [arXiv:1702.07347 [hep-ph]], M. Hoferichter, B. L. Hoid, B. Kubis, S. Leupold and S. P. Schneider, JHEP 10 (2018), 141 [arXiv:1808.04823 [hep-ph]], A. Gérardin, H. B. Meyer and A. Nyffeler, Phys. Rev. D 100 (2019) no.3, 034520 [arXiv:1903.09471 [hep-lat]], J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, Phys. Lett. B 798 (2019), 134994 [arXiv:1908.03331 [hep-ph]], G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, JHEP 03 (2020), 101 [arXiv:1910.13432 [hep-ph]], T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung and C. Lehner, Phys. Rev. Lett. 124 (2020) no.13, 132002 [arXiv:1911.08123 [hep-lat]], G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, Phys. Lett. B 735 (2014), 90-91 [arXiv:1403.7512 [hep-ph]].
  * (6) S. Borsanyi, Z. Fodor, J. N. Guenther, C. Hoelbling, S. D. Katz, L. Lellouch, T. Lippert, K. Miura, L. Parato and K. K. Szabo, et al. Nature 593, no.7857, 51-55 (2021) [arXiv:2002.12347 [hep-lat]].
  * (7) A. Bazavov, C. Davies, C. DeTar, A. X. El-Khadra, E. Gámiz, S. Gottlieb, W. I. Jay, H. Jeong, A. S. Kronfeld and S. Lahert, et al. [arXiv:2301.08274 [hep-lat]], T. Blum, P. A. Boyle, M. Bruno, D. Giusti, V. Gülpers, R. C. Hill, T. Izubuchi, Y. C. Jang, L. Jin and C. Jung, et al. [arXiv:2301.08696 [hep-lat]], C. Alexandrou, S. Bacchio, P. Dimopoulos, J. Finkenrath, R. Frezzotti, G. Gagliardi, M. Garofalo, K. Hadjiyiannakou, B. Kostrzewa and K. Jansen, et al. [arXiv:2206.15084 [hep-lat]], M. Cè, A. Gérardin, G. von Hippel, R. J. Hudspith, S. Kuberski, H. B. Meyer, K. Miura, D. Mohler, K. Ottnad and P. Srijit, et al. Phys. Rev. D 106, no.11, 114502 (2022) [arXiv:2206.06582 [hep-lat]], C. Aubin, T. Blum, M. Golterman and S. Peris, Phys. Rev. D 106, no.5, 054503 (2022) [arXiv:2204.12256 [hep-lat]], G. Wang et al. [chiQCD], [arXiv:2204.01280 [hep-lat]], C. Lehner and A. S. Meyer, Phys. Rev. D 101, 074515 (2020) [arXiv:2003.04177 [hep-lat]].
  * (8) G. Colangelo, A. X. El-Khadra, M. Hoferichter, A. Keshavarzi, C. Lehner, P. Stoffer and T. Teubner, Phys. Lett. B 833 (2022), 137313 doi:10.1016/j.physletb.2022.137313 [arXiv:2205.12963 [hep-ph]].
  * (9) K. Kannike, M. Raidal, D. M. Straub and A. Strumia, JHEP 02 (2012), 106 [erratum: JHEP 10 (2012), 136] [arXiv:1111.2551 [hep-ph]].
  * (10) A. Freitas, J. Lykken, S. Kell and S. Westhoff, JHEP 05 (2014), 145 [erratum: JHEP 09 (2014), 155] [arXiv:1402.7065 [hep-ph]].
  * (11) A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, Phys. Rev. D 98, no.11, 113002 (2018) [arXiv:1807.11484 [hep-ph]].
  * (12) L. Calibbi, R. Ziegler and J. Zupan, JHEP 07 (2018), 046 [arXiv:1804.00009 [hep-ph]].
  * (13) P. Athron, C. Balázs, D. H. J. Jacob, W. Kotlarski, D. Stöckinger and H. Stöckinger-Kim, JHEP 09 (2021), 080 [arXiv:2104.03691 [hep-ph]].
  * (14) M. Lindner, M. Platscher and F. S. Queiroz, Phys. Rept. 731 (2018), 1-82 [arXiv:1610.06587 [hep-ph]].
  * (15) R. Capdevilla, D. Curtin, Y. Kahn and G. Krnjaic, Phys. Rev. D 105 (2022) no.1, 015028 [arXiv:2101.10334 [hep-ph]].
  * (16) W. Yin and M. Yamaguchi, Phys. Rev. D 106 (2022) no.3, 033007 [arXiv:2012.03928 [hep-ph]].
  * (17) K. S. Babu, S. Jana, M. Lindner and V. P. K, JHEP 10 (2021), 240 [arXiv:2104.03291 [hep-ph]]; G. Mohlabeng, Phys. Rev. D 99 (2019) no.11, 115001 [arXiv:1902.05075 [hep-ph]]; D. W. P. Amaral, D. G. Cerdeno, A. Cheek and P. Foldenauer, Eur. Phys. J. C 81 (2021) no.10, 861 [arXiv:2104.03297 [hep-ph]].
  * (18) C. Arbeláez, R. Cepedello, R. M. Fonseca and M. Hirsch, Phys. Rev. D 102 (2020) no.7, 075005 [arXiv:2007.11007 [hep-ph]].
  * (19) I. Bigaran and R. R. Volkas, Phys. Rev. D 102 (2020) no.7, 075037 [arXiv:2002.12544 [hep-ph]]; I. Doršner, S. Fajfer and S. Saad, Phys. Rev. D 102 (2020) no.7, 075007 [arXiv:2006.11624 [hep-ph]]; D. Zhang, JHEP 07 (2021), 069 [arXiv:2105.08670 [hep-ph]].
  * (20) R. Dermisek and A. Raval, Phys. Rev. D 88 (2013), 013017 [arXiv:1305.3522 [hep-ph]], R. Dermisek, K. Hermanek, N. McGinnis and N. McGinnis, Phys. Rev. Lett. 126 (2021) no.19, 191801 [arXiv:2011.11812 [hep-ph]].
  * (21) A. Broggio, E. J. Chun, M. Passera, K. M. Patel and S. K. Vempati, JHEP 11 (2014), 058 [arXiv:1409.3199 [hep-ph]].
  * (22) V. Ilisie, JHEP 04 (2015), 077 [arXiv:1502.04199 [hep-ph]].
  * (23) T. Han, S. K. Kang and J. Sayre, JHEP 02 (2016), 097 [arXiv:1511.05162 [hep-ph]].
  * (24) A. Cherchiglia, P. Kneschke, D. Stöckinger and H. Stöckinger-Kim, JHEP 01 (2017), 007 [erratum: JHEP 10 (2021), 242] [arXiv:1607.06292 [hep-ph]].
  * (25) A. Cherchiglia, D. Stöckinger and H. Stöckinger-Kim, Phys. Rev. D 98 (2018), 035001 [arXiv:1711.11567 [hep-ph]].
  * (26) T. Abe, R. Sato and K. Yagyu, JHEP 07 (2017), 012 [arXiv:1705.01469 [hep-ph]].
  * (27) D. Stockinger, J. Phys. G 34 (2007), R45-R92 [arXiv:hep-ph/0609168 [hep-ph]].
  * (28) M. Endo and W. Yin, JHEP 08 (2019), 122 [arXiv:1906.08768 [hep-ph]].
  * (29) J. Hisano and K. Tobe, Phys. Lett. B 510 (2001), 197-204 [arXiv:hep-ph/0102315 [hep-ph]].
  * (30) P. Minkowski, Phys. Lett. B 67 (1977), 421-428; T. Yanagida, Conf. Proc. C7902131 (1979) 95; M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927 (1979), 315-321 [arXiv:1306.4669 [hep-th]]; R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44 (1980), 912.
  * (31) M. Magg and C. Wetterich, Phys. Lett. B 94 (1980), 61-64; J. Schechter and J. W. F. Valle, Phys. Rev. D 22 (1980), 2227; G. Lazarides, Q. Shafi and C. Wetterich, Nucl. Phys. B 181 (1981), 287-300; R. N. Mohapatra and G. Senjanovic, Phys. Rev. D 23 (1981), 165; E. Ma and U. Sarkar, Phys. Rev. Lett. 80 (1998), 5716-5719 [arXiv:hep-ph/9802445 [hep-ph]].
  * (32) R. Foot, H. Lew, X. G. He and G. C. Joshi, Z. Phys. C 44 (1989), 441; E. Ma, Phys. Rev. Lett. 81 (1998), 1171-1174 [arXiv:hep-ph/9805219 [hep-ph]]; E. Ma and D. P. Roy, Nucl. Phys. B 644 (2002), 290-302 [arXiv:hep-ph/0206150 [hep-ph]]; T. Hambye, Y. Lin, A. Notari, M. Papucci and A. Strumia, Nucl. Phys. B 695 (2004), 169-191 [arXiv:hep-ph/0312203 [hep-ph]].
  * (33) T. Fukuyama, H. Sugiyama and K. Tsumura, JHEP 03 (2010), 044 [arXiv:0909.4943 [hep-ph]].
  * (34) R. Coy and M. Frigerio, Phys. Rev. D 99 (2019) no.9, 095040 [arXiv:1812.03165 [hep-ph]].
  * (35) D. Wyler and L. Wolfenstein, Nucl. Phys. B 218 (1983), 205-214; R. N. Mohapatra and J. W. F. Valle, Phys. Rev. D 34 (1986), 1642; M. C. Gonzalez-Garcia and J. W. F. Valle, Phys. Lett. B 216 (1989), 360-366
  * (36) J. Kersten and A. Y. Smirnov, Phys. Rev. D 76 (2007), 073005 [arXiv:0705.3221 [hep-ph]].
  * (37) Y. Cai, J. Herrero-García, M. A. Schmidt, A. Vicente and R. R. Volkas, Front. in Phys. 5 (2017), 63 [arXiv:1706.08524 [hep-ph]].
  * (38) K. S. Babu and J. Julio, Nucl. Phys. B 841 (2010), 130-156 [arXiv:1006.1092 [hep-ph]]; T. Nomura, H. Okada and Y. Orikasa, Phys. Rev. D 94 (2016) no.5, 055012 [arXiv:1605.02601 [hep-ph]]; T. Nomura and H. Okada, Phys. Rev. D 94 (2016), 075021 [arXiv:1607.04952 [hep-ph]]; S. Lee, T. Nomura and H. Okada, Nucl. Phys. B 931 (2018), 179-191 [arXiv:1702.03733 [hep-ph]]; C. W. Chiang, H. Okada and E. Senaha, Phys. Rev. D 96 (2017) no.1, 015002 [arXiv:1703.09153 [hep-ph]]; S. Saad, Phys. Rev. D 102 (2020) no.1, 015019 [arXiv:2005.04352 [hep-ph]]; S. Jana, P. K. Vishnu, W. Rodejohann and S. Saad, Phys. Rev. D 102 (2020) no.7, 075003 [arXiv:2008.02377 [hep-ph]]; A. Alvarez, A. Banik, R. Cepedello, B. Herrmann, W. Porod, M. Sarazin and M. Schnelke, JHEP 06 (2023), 163 [arXiv:2301.08485 [hep-ph]].
  * (39) R. K. Barman, R. Dcruz and A. Thapa, JHEP 03 (2022), 183 [arXiv:2112.04523 [hep-ph]].
  * (40) D. A. Dicus, H. J. He and J. N. Ng, Phys. Rev. Lett. 87 (2001), 111803 [arXiv:hep-ph/0103126 [hep-ph]].
  * (41) E. Ma, Phys. Rev. Lett. 86 (2001), 2502-2504 [arXiv:hep-ph/0011121 [hep-ph]].
  * (42) E. Ma and M. Raidal, Phys. Rev. Lett. 87 (2001), 011802 [erratum: Phys. Rev. Lett. 87 (2001), 159901] [arXiv:hep-ph/0102255 [hep-ph]].
  * (43) S. Gabriel and S. Nandi, Phys. Lett. B 655 (2007), 141-147 [arXiv:hep-ph/0610253 [hep-ph]]. F. Wang, W. Wang and J. M. Yang, EPL 76 (2006), 388-394 [arXiv:hep-ph/0601018 [hep-ph]];
  * (44) S. M. Davidson and H. E. Logan, Phys. Rev. D 80 (2009), 095008 [arXiv:0906.3335 [hep-ph]]; S. Kanemura, T. Matsui and H. Sugiyama, Phys. Lett. B 727 (2013), 151-156 [arXiv:1305.4521 [hep-ph]]; P. A. N. Machado, Y. F. Perez, O. Sumensari, Z. Tabrizi and R. Z. Funchal, JHEP 12 (2015), 160 [arXiv:1507.07550 [hep-ph]].
  * (45) N. Haba and K. Tsumura, JHEP 06 (2011), 068 [arXiv:1105.1409 [hep-ph]]; N. Haba and M. Hirotsu, Eur. Phys. J. C 69 (2010), 481-492 [arXiv:1005.1372 [hep-ph]].
  * (46) N. Haba, T. Omija and T. Yamada, PTEP 2020 (2020) no.7, 073B08 [arXiv:2004.01896 [hep-ph]].
  * (47) B. Fu and S. F. King, Phys. Rev. D 105 (2022) no.9, 095001 [arXiv:2107.01486 [hep-ph]; J. Hernandez-Garcia and S. F. King, JHEP 05 (2019), 169 [arXiv:1903.01474 [hep-ph]].
  * (48) P. Escribano, J. Terol-Calvo and A. Vicente, Phys. Rev. D 103 (2021) no.11, 115018 [arXiv:2104.03705 [hep-ph]]; C. H. Chen, C. W. Chiang and T. Nomura, Phys. Rev. D 104 (2021) no.5, 055011 [arXiv:2104.03275 [hep-ph]]; A. E. C. Hernández, S. F. King and H. Lee, Phys. Rev. D 103 (2021) no.11, 115024 [arXiv:2101.05819 [hep-ph]]; A. de Giorgi, L. Merlo and S. Pokorski, [arXiv:2211.03797 [hep-ph]]; T. Mondal and H. Okada, Nucl. Phys. B 976 (2022), 115716 [arXiv:2103.13149 [hep-ph]].
  * (49) A. Zee, Phys. Lett. B 93 (1980), 389 [erratum: Phys. Lett. B 95 (1980), 461].
  * (50) G. Abbiendi et al. [ALEPH, DELPHI, L3, OPAL and LEP], Eur. Phys. J. C 73 (2013), 2463 [arXiv:1301.6065 [hep-ex]].
  * (51) P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020)
  * (52) H. Augustin et al. [Mu3e], IEEE Trans. Nucl. Sci. 68 (2021), 1833-1840 [arXiv:2010.15648 [physics.ins-det]].
  * (53) Y. S. Amhis et al. [HFLAV], Eur. Phys. J. C 81, no.3, 226 (2021) [arXiv:1909.12524 [hep-ex]].
  * (54) E. Kou et al. [Belle-II], PTEP 2019 (2019) no.12, 123C01 [erratum: PTEP 2020 (2020) no.2, 029201] [arXiv:1808.10567 [hep-ex]].
  * (55) A. M. Baldini et al. [MEG II], Eur. Phys. J. C 78 (2018) no.5, 380.
  * (56) W. H. Bertl et al. [SINDRUM II], Eur. Phys. J. C 47 (2006), 337-346.
  * (57) M. Yucel et al. [Mu2e], PoS ICHEP2020 (2021), 439; R. J. Abrams et al. [Mu2e], [arXiv:1211.7019 [physics.ins-det]].
  * (58) A. Kurup [COMET], Nucl. Phys. B Proc. Suppl. 218 (2011), 38-43; R. Abramishvili et al. [COMET], PTEP 2020 (2020) no.3, 033C01 [arXiv:1812.09018 [physics.ins-det]]; M. Lee et al. [COMET and Mu2e], “Experimental Searches for Muon to Electron Conversion in a Nucleus: COMET, DeeMe, and Mu2e. A Contributed paper for Snowmass 21,” [arXiv:2203.07089 [hep-ex]].
  * (59) R. J. Barlow, Nucl. Phys. B Proc. Suppl. 218 (2011), 44-49
  * (60) Y. Kuno and Y. Okada, Rev. Mod. Phys. 73 (2001), 151-202
  * (61) S. Weinberg and G. Feinberg, Phys. Rev. Lett. 3, 111-114 (1959).
  * (62) V. Cirigliano, R. Kitano, Y. Okada and P. Tuzon, Phys. Rev. D 80 (2009), 013002 [arXiv:0904.0957 [hep-ph]].
  * (63) E. Bertuzzo, Y. F. Perez G., O. Sumensari and R. Zukanovich Funchal, JHEP 01 (2016), 018 [arXiv:1510.04284 [hep-ph]].
  * (64) E. Arganda, M. J. Herrero and A. M. Teixeira, JHEP 10 (2007), 104 [arXiv:0707.2955 [hep-ph]].
  * (65) R. Kitano, M. Koike and Y. Okada, Phys. Rev. D 66 (2002), 096002 [erratum: Phys. Rev. D 76 (2007), 059902] [arXiv:hep-ph/0203110 [hep-ph]].
  * (66) F. Feruglio and A. Romanino, Rev. Mod. Phys. 93 (2021) no.1, 015007 [arXiv:1912.06028 [hep-ph]].
  * (67) J. Baron et al. [ACME], Science 343 (2014), 269-272 [arXiv:1310.7534 [physics.atom-ph]]; V. Andreev et al. [ACME], Nature 562 (2018) no.7727, 355-360.
  * (68) G. W. Bennett et al. [Muon (g-2)], Phys. Rev. D 80 (2009), 052008 [arXiv:0811.1207 [hep-ex]].
  * (69) P. F. de Salas, D. V. Forero, C. A. Ternes, M. Tortola and J. W. F. Valle, Phys. Lett. B 782 (2018), 633-640 [arXiv:1708.01186 [hep-ph]].
  * (70) S. S. D. Willenbrock, Phys. Rev. D 35 (1987), 173; O. Brein and W. Hollik, Eur. Phys. J. C 13 (2000), 175-184 [arXiv:hep-ph/9908529 [hep-ph]].
  * (71) F. F. Deppisch, P. S. Bhupal Dev and A. Pilaftsis, New J. Phys. 17 (2015) no.7, 075019 [arXiv:1502.06541 [hep-ph]].
  * (72) M. S. Bilenky and A. Santamaria, Nucl. Phys. B 420 (1994), 47-93 [arXiv:hep-ph/9310302 [hep-ph]].
  * (73) M. Raidal and A. Santamaria, Phys. Lett. B 421 (1998), 250-258 [arXiv:hep-ph/9710389 [hep-ph]].