# Flavor Violating Di-Higgs Couplings

Fayez Abu-Ajamieh<EMAIL_ADDRESS>Centre for High Energy Physics,
Indian Institute of Science, Bangalore 560012, India  Marco Frasca
<EMAIL_ADDRESS>Rome, Italy  Sudhir K. Vempati<EMAIL_ADDRESS>Centre for High Energy Physics, Indian Institute of Science, Bangalore 560012,
India

###### Abstract

Di-Higgs couplings to fermions of the form $h^{2}\overline{f}f$ are absent in
the Standard Model, however, they are present in several physics Beyond
Standard Model (BSM) extensions, including those with vector-like fermions. In
Effective Field Theories (EFTs), such as the Standard Model Effective Field
Theory (SMEFT) and the Higgs Effective Field Theory (HEFT), these couplings
appear at dimension 6 and can in general, be flavour-violating (FV). In the
present work, we employ a bottom-up approach to investigate the FV in the
lepton and quarks sectors through the di-Higgs effective couplings. We assume
that all FV arises from this type of couplings and assume that the Yukawa
couplings $Y_{ij}$ are given by their SM values, i.e.
$Y_{ij}=\sqrt{2}m_{i}\delta_{ij}/v$. In the lepton sector, we set upper limits
on the Wilson coefficients $C_{ll^{\prime}}$ from $l\rightarrow 3l^{\prime}$
decays, $l\rightarrow l\gamma$ decays, muonium oscillations, the $(g-2)_{\mu}$
anomaly, LEP searches, muon conversion in nuclei, FV Higgs decays, and $Z$
decays. We also make projections on some of these coefficients from Belle II,
the Mu2e experiment and the LHC’s High Luminosity (HL) run. In the quark
sector, we set upper limits on the Wilson coefficients $C_{qq^{\prime}}$ from
meson oscillations and from $B$-physics searches. A key takeaway from this
study is that current and future experiments should set out to measure the
effective di-Higgs couplings $C_{ff^{\prime}}$, whether these couplings are FV
or flavour-conserving. We also present a matching between our formalism and
the SMEFT operators and show the bounds in both bases.

##  I Introduction

Flavour physics provides an essential probe for the Standard Model (SM) and
for new physics BSM. In the SM, flavor violation (FV) arises entirely through
the fermionic couplings to the Higgs bosons, i.e. through the Yukawa matrices.
These Yukawa matrices encode FV in the CKM matrix in the hadronic sector, and
in the UPMNS matrix in the leptonic sector. In physics BSM, any new source of
flavour violation is severely constrained. FV processes are well measured in
$\Delta F=1$ and $\Delta F=2$ transitions. Some of the most robust constraints
are obtained from $K^{0}-\bar{K}^{0}$ system in the quark sector, and from
$\mu\to e+\gamma$ in the leptonic sector. Other processes which are not flavor
violating ($\Delta F=0$) but still play an essential role in constraining new
physics, are the magnetic and electric dipole moments of leptons, nucleons,
atoms and molecules. To avoid strong constraints on new physics from flavor
physics, typically it is assumed to follow the paradigm of Minimal Flavor
ViolationDAmbrosio:2002vsn .

An interesting scenario would arise when non-minimal FV is induced through the
effective Higgs couplings to fermions. There are many new physics scenarios
where non-minimal FV can arise through the Higgs couplings, such as the multi-
Higgs models, the Randall-Sundrum models and so on. The case of FV couplings
with a single Higgs has been studied in Ref. Dery:2013rta ; Harnik:2012pb . FV
can be understood in terms of deviations of the SM Yukawa couplings from their
SM values in the generation space. A complete global analysis of flavor
observables was performed and the limits on the FV Yukawa couplings were
derived. This work is similar in theme with the analysis conducted in
Dery:2013rta ; Harnik:2012pb , and extends it to the case of FV through the
di-Higgs couplings to fermions.

Di-Higgs-fermion-fermion couplings are absent in the SM, however, they can be
generated in a way similar to the single Higgs couplings in many new physics
scenarios. A simple example of this are extensions of the SM with extra
vector-like fermions. In the limit of heavy vector-like fermions, integrating
them out would lead to operators with di-Higgs couplings to the SM
fermions111These are not the only set of operators after integrating the heavy
fermions. But we focus on these operators for the present discussion.. These
operators can be mapped to EFT frameworks, such as the HEFT and the SMEFT, at
the level of dimension six operator (see for example, deBlas:2017xtg ;
delAguila:2000aa ; Chen:2017hak ; Batell:2012ca  and the references therein).
The study of FV in EFTs has been performed in many works in the literature,
see for example Silvestrini:2018dos ; Descotes-Genon:2018foz ;
Aebischer:2018iyb ; Greljo:2023adz ; Greljo:2022cah ; Bruggisser:2021duo ;
Aoude:2020dwv ; Hurth:2019ula . To the best of our knowledge, non-minimal FV
di-Higgs couplings have never been studied previously in the literature, as in
most cases, these $h^{2}\overline{f}f$ operators are either avoided entirely,
or assumed to be proportional the Yukawa couplings by imposing (minimal)
flavor symmetries Aebischer:2020lsx ; Faroughy:2020ina .

Non-minimal di-Higgs couplings are interesting, as they have unique
signatures, and can be probed by future colliders, especially the muon
collider. A non-minimal di-Higgs coupling could even explain the discrepancy
of the muon $g-2$ anomaly Abu-Ajamieh:2022nmt . In studying these couplings in
the present work, we find it suitable to follow the framework proposed in
Chang:2019vez ; Abu-Ajamieh:2020yqi ; Abu-Ajamieh:2021egq ; Abu-
Ajamieh:2022ppp ; Abu-Ajamieh:2021vnh . We call this framework the Weak Scale
Deviations framework (WSD). This formalism is model-independent and bottom-up,
as it considers all possible deviations from the SM Lagrangian. The FV di-
Higgs couplings appear naturally in the expansion of the Higgs operator in
this formalism, along with deviations in the Yukawa couplings. While one could
choose to work within either the SMEFT or the HEFT, we find the WSD framework
to be more convenient and advantageous, as it has fewer assumptions compared
to either the SMEFT or the HEFT, and is more closely-linked to experiment as
we show later on. Nonetheless, we shall present the mapping of the WSD to the
SMEFT, and present the SMEFT cutoff scale that corresponds to the upper limits
on the FV di-Higgs Wilson coefficients for convenience.

Focusing on the di-Higgs couplings, we provide a complete analysis of the
flavor physics constraints for both the quark and the lepton sectors. Our
analysis follows similar lines as the analysis performed in Harnik:2012pb  for
FV Higgs Yukawa couplings. The results for the di-Higgs couplings are
presented in terms of the bounds on the Wilson coefficient of the
$h^{2}\overline{f}f$ operators and also to the corresponding UV scale in the
SMEFT. The bounds on the SMEFT operators are competitive and are similar to
those on new physics. For example, assuming the Wilson coefficients to be
$\mathcal{O}(1)$ in the SMEFT, the bounds on the UV scale $\Lambda$ range from
$\sim 1-10$ TeV in the leptonic sector, and can exceed $100$ TeV in the
$K^{0}-\bar{K^{0}}$ oscillations in the quark sector.

This paper is organised as follows: In Section II, we briefly review the WSD
formalism we utilize in this paper. In Section III, we present our complete
analysis on the FV through the di-Higgs couplings in the leptonic sector,
whereas in Section IV, we do the same analysis in the quark sector. In Section
V, and show the how this formalism can be mapped to the SMEFT framework, and
in particular derive the UV scale that corresponds to the upper limit on the
FV Wilson coefficient. Finally, we present our conclusions in Section VI. We
relegated much of the calculational details to the appendices A - D.

##  II Framework

We begin by introducing our FV framework, which is essentially based on the
phenomenological bottom-up WSD approach introduced in Chang:2019vez ; Abu-
Ajamieh:2020yqi ; Abu-Ajamieh:2021egq ; Abu-Ajamieh:2022ppp ; Abu-
Ajamieh:2021vnh , generalized to the case of FV couplings and Wilson
coefficients. In this framework, we avoid power expansion in writing down
higher-dimensional operators, as the case in the SMEFT. Instead, we
parameterize New Physics (NP) as deviations from the SM predictions without
making any references to any UV scale. Therefore, we write the most general FV
effective Lagrangian of the Yukawa interaction as follows

$\mathcal{L}_{\text{eff}}=-\frac{v}{\sqrt{2}}\Big{(}\overline{L}_{l}^{i}\tilde{\hat{H}}l_{R}^{j}+\text{h.c.}\Big{)}\Bigg{[}Y^{l}_{ij}\frac{X}{v}+C^{l}_{ij}\frac{X^{2}}{2!v^{2}}+\dots\Bigg{]}-\frac{v}{\sqrt{2}}\Big{(}\overline{Q}_{l}^{i}\tilde{\hat{H}}q_{R}^{j}+\text{h.c.}\Big{)}\Bigg{[}Y^{q}_{ij}\frac{X}{v}+C^{q}_{ij}\frac{X^{2}}{2!v^{2}}+\dots\Bigg{]},$
(1)

where $Y^{l}_{ij}$ and $Y^{q}_{ij}$ are the Yukawa coupling matrices for the
leptons and the quarks, respectively, whereas $C^{l}_{ij}$ and $C^{q}_{ij}$
are matrices containing FV Wilson coefficients that do not have SM
counterparts. Also notice that in the SM we have
$Y^{l,q}_{ij}=\delta_{ij}\sqrt{2}m_{i}/v$, and $C^{l,q}_{ij}=0$. The field $X$
is defined in terms of the Higgs doublet $H$ as

$X=\sqrt{2H^{\dagger}H}-v,$ (2)

whereas we define the projector $\tilde{\hat{H}}=\epsilon\hat{H}^{\ast}$, with

$\epsilon=\begin{pmatrix}0&&1\\\ -1&&0\end{pmatrix},\hskip
28.45274pt\hat{H}=\frac{H}{\sqrt{H^{\dagger}H}}=\begin{pmatrix}0\\\
1\end{pmatrix}+O(\vec{G}),$ (3)

where $\vec{G}$ are the Goldstone bosons. Notice that $X$ has the same quantum
numbers as the Higgs field, and in the unitary gauge we have $X\rightarrow h$.
Before we proceed, a few of remarks are in order.

  * •

Notice that in Eq. 1, we are dividing the field $X$ by appropriate powers of
$v$ in order to keep Wilson coefficients dimensionless, i.e., $v$ should not
be interpreted as an expansion scale as the case in the HEFT Grinstein:2007iv
, and the Wilson coefficients could in principle assume any value allowed by
unitarity and experiment,

  * •

We are assuming that $v$ is the minimum of Higgs potential including all
higher-order corrections. Therefore, $v=246$ GeV. In addition, the value Higgs
mass remains equal to the measured one, i.e. $125$ GeV,

  * •

Although Eq. (1) appear to be similar to the HEFT, we should keep in mind that
secretly we are using the Higgs doublet in our expansion, and one can easily
demonstrate that the effective Lagrangian in Eq. (1) can be mapped to either
the SMEFT or the HEFT, depending on the chosen expansion, i.e. Eq. (1) can be
mapped to SMEFT when $X\rightarrow H$, and can be mapped to HEFT when
$X\rightarrow h$, as the case when the unitary gauge is chosen, AND when $v$
is interpreted as a true expansion scale. In either the SMEFT or the HEFT
frameworks, the deviations and Wilson coefficients in eqs. (1) can receive
corrections from a tower of higher-order operators, which might be different
depending on the order at which we truncate the expansion. We will present the
matching to the SMEFT in Section V below and show the corresponding scale of
new physics. The interested reader is instructed to refer to Abu-
Ajamieh:2020yqi ; Abu-Ajamieh:2021egq ; Abu-Ajamieh:2022ppp ; Abu-
Ajamieh:2021vnh  for more details on mapping the operators into the SMEFT and
the HEFT.

  * •

There are two advantages to this construction: First, there are fewer
assumptions in this framework compared to either the SMEFT or the HEFT.
Namely, we are only assuming that there are no light degrees of freedom below
the energy scale at which the EFT breaks down, and that the deviations and
Wilson coefficients are compatible with experimental measurements. The second
benefit lies in the fact that parameterizing NP this way is more transparent
phenomenologically, and more closely-linked to experiment, as these deviations
and Wilson coefficients are what is measured experimentally as opposed to any
expansion scale.

Figure 1: Induced FV Yukawa couplings $Y^{\text{ind}}_{ij}$ through non-zero
$C_{ij}$.

It is commonly assumed in the literature that $Y_{ij}$ are the main source of
FV, and studies that investigate limits on $Y_{ij}$ abound (see for instance
Dery:2013rta ; Harnik:2012pb ; Zhang:2021nzv ; Vicente:2019ykr ; Soreq:2016rae
; Buschmann:2016uzg . In this paper however, we are more interested in the
case where the effective couplings $C_{ij}$ are the main source of FV.
Therefore, we assume

$Y_{ij}\simeq Y^{\text{SM}}_{ij}=\frac{\sqrt{2}m_{i}}{v}\delta_{ij},$ (4)

for both the quarks and the leptons. We call FV through the couplings
$h^{2}\overline{f}f$ the next-to-minimum FV through di-Higgs effective
couplings. The reason why it is not possible to make
$Y_{ij}=Y^{\text{SM}}_{ij}$ exactly, is that it is not possible to
simultaneously diagonalize both $Y_{ij}$ and $C_{ij}$, as non-zero $C_{ij}$
will induce corrections to $Y_{ij}$ at 2-loops as shown in Fig. 1. Let’s call
this part of the Yukawas $Y^{\text{ind}}_{ij}$ to distinguish it from any
corrections arising from any other source. We can estimate the size of
$Y^{\text{ind}}_{ij}$ as follows

$Y^{\text{ind}}_{ij}\sim\Big{(}\frac{1}{16\pi^{2}}\Big{)}^{2}\times
Y_{ff}C_{ff_{i}}C_{ff_{j}},$ (5)

which for $C_{ff_{i}},C_{ff_{j}}\sim O(1)$ implies that
$Y^{\text{ind}}_{ij}\sim O(10^{-5})$ at best, i.e. the FV contributions from
$Y^{\text{ind}}_{ij}$ are always suppressed compared to those arising from
$C_{ij}$ and are thus negligible. We will not concern ourselves with these
corrections in the remainder of this paper.

In the unitary gauge, the FV part of Eq. 1 reads

$\mathcal{L}_{\text{FV}}\supset-\frac{C^{l}_{ij}}{2\sqrt{2}v}\Big{(}\overline{l}^{i}_{L}l^{j}_{R}+\text{h.c.}\Big{)}h^{2}-\frac{C^{q}_{ij}}{2\sqrt{2}v}\Big{(}\overline{q}^{i}_{L}q^{j}_{R}+\text{h.c.}\Big{)}h^{2}.$
(6)

In general, the matrices $C^{l,q}_{ij}$ could be complex and needn’t be
symmetric. However, in this paper, we will simplify by assuming that they are
both real and symmetric, i.e. $\text{Im}(C^{l,q}_{ij})=0$ and
$C^{l,q}_{ij}=C^{l,q}_{ji}$.

##  III The Lepton Sector

We focus first on FV in the lepton sector. Explicitly, the lepton part of Eq.
(6) reads

$\mathcal{L}^{l}_{\text{FV}}\supset-\frac{1}{2\sqrt{2}v}\Big{[}C_{ee}\overline{e}e+C_{\mu\mu}\overline{\mu}\mu+C_{\tau\tau}\overline{\tau}{\tau}+C_{\mu
e}(\overline{e}\mu+\overline{\mu}e)+C_{\tau\mu}(\overline{\mu}\tau+\overline{\tau}\mu)+C_{\tau
e}(\overline{e}\tau+\overline{\tau}e)\Big{]}h^{2}.$ (7)

Notice that $C_{ll}$ are not FV, however, they will enter into the calculation
and bounds along with the FV couplings $C_{ll^{\prime}}$. The bounds on are
summarized in Table 1 and shown in Figures 9 and 10. Below, we discuss these
bounds in more detail.

Channel | Couplings | Bounds ($\Lambda$ TeV) | Projections ($\Lambda$ TeV)
---|---|---|---
$\tau\rightarrow 3\mu$ | $|C_{\tau\mu}||C_{\mu\mu}|$ |  $<2.54\times 10^{-2}$ $(>1.07)$ |  $<3.92\times 10^{-3}$ $(>1.7)$
$\mu\rightarrow 3e$ | $|C_{\mu e}||C_{ee}|$ |  $<4.41\times 10^{-5}$ $(>5.23)$ |  $<4.41\times 10^{-7}$ $(>16.53)$
$\tau\rightarrow 3e$ | $|C_{\tau e}||C_{ee}|$ |  $<2.88\times 10^{-2}$ $(>1.03)$ |  $<3.92\times 10^{-3}$ $(>1.7)$
$\tau^{-}\rightarrow e^{+}\mu^{-}\mu^{-}$ | $|C_{\tau\mu}||C_{\mu e}|$ |  $<2.29\times 10^{-2}$ $(>1.1)$ |  $<2.83\times 10^{-3}$ $(>1.85)$
$\tau^{-}\rightarrow\mu^{+}e^{-}e^{-}$ | $|C_{\tau e}||C_{\mu e}|$ |  $<2.15\times 10^{-2}$ $(>1.11)$ |  $<2.66\times 10^{-3}$ $(>1.88)$
$\tau^{-}\rightarrow\mu^{+}\mu^{-}e^{-}$ |  $|C_{\tau\mu}||C_{\mu e}|$, $|C_{\tau e}||C_{\mu\mu}|$ |  $<2.88\times 10^{-2}$ $(>1.03)$ |  $<3.72\times 10^{-3}$ $(>1.73)$
$\tau^{-}\rightarrow\mu^{-}e^{+}e^{-}$ |  $|C_{\tau\mu}||C_{ee}|$, $|C_{\tau e}||C_{\mu e}|$ |  $<2.35\times 10^{-2}$ $(>1.09)$ |  $<2.99\times 10^{-3}$ $(>1.82)$
$\mu\rightarrow e\gamma$ ($\tau$ in loop) | $|C_{\tau\mu}||C_{\tau e}|$ |  $<7.83\times 10^{-5}$ $(>4.53)$ |  $<2.7\times 10^{-5}$ $(>5.91)$
$\mu\rightarrow e\gamma$ ($\mu$ in loop) | $|C_{\mu\mu}||C_{\mu e}|$ |  $<4.4\times 10^{-4}$ $(>2.94)$ |  $<1.52\times 10^{-4}$ $(>3.84)$
$\mu\rightarrow e\gamma$ ($e$ in loop) | $|C_{\mu e}||C_{ee}|$ |  $<8.28\times 10^{-4}$ $(>2.51)$ |  $<2.86\times 10^{-4}$ $(>3.28)$
$\tau\rightarrow\mu\gamma$ ($\tau$ in loop) | $|C_{\tau\tau}||C_{\tau\mu}|$ |  $<0.66$ $(>0.47)$ |  $<9.92\times 10^{-2}$ $(>0.76)$
$\tau\rightarrow\mu\gamma$ ($\mu$ in loop) | $|C_{\tau\mu}||C_{\mu\mu}|$ |  $<1.12$ $(>0.41)$ |  $<0.17$ $(>0.66)$
$\tau\rightarrow\mu\gamma$ ($e$ in loop) | $|C_{\tau e}||C_{\mu e}|$ |  $<0.64$ $(>0.48)$ |  $<9.66\times 10^{-2}$ $(>0.76)$
$\tau\rightarrow e\gamma$ ($\tau$ in loop) | $|C_{\tau\tau}||C_{\tau e}|$ |  $<0.57$ $(>0.49)$ |  $<0.22$ $(>0.62)$
$\tau\rightarrow e\gamma$ ($\mu$ in loop) | $|C_{\tau\mu}||C_{\mu e}|$ |  $<0.97$ $(>0.43)$ |  $<0.38$ $(>0.54)$
$\tau\rightarrow e\gamma$ ($e$ in loop) | $|C_{\tau e}||C_{ee}|$ |  $<0.55$ $(>0.49)$ |  $<0.22$ $(>0.62)$
$M-\overline{M}$ oscillations | $|C_{\mu e}|$ |  $<0.39$ $(>0.68)$ | -
$(g-2)_{\mu}$ | $|C_{\tau\mu}|$ |  $0.26\pm 0.03$ $(>0.84)$ | -
$(g-2)_{\mu}$ | $|C_{\mu\mu}|$ |  $0.79\pm 0.1$ $(>0.48)$ | -
$(g-2)_{\mu}$ | $|C_{\mu e}|$ |  $6.34\pm 0.8$ $(>0.17)$ | -
LEP | $|C_{\tau e}|$ |  $<9.52$ $(>0.14)$ | -
LEP | $|C_{\mu e}|$ |  $<9.0$ $(>0.14)$ | -
LEP | $|C_{ee}|$ |  $<13.25$ $(>0.12)$ | -
$\mu\rightarrow e$ conversion in nuclei | $|C_{\mu e}|$ |  $<0.34$ $(>0.73)$ |  $<4.56\times 10^{-3}$ $(>6.31)$
$h\rightarrow\tau\mu$ | $|C_{\tau\mu}|$ |  $<0.67$ $(>0.52)$ |  $<0.23$ $(>0.89)$
$h\rightarrow\tau e$ | $|C_{\tau e}|$ |  $<1.04$ $(>0.42)$ |  $<0.23$ $(>0.89)$
$h\rightarrow\mu e$ | $|C_{\mu e}|$ |  $<0.25$ $(>0.85)$ |  $<7.3\times 10^{-2}$ $(>1.58)$
$h\rightarrow ee$ | $|C_{ee}|$ |  $<0.58$ $(>0.56)$ | -
$Z\rightarrow\tau^{+}\tau^{-}$ |  $|C_{\tau\tau}|$, $|C_{\tau\mu}|$, $|C_{\tau e}|$ |  $<7.9$ $(>0.15)$ | -
$Z\rightarrow\mu^{+}\mu^{-}$ |  $|C_{\tau\mu}|$, $|C_{\mu\mu}|$, $|C_{\mu e}|$ |  $<7.04$ $(>0.16)$ | -
$Z\rightarrow e^{+}e^{-}$ |  $|C_{\tau e}|$, $|C_{\mu e}|$, $|C_{ee}|$ |  $<5.62$ $(>0.18)$ | -
$Z\rightarrow\tau^{\pm}\mu^{\mp}$ | $|C_{\tau\mu}|$ |  $<0.11$ $(>1.28)$ | -
$Z\rightarrow\tau^{\pm}e^{\mp}$ | $|C_{\tau e}|$ |  $<9.65\times 10^{-2}$ $(>1.37)$ | -
$Z\rightarrow\mu^{\pm}e^{\mp}$ | $|C_{\mu e}|$ |  $<1.59\times 10^{-3}$ $(>10.69)$ | -
Table 1:  $90\%$ CL bounds and projections on the leptonic next-to-minimal FV
di-Higgs couplings and the corresponding lower limit on the scale of NP
$\Lambda$ from matching to the SMEFT.

###  III.1 Bounds from $l\rightarrow l_{1}l_{2}l_{3}$ decays

The $l\rightarrow l_{1}l_{2}l_{3}$ decay through the di-Higgs couplings
proceeds at one loop as in Figure 2. Here, the $h^{2}ll^{\prime}$ vertices
should be viewed as effective interactions of some heavy degree(s) of freedom
that has been integrated out. In the limit $M_{h}\gg m_{l}$, the decay width
can be approximated as

Figure 2: The $l\rightarrow l_{1}l_{2}l_{3}$ decay through the di-Higgs
effective couplings $C_{ij}$.

$\Gamma(l\rightarrow
l_{1}l_{2}l_{3})\simeq\frac{m_{l}^{5}}{v^{4}}\Bigg{[}\frac{C_{ll_{1}}C_{l_{2}l_{3}}}{512\pi^{3}\sqrt{6\pi}}\log{\Big{(}\frac{M_{h}^{2}}{m_{l}^{2}}\Big{)}}\Bigg{]}^{2}.$
(8)

The detailed calculation is given in Appendix A. Before we proceed with
extracting the bounds, we should note that the 2-loop diagram (similar to the
bottom diagram in Figure 3, with the photon decaying to $l_{2}l_{3}$) is
suppressed relative to the 1-loop diagram and can be neglected.

The relevant processes are $\tau\rightarrow 3\mu$, $\mu\rightarrow 3e$,
$\tau\rightarrow\mu\mu e$, $\tau\rightarrow\mu ee$ and $\tau\rightarrow 3e$.
The latest bounds on the branching rations of these processes can be found in
ParticleDataGroup:2018ovx , and all of them are given $@$ $90\%$ C.L., which
we stick to throughout this paper. For the first process, the experimental
bound is $\text{Br}(\tau\rightarrow 3\mu)<2.1\times 10^{-8}$, which translates
into the bound $|C_{\tau\mu}||C_{\mu\mu}|<2.54\times 10^{-2}$. Notice that the
FV coupling $|C_{\tau\mu}|$ cannot be isolated from the non-FV one
$|C_{\mu\mu}|$. This is a common feature of these types of couplings. The
second experimental limit is given by $\text{Br}(\mu\rightarrow 3e)<1\times
10^{-12}$, which translates to the $|C_{\mu e}||C_{ee}|<4.41\times 10^{-5}$.
The limit on the third process is
$\text{Br}(\tau^{-}\rightarrow\mu^{-}\mu^{-}e^{+})<1.7\times 10^{-8}$,
translating into $|C_{\tau\mu}||C_{\mu e}|<2.29\times 10^{-2}$. The limit on
the fourth decay is $\text{Br}(\tau\rightarrow 3e)<2.7\times 10^{-8}$,
yielding the bound $|C_{\tau e}||C_{ee}|<2.88\times 10^{-2}$. the bounds on
the fifth process read
$\text{Br}(\tau^{-}\rightarrow\mu^{+}e^{-}e^{-})<1.5\times 10^{-8}$ and
translate to the limit $|C_{\tau e}||C_{\mu e}|<2.15\times 10^{-2}$.

The last 2 decays are more subtle as they involve two Feynman diagrams instead
of one. The decay width is obtained by summing two matrix element which have
different FV couplings. For the decay
$\tau^{-}\rightarrow\mu^{+}\mu^{-}e^{-}$, in the first diagram, we have
$l=\tau^{-}$, $l_{1}=e^{-}$, $l_{2}=\mu^{+}$, $l_{3}=\mu^{-}$, whereas in the
second we have $l=\tau^{-}$, $l_{1}=\mu^{-}$, $l_{2}=\mu^{+}$, $l_{3}=e^{-}$.
The experimental limit is
$\text{Br}(\tau^{-}\rightarrow\mu^{+}\mu^{-}e^{-})<2.7\times 10^{-8}$, which
translates into the bound $[2C_{\tau\mu}^{2}C_{\mu
e}^{2}+2C_{\mu\mu}^{2}C_{\tau e}^{2}-C_{\tau\mu}C_{\tau e}C_{\mu\mu}C_{\mu
e}]^{1/2}<4.07\times 10^{-2}$. Upper bounds can be obtained by setting
$C_{\tau\mu}=C_{\mu e}=0$ ($C_{\tau e}=C_{\mu\mu}=0$) in the first (second)
diagrams, which yields the bounds $|C_{\tau\mu}||C_{\mu e}|,|C_{\tau
e}||C_{\mu\mu}|<2.88\times 10^{-2}$. In the final process
$\tau^{-}\rightarrow\mu^{-}e^{+}e^{-}$, the two Feynman diagrams are given by
$l=\tau^{-}$, $l_{1}=e^{-}$, $l_{2}=e^{+}$, $l_{3}=\mu^{-}$ in the first
diagram, and $l=\tau^{-}$, $l_{1}=\mu^{-}$, $l_{2}=e^{+}$, $l_{3}=e^{-}$. The
experimental bound for this process is
$\text{Br}(\tau^{-}\rightarrow\mu^{-}e^{+}e^{-})<1.8\times 10^{-8}$, which
translates into the mixed bound $[2C_{\tau\mu}^{2}C_{ee}^{2}+2C_{\tau
e}^{2}C_{\mu e}^{2}-C_{\tau\mu}C_{\tau e}C_{\mu e}C_{ee}]^{1/2}<3.33\times
10^{-3}$, from which the upper bounds $|C_{\tau\mu}||C_{ee}|,|C_{\tau
e}||C_{\mu e}|<2.35\times 10^{-2}$ are obtained.

Better bounds can be obtained from future experiments. In particular, the
Belle II experiment Aushev:2010bq ; Belle-II:2022cgf  is expected to collect
$50\hskip 2.84526pt\text{ab}^{-1}$ over the next decade, and the bounds on the
branching rations of the above processes are projected to be $\sim
O(10^{-10})$ (see also Calibbi:2017uvl ; Banerjee:2022vdd 222The projections
provided in these two references are slightly different. For our projected
limits, we use the stronger of the two.). This leads to bounds that are 1-2
orders of magnitude stronger that what is currently available. For instance,
the projected bound from Belle II for $\text{Br}(\tau\rightarrow 3\mu)$ is
$5\times 10^{-10}$. This yields the projected bound
$|C_{\tau\mu}||C_{\mu\mu}|<3.92\times 10^{-3}$. The rest of the projections
are summarized in Table 1.

###  III.2 Bounds from $l_{i}\rightarrow l_{k}\gamma$

Stringent constraints can be obtained from the bounds on the FV decays
$\tau\rightarrow\mu\gamma$, $\tau\rightarrow e\gamma$ and $\mu\rightarrow
e\gamma$. The Feynman diagrams of these processes are shown in Figure 3. The
1-loop contributions are shown on the top row of the figure, where the photon
could be emitted from the initial or final state lepton. The two contributions
cancel one another and the contribution at one loop vanishes. Thus the leading
contribution arises at 2-loops 333Notice that there are two more 2-loop
diagrams where the photon is emitted from the initial and final states,
however, these two contribution cancel each other in exactly the same manner
as in the 1-loop case.

Figure 3: FV decays $l_{i}\rightarrow\l_{k}\gamma$ through one (top) and two
loops (bottom).

Calculating the 2-loop diagram is somewhat subtle and we show the details in
Appendix B. For each decay process, the structure of the matrix element and
the corresponding Wilson coefficients depend on the lepton inside the loop,
i.e., each decay will have 3 contributions corresponding to setting the
particle in the loop $j=\\{\tau,\mu,e\\}$. In order to set upper bounds on the
Wilson coefficients, we isolate each contribution individually. This will lead
to 9 different decay processes. For example, the decay width $\Gamma_{\tau\mu
e}$ refers to the decay $\tau\rightarrow e\gamma$ with $\mu$ running in the
loop.

Utilizing the results in Appendix B, assuming $m_{\tau}\gg m_{\mu}\gg m_{e}$,
and setting the renormalization scale $\mu=m_{j}$, the decay widths are given
by

$\displaystyle\Gamma_{\mu\tau e}$
$\displaystyle\simeq\frac{\alpha|C_{\tau\mu}|^{2}|C_{\tau
e}|^{2}}{16(4\pi)^{8}v^{4}}m_{\tau}^{2}m_{\mu}^{3}\Bigg{[}\log{\Big{(}\frac{M_{h}^{2}}{m_{\tau}^{2}}\Big{)}}-\frac{\pi^{2}}{3}\Bigg{]}^{2},$
(9) $\displaystyle\Gamma_{\mu\mu e}$
$\displaystyle\simeq\frac{\alpha|C_{\mu\mu}|^{2}|C_{\mu
e}|^{2}}{9(4\pi)^{8}v^{4}}m_{\mu}^{5}\Bigg{[}\log{\Big{(}\frac{M_{h}^{2}}{m_{\mu}^{2}}\Big{)}}-\frac{\pi^{2}}{4}\Bigg{]}^{2},$
(10) $\displaystyle\Gamma_{\mu ee}$ $\displaystyle\simeq\frac{\alpha|C_{\mu
e}|^{2}|C_{ee}|^{2}}{144(4\pi)^{8}v^{4}}m_{\mu}^{5}\log^{2}{\Big{(}\frac{M_{h}^{2}}{m_{e}^{2}}\Big{)}},$
(11) $\displaystyle\Gamma_{\tau\tau\mu}$
$\displaystyle\simeq\frac{\alpha|C_{\tau\tau}|^{2}|C_{\tau\mu}|^{2}}{9(4\pi)^{8}v^{4}}m_{\tau}^{5}\Bigg{[}\log{\Big{(}\frac{M_{h}^{2}}{m_{\tau}^{2}}\Big{)}}-\frac{\pi^{2}}{4}\Bigg{]}^{2},$
(12) $\displaystyle\Gamma_{\tau\mu\mu}$
$\displaystyle\simeq\frac{\alpha|C_{\tau\mu}|^{2}|C_{\mu\mu}|^{2}}{144(4\pi)^{8}v^{4}}m_{\tau}^{5}\log^{2}{\Big{(}\frac{M_{h}^{2}}{m_{\mu}^{2}}\Big{)}},$
(13) $\displaystyle\Gamma_{\tau e\mu}$
$\displaystyle\simeq\frac{\alpha|C_{\tau e}|^{2}|C_{\mu
e}|^{2}}{144(4\pi)^{8}v^{4}}m_{\tau}^{5}\log^{2}{\Big{(}\frac{M_{h}^{2}}{m_{e}^{2}}\Big{)}},$
(14) $\displaystyle\Gamma_{\tau\tau e}$
$\displaystyle\simeq\frac{\alpha|C_{\tau\tau}|^{2}|C_{\tau
e}|^{2}}{9(4\pi)^{8}v^{4}}m_{\tau}^{5}\Bigg{[}\log{\Big{(}\frac{M_{h}^{2}}{m_{\tau}^{2}}\Big{)}}-\frac{\pi^{2}}{4}\Bigg{]}^{2},$
(15) $\displaystyle\Gamma_{\tau\mu e}$
$\displaystyle\simeq\frac{\alpha|C_{\tau\mu}|^{2}|C_{\mu
e}|^{2}}{144(4\pi)^{8}v^{4}}m_{\tau}^{5}\log^{2}{\Big{(}\frac{M_{h}^{2}}{m_{\mu}^{2}}\Big{)}},$
(16) $\displaystyle\Gamma_{\tau ee}$ $\displaystyle\simeq\frac{\alpha|C_{\tau
e}|^{2}|C_{ee}|^{2}}{144(4\pi)^{8}v^{4}}m_{\tau}^{5}\log^{2}{\Big{(}\frac{M_{h}^{2}}{m_{e}^{2}}\Big{)}}.$
(17)

The experimental limits $\text{Br}(\mu\rightarrow e\gamma)<4.2\times 10^{-13}$
ParticleDataGroup:2012pjm  can be used in the decays (9), (10) and (11). The
decay $\Gamma_{\mu\tau e}$ yields bounds $|C_{\tau\mu}||C_{\tau e}|<7.83\times
10^{-5}$, whereas $\Gamma_{\mu\mu e}$ yields $|C_{\mu\mu}||C_{\mu
e}|<4.4\times 10^{-4}$, and $\Gamma_{\mu ee}$ translates to $|C_{\mu
e}||C_{ee}|<8.28\times 10^{-4}$. On the other hand, the limit
$\text{Br}(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$
ParticleDataGroup:2012pjm  can be used in the decays (12), (13) and (14), with
the decay $\Gamma_{\tau\tau\mu}$ leading to the bound
$|C_{\tau\tau}||C_{\tau\mu}|<0.66$ and the decay $\Gamma_{\tau\mu\mu}$ leading
to the bound $|C_{\tau\mu}||C_{\mu\mu}|<1.12$, whereas the decay $\Gamma_{\tau
e\mu}$ leads to the bound $|C_{\tau e}||C_{\mu e}|<0.64$. Finally, the
experimental limits $\text{Br}(\tau\rightarrow e\gamma)<3.3\times 10^{-8}$ is
used in last 3 decays in Eqs. (15), (16) and (17), with $\Gamma_{\tau\tau e}$
yielding the bound $|C_{\tau\tau}||C_{\tau e}|<0.57$, $\Gamma_{\tau\mu e}$
yielding the bound $|C_{\tau\mu}||C_{\mu e}|<0.97$ and finally $\Gamma_{\tau
ee}$ yielding the bound $|C_{\tau e}||C_{ee}|<0.55$.

Notice that bounds obtained here are roughly an order of magnitude weaker than
the bounds obtained from $l\rightarrow 3l^{\prime}$ decays. The reason for
this is that the former case proceeds through two loops, whereas the latter
proceeds through one loop.

As the case with the decays $l\rightarrow l_{1}l_{2}l_{3}$, the Belle II
experiment is projected to provide stronger bounds Aushev:2010bq ; Belle-
II:2022cgf ; Calibbi:2017uvl ; Calibbi:2017uvl , with projected branching
ratios that are about an order of magnitude stronger than the current limits.
For example, the projected Belle II constraints on the decay
$\tau\rightarrow\mu\gamma$ are $\text{Br}(\tau\rightarrow\mu\gamma)<1\times
10^{-9}$. This can be used in $\Gamma_{\tau\tau\mu}$, $\Gamma_{\tau\mu\mu}$
and $\Gamma_{\tau e\mu}$ to yield the projections
$|C_{\tau\tau}||C_{\tau\mu}|<9.92\times 10^{-2}$,
$|C_{\tau\mu}||C_{\mu\mu}|<0.17$ and $|C_{\tau e}||C_{\mu e}|<9.66\times
10^{-2}$, respectively. The projected limits are summarized in Table 1.

###  III.3 Constraints from muonium-antimuonium oscillations

$\mu^{+}$ and $e^{-}$ can form a bound state called muonium. This bound state
can oscillate to antimuonium $\mu^{-}e^{+}$ through the diagrams shown in
Figure 4, with $f_{i}=e^{-}$, $\overline{f}_{j}=\mu^{+}$, $f_{k}=\mu^{-}$ and
$\overline{f}_{l}=e^{+}$. The time-integrated $M-\overline{M}$ conversion
probability is constrained by the MACS experiment at PSI Willmann:1998gd

$P(M\rightarrow\overline{M})<8.3\times 10^{-11}/S_{B},$ (18)

where $S_{B}$ accounts for the splitting of muonium in the magnetic field of
the detectors, and is given by $S_{B}=0.35$ for $(S\pm P)\times(S\pm P)$
operators and $S_{B}=0.9$ for $P\times P$ operators. In this paper, we chose
to be conservative and set $S_{B}=0.35$. The loops in the s- and t-channels in
Figure 4 are given by Eq. 67, which can be integrated out in the non-
relativistic limit, yielding the following effective Lagrangian

$\mathcal{L}_{\text{eff}}=\frac{C_{\mu
e}^{2}}{32\pi^{2}v^{2}}\log{\Big{(}\frac{m_{\mu}^{2}}{M_{h}^{2}}\Big{)}}[\overline{\mu}e][\overline{e}\mu],$
(19)

Figure 4: The $f_{i}\overline{f}_{j}\rightarrow f_{k}\overline{f}_{l}$
scattering through the $h^{2}ff$ couplings. The left diagram is the s-channel,
whereas the right diagram is the t-channel.

where we have set the renormalization scale $\mu^{2}=m_{\mu}^{2}$. The
theoretical prediction for the conversion rate is governed by the matrix
element

$\mathcal{M}_{M}=\bra{\uparrow_{\mu}\downarrow_{\overline{e}}-\downarrow_{\mu}\uparrow_{\overline{e}}}\frac{1}{2}\mathcal{L}_{\text{eff}}\ket{\uparrow_{e}\downarrow_{\overline{\mu}}-\downarrow_{e}\uparrow_{\overline{\mu}}},$
(20)

where the factor of $1/2$ arises from the normalization of the initial and
final states. Following the argument in Clark:2003tv , the mass splitting
between two states is given by

$|\Delta M|=2|\mathcal{M}_{M}|=\frac{1}{\tau_{M}},$ (21)

where ${\tau_{M}}$ is the muonium oscillation time. A non-relativistic
reduction of the effective Lagrangian in Eq. 19 yields the following effective
potential in position space

$V_{\text{eff}}(\vec{r})=\frac{C_{\mu
e}^{2}}{64\pi^{2}v^{2}}\log{\Big{(}\frac{m_{\mu}^{2}}{M_{h}^{2}}\Big{)}}\delta^{3}(\vec{r}).$
(22)

We can assume that both $M$ and $\overline{M}$ are in the Coulombic ground
state, such that their wavefunctions are
$\phi_{100}=\exp(-\vec{r}/a_{M})/\sqrt{\pi a_{M}^{3}}$, with $a_{M}=1/\alpha
m_{\text{red}}$ being the muonium Bohr radius, and
$m_{\text{red}}=m_{\mu}m_{e}/(m_{\mu}+m_{e})\simeq m_{e}$ being the muonium
reduced mass. Therefore, the mass splitting can easily be calculated as

$|\Delta M|\simeq 2\int
d^{3}\vec{r}\phi_{100}^{*}(\vec{r})V_{\text{eff}}(\vec{r})\phi_{100}(\vec{r})\simeq\frac{C_{\mu
e}^{2}}{32\pi^{3}v^{2}a_{M}^{3}}\log{\Big{(}\frac{m_{\mu}^{2}}{M_{h}^{2}}\Big{)}},$
(23)

and the conversion rate readily follows

$P(M\rightarrow\overline{M})=\int_{0}^{\infty}dt\Gamma_{\mu}\sin^{2}{(\Delta
Mt)}e^{-\Gamma_{\mu}t}=\frac{2}{4+\Gamma_{\mu}^{2}/(\Delta M)^{2}}.$ (24)

Given the bound in Eq. (18), we find the upper limit on $|C_{\mu e}|<0.39$.

###  III.4 Constraints from the magnetic dipole moment and the $g-2$ anomaly

It was first shown by the E821 experiment at BNL Muong-2:2006rrc  and later
confirmed by the E989 experiment at Fermilab Muong-2:2021ojo ; Muong-2:2021ovs
; Muong-2:2021vma , that there is a discrepancy between the measured and
predicted Aoyama:2020ynm  magnetic dipole moment of the muon. This
discrepancy, known as the $g-2$ anomaly, currently stands at

$\Delta a_{\mu}=a_{\mu}^{\text{Exp}}-a_{\mu}^{\text{SM}}=(251\pm 59)\times
10^{-11},$ (25)

with a significance of $4.2\sigma$. On the other hand, several lattice QCD
groups have recently reported higher theoretical predictions compared to the
data-driven approach, and seem to agree with experiment Borsanyi:2020mff ;
Ce:2022kxy ; Alexandrou:2022amy . For the purposes of extracting the relevant
bounds, we shall assume that the $g-2$ anomaly exists and that it is given by
Eq. (25) above, and if future studies show that indeed the theory and
experiment agree, then the bounds are simply ignored.

Figure 5: FV contribution to the muon magnetic dipole moment through the
couplings $h^{2}\overline{l}l$. Here, $l=\\{e,\mu,\tau\\}$.

The possibility of the effective coupling $h^{2}\overline{\mu}\mu$ solving the
$g-2$ anomaly was considered in Abu-Ajamieh:2022nmt , where it was shown that
this type of coupling can accommodate the anomaly if this coupling is large
enough. It was also shown that such a deviation from the SM would point to a
scale of NP $\sim 10-18$ TeV through unitarity arguments, which can be lowered
to $\sim 5$ TeV if the Higgs couplings to $W/Z$ conform to the SM predictions.

Here we generalize the situation to FV $h^{2}\overline{l}l$ couplings. These
couplings contribute to the muon magnetic dipole moment at 2 loops as shown in
the diagrams in Figure 5. Notice that the FV case corresponds to $l=e,\tau$.
These diagrams can be evaluated using the same techniques illustrated in the
appendices and in Abu-Ajamieh:2022nmt , and they are found to provide the
following contribution to $(g-2)_{\mu}$

$\Delta a_{\mu}^{l}\simeq\frac{C_{\mu
l}^{2}}{2(4\pi)^{4}v^{2}}m_{\mu}m_{l}\Bigg{[}2\log^{2}\Big{(}\frac{m_{l}^{2}}{\Lambda^{2}}\Big{)}-\Big{(}1+\frac{2m_{\mu}}{3m_{l}}\Big{)}\log\Big(\frac{M_{h}^{2}}{\Lambda^{2}}\Big{missing})+\frac{\pi^{2}}{3}\Bigg{]},$
(26)

where the UV cutoff $\Lambda\gg M_{h}$. Setting $\Lambda=10$ TeV, the $g-2$
anomaly in Eq. (25) can be explained with the following values444The coupling
$C_{\mu\mu}$ defined here is rescaled compared to the coupling $C_{\mu 2}$
defined in Abu-Ajamieh:2022nmt . The two couplings are related as follows:
$C_{\mu 2}=\frac{vC_{\mu\mu}}{\sqrt{2}m_{\mu}}$. With this rescaling, both
results are consistent.

$\displaystyle|C_{\tau\mu}|$ $\displaystyle\simeq 0.26\pm 0.03,$ (27)
$\displaystyle|C_{\mu\mu}|$ $\displaystyle\simeq 0.79\pm 0.1,$ (28)
$\displaystyle|C_{\mu e}|$ $\displaystyle\simeq 6.34\pm 0.8.$ (29)

###  III.5 Constraints form electric dipole moment

In general, FV coupling of the form $C_{ij}$ can contribute to the Electric
Dipole Moment (EDM) of electrons and muon if such couplings are complex. In
such case, the EDM will be proportional to the imaginary parts of $C_{ij}$,
however, as we are assuming real couplings, there will be no constraints from
the EDM of the electron or the muon.

###  III.6 LEP constraints

Constraints can be obtained from LEP from the processes
$e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$. These processes are
shown in Figure 4. The s-channel involves the couplings $C_{ee}$, $C_{\mu\mu}$
and $C_{\tau\tau}$ and thus does not lead to any FV. Therefore, we ignore it
by setting these couplings to $0$. On the other hand, the t-channel involves
the FV couplings $C_{\mu e}$ and $C_{\tau e}$. Details for calculating the
loop are given in Appendix C. Using the explicit expression of the loop
integral in Eq. (68), it is a simple exercise to calculate the cross-section
of the above processes. Neglecting the masses of the initial and final state
leptons, and using $\sqrt{s}=207$ GeV555Although the COM energy of LEP is
$209$ GeV, the relevant COM energy for the processes
$e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$ quoted in
ALEPH:2006bhb  is actually $207$ GeV. and a UV cutoff $\Lambda=10^{4}$ GeV, we
find

$\sigma(e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}(\tau^{+}\tau^{-}))\simeq 2.2\times
10^{-2}C_{\mu(\tau)e}^{4}\hskip 5.69054pt\text{fb}.$ (30)

The $1\sigma$ uncertainties on
$\sigma(e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}(\tau^{+}\tau^{-}))$ are given by
$0.088$ ($0.11$) pb ALEPH:2006bhb , which can be translated into the rather
weak bounds $|C_{\mu e}|<9$ ($|C_{\tau e}|<9.52$). This is expected as these
processes are proportional to four powers of the couplings and thus cannot
compete with decay processes, which are proportional to only two powers of the
coupling. This is consistent with the case of FV from Yukawa couplings, see
for instance Harnik:2012pb .

###  III.7 Constraints from $\mu$ conversion in nuclei

The experimental searches for the conversion of $\mu\rightarrow e$ in nuclei
can be used to set limits of the leptonic effective FV couplings
$h^{2}\overline{l}l$. This process can proceed at one and two loops as shown
in Figure 6. In the notation of Kitano:2002mt , the diagram on the left
(right) is called the scalar (tensor) contribution.

Figure 6: FV contribution to $\mu\rightarrow e$ conversion in nuclei.

The scalar contribution can set limits on the coupling $C_{\mu e}$. On the
other hand, the tensor diagram can provide constraints on the combinations
$|C_{\tau\mu}||C_{\tau e}|$, $|C_{\mu\mu}||C_{\mu e}|$ or $|C_{\mu
e}||C_{ee}|$ depending on the lepton running inside the loop. However, the
tensor contribution is not expected to compete with the bounds from
$l_{i}\rightarrow l_{k}\gamma$ and therefore we neglect it here. We present
the detailed calculation in Appendix D.

From Eq. (77), the bound $\text{Br}(\mu\rightarrow e)<7\times 10^{-13}$ $@$
$90\%$ C.L. SINDRUMII:2006dvw  translates to the upper bound $|C_{\mu
e}|<0.34$. On the other hand, the Mu2e experiment is planning on improving
this limit to $\text{Br}(\mu\rightarrow e)<10^{-16}$ Kargiantoulakis:2019rjm .
This would better the bound to become $|C_{\mu e}|<4.56\times 10^{-3}$.

###  III.8 Higgs FV decays

Higgs FV decays can be used to set constraints on the leptonic couplings
$C_{ij}$. These decays proceed at one loop as shown in Figure 7. The diagram
is easily evaluated using Dimensional Regularization (DR), and the decay width
is given by

$\Gamma(h\rightarrow\overline{l}_{i}l_{j})\simeq\frac{9C_{ij}^{2}M_{h}^{5}}{4(4\pi)^{5}v^{4}},$
(31)

where we have set the renormalization scale
$\mu^{2}=M_{h}^{2}e^{-\pi/\sqrt{3}}$ and neglected the masses of the final
state. The latest bounds on these decays can be obtained from
ParticleDataGroup:2018ovx . In specific, we have the following bounds:666The
quoted bounds are $@$ $95\%$ CL. Therefore, we rescale them to $90\%$ to be
consistent with the other bounds.  $\text{Br}(h\rightarrow\mu e)<3.5\times
10^{-4}$, $\text{Br}(h\rightarrow\tau e)<6.1\times 10^{-3}$ and
$\text{Br}(h\rightarrow\tau\mu)<2.5\times 10^{-3}$. These bounds translate
into the constraints: $|C_{\mu e}|<0.25$, $|C_{\tau e}|<1.04$ and
$|C_{\tau\mu}|<0.67$ respectively. For completeness, ParticleDataGroup:2018ovx
also provides the upper bound $\text{Br}(h\rightarrow e^{+}e^{-})<1.9\times
10^{-3}$, which provides the constraint $|C_{ee}|<0.58$.

The High-Luminosity (HL) LHC is expected to yield stronger bounds on the Higgs
FV decays Banerjee:2016foh . The projected limited on the decay $h\rightarrow
e\mu$ is $3\times 10^{-5}$, which translates into a projected bound of
$|C_{\mu e}|<7.3\times 10^{-2}$. On the other hand, the project limit on the
decays $h\rightarrow\mu\tau,\tau e$ is $3\times 10^{-4}$, which leads to the
upper bound of $|C_{\tau\mu}|,|C_{\tau e}|<0.23$.

Figure 7: FV Higgs decays to $l_{i}l_{j}$.

###  III.9 Constraints from $Z\rightarrow\overline{l}l$

The excellent measurements of the $Z$ branching ratios suggest that they can
be used to extract bounds on the leptonic FV couplings. The FV couplings
$C_{ij}$ can contribute to the $Z$ decay through a process similar to the
bottom diagram of Figure 5, with the photon being replaced with $Z$, and the
external particles being leptons of the same flavor, whereas the internal
leptons being of a different flavor. Using DR, we express the corrections to
decay width of the $Z$ boson as

$\delta\text{Br}(Z\rightarrow\overline{l}_{i}l_{i})\simeq\frac{(g_{V_{l}}^{2}+g_{A_{l}}^{2})M_{Z}}{48\Gamma_{Z}}\Big{(}\frac{C_{ij}M_{h}}{16\pi^{2}v}\Big{)}^{4}\Bigg{[}\log^{2}\Big{(}\frac{M_{h}^{2}}{M_{Z}^{2}}\Big{)}-3\log\Big(\frac{M_{h}^{2}}{M_{Z}^{2}}\Big{missing})+\frac{\pi^{2}}{12}+\frac{7}{2}\Bigg{]}^{2},$
(32)

where
$g_{V_{l}}=\frac{g}{2\cos{\theta_{W}}}(T_{l_{L}}^{3}-2\sin^{2}{\theta_{W}}Q_{l})$
and $g_{A_{l}}=\frac{g}{2\cos{\theta_{W}}}(T_{l_{L}}^{3})$ are the vector and
axial couplings of the lepton $l$ to the $Z$ boson in the standard notation,
and $\Gamma_{Z}=2.4952$ GeV is the measured decay width of the $Z$. The limits
on non-FV leptonic $Z$ decays are given by ParticleDataGroup:2018ovx

$\displaystyle\text{Br}(Z\rightarrow e^{+}e^{-})$ $\displaystyle=(3363.2\pm
4.2)\times 10^{-3}\%,$ (33)
$\displaystyle\text{Br}(Z\rightarrow\mu^{+}\mu^{-})$ $\displaystyle=(3366.2\pm
6.6)\times 10^{-3}\%,$ (34)
$\displaystyle\text{Br}(Z\rightarrow\tau^{+}\tau^{-})$
$\displaystyle=(3369.8\pm 8.3)\times 10^{-3}\%.$ (35)

Given these bounds, we can extract $90\%$ C.L. constraints on the FV couplings
$C_{ij}$ by demanding that
$(\delta\text{Br})^{\text{FV}}<1.645(\delta\text{Br})^{\text{Exp}}$. Each
bound can help constrain 3 different couplings depending on the flavor of the
internal lepton, two of which are FV whereas one is flavor-conserving. Apart
from the coupling $C_{ij}$, the correction in Eq. (32) is identical for all
lepton flavors. This means that for each decay mode, the upper limit for all
three FV couplings will be identical.

The experimental limits in Eq. (33) lead to the constraints $|C_{\tau
e}|,|C_{\mu e}|,|C_{ee}|<5.62$, whereas the limits in Eq. (34) translate into
constraints $|C_{\tau\mu}|,|C_{\mu\mu}|,|C_{\mu e}|<7.04$ and the limits in
Eq. (35) yield the bounds $|C_{\tau\tau}|,|C_{\tau\mu}|,|C_{\tau e}|<7.9$.
These limits are comparable to the ones obtained from the LEP measurements
above (see subsection III.6), which is expected, as the experimental limits
shown in Eqs. (33) - (35) are essentially obtained from LEP data. However,
improved $Z$ decay measurements in future experiment, such as in the ILC
Behnke:2013xla ; can improve the these limits through its proposed ultra-
precision electroweak measurements.

###  III.10 Constraints from $Z\rightarrow\overline{l}_{i}l_{j}$

Better constraints can be obtained from the bounds on the FV decays of the $Z$
boson, because unlike the decays $Z\rightarrow\overline{l}l$ which starts at 2
loops, the decays $Z\rightarrow\overline{l}_{i}l_{j}$ start at 1 loop as shown
in Figure 8. In addition, the experimental bounds on FV final states are more
stringent compared to the flavor-conserving ones.

Figure 8: FV $Z$ decays to $\overline{l}_{i}l_{j}$.

The corrections of the diagrams in Figure 8 are easy to calculate by first
integrating out the Higgs loop then calculating the tree-level diagram. Using
DR in the $\overline{\text{MS}}$ scheme, and setting the renormalization scale
$\mu=M_{h}$, we obtain

$\delta\Gamma(Z\rightarrow\overline{l}_{i}l_{j})\simeq\frac{C_{ij}^{2}(g_{V_{l}}^{2}+g_{A_{l}}^{2})}{6(4\pi)^{5}v^{2}}\frac{m_{i}^{2}M_{Z}M_{h}^{4}}{(m_{i}^{2}-m_{j}^{2})^{2}}.$
(36)

The limits on FV leptonic $Z$ decays are given by ParticleDataGroup:2018ovx
777Here too, the bounds quoted are @ $95\%$ C.L., and we rescale them to be @
$90\%$ C.L. to be consistent with the previous results.

$\displaystyle\text{Br}(Z\rightarrow e^{\pm}\mu^{\mp})$
$\displaystyle=7.5\times 10^{-7},$ (37) $\displaystyle\text{Br}(Z\rightarrow
e^{\pm}\tau^{\mp})$ $\displaystyle=9.8\times 10^{-6},$ (38)
$\displaystyle\text{Br}(Z\rightarrow\mu^{\pm}\tau^{\mp})$
$\displaystyle=1.2\times 10^{-5}.$ (39)

Notice that for each decay, the corresponding $C_{ij}$ will have 2 possible
upper limits depending on which particle is identified as $i$ and which one is
identified as $j$. For example, for the first decay, we will have a different
bound when we identify $i$ as $e$ and $j$ as $\mu$ compared to when these
particles are flipped. As can clearly seen from Eq. (36), the strongest bound
is obtained when $i$ is identified with the heavier of the two leptons. In the
following, we quote the stronger of the two bounds. Specifically, Eqs. (37),
(38) and (39) lead to the constraints $|C_{\mu e}|<1.59\times 10^{-3}$,
$|C_{\tau e}|<9.65\times 10^{-2}$ and $|C_{\tau\mu}|<0.11$ respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: The current experimental constraints and the future projections on
the lepton FV through di-Higgs effective couplings $C_{ll^{\prime}}$. The
solid regions correspond to constraints, whereas the dashed lines represent
future projections. In (a), the green color corresponds bounds/ projections
from $\tau\rightarrow 3\mu$, whereas the blue corresponds to
$\tau\rightarrow\mu\gamma$. In (b), the green corresponds to $\mu\rightarrow
3e$, and the blue to $\mu\rightarrow e\gamma$. In (c), the green corresponds
to $\tau\rightarrow 3e$, and the blue to $\tau\rightarrow e\gamma$. In (d),
the green corresponds to $\tau^{-}\rightarrow e^{+}\mu^{-}\mu^{-}$, whereas
the blue to $\tau^{-}\rightarrow\mu^{+}\mu^{-}e^{-}$, and the red to
$\tau\rightarrow e\gamma$. In (e), the green arises from
$\tau^{-}\rightarrow\mu^{+}e^{-}e^{-}$, the blue from
$\tau^{-}\rightarrow\mu^{-}e^{+}e^{-}$, and the red from
$\tau\rightarrow\mu\gamma$. In (f), the green arises from
$\tau^{-}\rightarrow\mu^{+}\mu^{-}e^{-}$.

(a)

(b)

(c)

(d)

(e)

Figure 10: (Cont.) The current experimental constraints and the future
projections on the lepton FV through di-Higgs effective couplings
$C_{ll^{\prime}}$. The solid regions correspond to constraints, whereas the
dashed lines represent future projections. In (a), the bounds and projections
arise from $\tau^{-}\rightarrow\mu^{-}e^{+}e^{-}$, in (b) and (c) from
$\mu\rightarrow e\gamma$, in (d) from $\tau\rightarrow\mu\gamma$ and in (e)
from $\tau\rightarrow e\gamma$.

###  III.11 Fine-tuning and lepton mass corrections

None-zero $C_{ll}$ can give rise to corrections to the mass of the lepton when
the Higgs loop is closed. These corrections need to be suppressed in order to
avoid the stringent bounds on the leptons’ masses, which could lead to fine-
tuning. We can easily estimate the level of fine-tuning associated with
$C_{ll}$ as

$\frac{\delta m_{l}}{m_{l}}\sim\frac{C_{ll}M_{h}^{2}}{32\pi^{2}v^{2}}\sim
O(10^{-3})\times C_{ll},$ (40)

which is negligible for the range of $C_{ll}$ required by FV constraints.
Therefore, FV through $C_{ll}$ does not require any fine-tuning.

##  IV Quark Sector

We now turn our attention to investigating the next-to-minimal FV couplings in
the quark sector. We first discuss the constraints on the couplings $C_{ij}$
that arise from meson oscillations, then we investigate the bounds that can be
extracted from $B$-physics searches. The constraints are summarized in Table
2.

Channel | Couplings | Bounds |  $\Lambda$ (TeV)
---|---|---|---
$D^{0}$ Oscillations | $|C_{uc}|$ | $<7.73\times 10^{-4}$ | 15.3
$B^{0}_{d}$ Oscillations | $|C_{db}|$ | $<1.73\times 10^{-3}$ | 10.2
$B^{0}_{s}$ Oscillations | $|C_{sb}|$ | $<1.50\times 10^{-2}$ | 3.5
$K^{0}$ Oscillations | $|C_{sd}|$ | $<1.20\times 10^{-5}$ | 123
$R_{K^{+}}$ | $|C_{\mu\mu}|/|C_{ee}|$ | $[0.93,1.01]$ | -
$B^{0}_{d}\rightarrow\mu^{+}\mu^{-}$ | $|C_{\mu\mu}||C_{db}|$ | $<4.17\times 10^{-5}$ | 66
$B^{s}_{d}\rightarrow\mu^{+}\mu^{-}$ | $|C_{\mu\mu}||C_{sb}|$ | $<9.64\times 10^{-5}$ | 43.4
Table 2:  $90\%$ CL bounds on the FV di-Higgs effective couplings in the quark
sector and the corresponding lower limit on the scale of NP from matching to
the SMEFT. The bounds are obtained from meson oscillations and $B$-physics
searches.

###  IV.1 Constraints from meson oscillations

Constraints on the couplings $C_{ij}$ can be obtained from meson oscillations,
which include in particular $D^{0}-\bar{D}^{0}$,
$B_{s,d}^{0}-\bar{B}_{s,d}^{0}$, and $K^{0}-\bar{K}^{0}$ oscillations. These
oscillations can proceed through the di-Higgs couplings $C_{ij}$ via diagrams
identical to the ones shown in Figure 4. The effective Hamiltonian of these
diagrams can be written as UTfit:2007eik

$\mathcal{H}_{\text{eff}}=C_{2,ij}(\overline{q}_{jR}^{\alpha}q_{iL}^{\alpha})(\overline{q}_{jR}^{\beta}q_{iL}^{\beta})+C_{5,ij}(\overline{q}_{jR}^{\alpha}q_{iL}^{\beta})(\overline{q}_{jL}^{\beta}q_{iR}^{\alpha}),$
(41)

where $C_{2,ij}$ arises from integrating out the t-channel, whereas $C_{5,ij}$
arises from integrating out the s-channel in Figure 4. The detailed
calculation of these loops is presented in Appendix C. In particular, the loop
factor $V(P^{2})$ is given is Eq. (67), and in the non-relativistic limit
where we can assume that $M_{h}^{2}\gg P^{2}$, $V(P^{2})$ is approximately
given in Eq. (69). Thus, identifying the renormalization scale with the mass
of the meson $m$, we can relate $C_{2,ij}$ and $C_{5,ij}$ defined in
UTfit:2007eik  to the FV di-Higgs couplings $C_{ij}$ as follows

$C_{2,ij}\simeq
C_{5,ij}\simeq-\frac{iC_{ij}^{2}}{64\pi^{2}v^{2}}\log{\Big{(}\frac{M_{h}^{2}}{m^{2}}\Big{)}}.$
(42)

Using Eq. (42) above, we can translate the bounds on $C_{2,ij}$ and $C_{5,ij}$
presented in UTfit:2007eik  into bounds on $C_{ij}$.888The bounds presented in
UTfit:2007eik  are @ $95\%$ CL. So, we rescale them to a $90\%$ C.L. as usual
$D^{0}-\bar{D}^{0}$ oscillations place constraints on the coupling $C_{uc}$.
The stronger bound arises from $|C_{2,uc}|$ with an upper limit of $1.6\times
10^{-13}$, which translated to the constraint $|C_{uc}|<7.73\times 10^{-4}$.
$B_{d}^{0}-\bar{B}_{d}^{0}$, oscillations can set limits on the coupling
$C_{db}$, where here, the stronger of the two bounds is $|C_{5,db}|<6\times
10^{-13}$, which translated to $|C_{db}|<1.73\times 10^{-3}$. On the other
hand, $B_{s}^{0}-\bar{B}_{s}^{0}$, oscillations constrain the coupling
$C_{sb}$, with $|C_{5,sb}|<4.5\times 10^{-11}$ being the more stringent bound,
which leads to $|C_{sb}|<1.5\times 10^{-2}$. Finally, $K^{0}-\bar{K}^{0}$
oscillations place bounds on $C_{ds}$. These bounds however, only constrain
the imaginary parts of $|C_{2,sb}|$ and $|C_{5,sb}|$. Specifically, the bounds
read

$\displaystyle\text{Im}(C_{2,ds})$ $\displaystyle=[-5.1,9.3]\times 10^{-17},$
(43) $\displaystyle\text{Im}(C_{5,ds})$ $\displaystyle=[-5.2,2.8]\times
10^{-17}.$ (44)

Given Eq. (42) and our assumption that $C_{ij}$ are real, it’s not hard to see
that only the negative part of bounds in Eqs. (43) and (44) will be translated
into a bound on $C_{ds}$. In addition, it’s easy to see that $C_{2,ds}$ leads
to a stronger bound, which translates to $|C_{ds}|<1.2\times 10^{-5}$.

###  IV.2 Bounds from $B-$physics

Historically, $B$-physics attracted a lot of attention because experimental
searches revealed several discrepancies between their findings and the SM
predictions. These flavor anomalies have stirred intensive research in
$B$-physics (see London:2021lfn  for a recent review), however, recent
experimental searches seem to eliminate most of these anomalies. In
particular, the recent results from the LHCb LHCb:2022zom , reveal that lepton
universality ratios $R_{K^{+}}$ and $R_{K^{*}}$ are consistent with the SM
model. Furthermore, the latest CMS results CMS:2022mgd  also show consistency
of the decays $B_{d,s}\rightarrow\mu^{+}\mu^{-}$ with the SM model. Thus, we
can use these searches to constraint the quark FV couplings.

The lepton universality ratio $R_{K^{+}}$ is defined as

$R_{K^{+}}=\frac{\text{Br}(B^{+}\rightarrow
K^{+}\mu^{+}\mu^{-})}{\text{Br}(B^{+}\rightarrow K^{+}e^{+}e^{-})}.$ (45)

At the quark level, the decay of the $B^{+}$ meson to $K^{+}l^{+}l^{-}$ with
two leptons involves the decay
$\overline{b}\rightarrow\overline{s}l^{+}l^{-}$, which can can proceed via di-
Higgs couplings through a diagram similar to the ones in Figure 4. Given the
results in Appendix C, it’s easy to see that within our framework,
$R_{K^{+}}=\frac{|C_{\mu\mu}|^{2}}{|C_{ee}|^{2}}$. The strongest bound on
$R_{K^{+}}$ arises from the central $q^{2}$ region LHCb:2022zom , with
$R^{\text{Exp}}_{K^{+}}=0.949^{+0.048}_{-0.047}$, which translates into the
bound

$\frac{|C_{\mu\mu}|}{|C_{ee}|}=[0.93,1.01]\hskip 28.45274pt\text{@ 90\% C.L.}$
(46)

The bound is shown in part (a) of Figure 11 The Measurements of the decays
$B_{0},B_{s}\rightarrow\mu^{+}\mu^{-}$ can also furnish stringent bounds. At
the quark level, the decays $d(s)\overline{b}\rightarrow\mu^{+}\mu^{+}$
proceed via the $s$-channel in Figure 4. Utilizing the results in Appendix C,
one can show that in the limit $m_{B_{s}},m_{b}\gg m_{s},m_{\mu}$, the decay
width is given by

$\Gamma(B_{x}\rightarrow\mu^{+}\mu^{-})\simeq\frac{|C_{bx}|^{2}|C_{\mu\mu}|^{2}}{(16\pi
v)^{4}\pi}m^{3}_{B_{x}}\Big{(}1-\frac{m_{b}^{2}}{m_{B_{x}}^{2}}\Big{)}\Bigg{[}\log{\Big{(}\frac{m_{B_{x}}^{2}}{M_{h}^{2}}\Big{)}}-1\Bigg{]}^{2},$
(47)

where $x=d(s)$ for $B_{d(s)}$. The recent measurement from CMS:2022mgd  can be
used to set the following upper bound of the branching ratio of the former
decay. The bound comes out to be
$\text{Br}(B_{d}\rightarrow\mu^{+}\mu^{-})<1.61\times 10^{-10}$, which using
Eq. (47) translates into the $|C_{\mu\mu}||C_{db}|<4.17\times^{-5}$. This
bound is shown in part (b) of Figure 11.

As for the decay $B_{s}\rightarrow\mu^{+}\mu^{-}$, the SM prediction is
$\text{Br}^{\text{SM}}(B_{s}\rightarrow\mu^{+}\mu^{-})=(3.66\pm 0.14)\times
10^{-9}$, whereas the CMS measurement is given by
$\text{Br}^{\text{CMS}}(B_{s}\rightarrow\mu^{+}\mu^{-})=(3.83^{+0.38+0.24}_{-0.36-0.11})\times
10^{-9}$. Thus, we can set an upper limit on the deviation between the SM
prediction and the CMS measurement:
$\delta\text{Br}(B_{s}\rightarrow\mu^{+}\mu^{-})<9.44\times 10^{-10}$, where
we have added the theoretical and experimental uncertainties in quadrature and
rescaled the bound to be $@$ $90\%$ CL. Using Eq. (47), this translates to
$|C_{\mu\mu}||C_{sb}|<9.64\times 10^{-5}$. The bound is shown in part (c) of
Figure 11.

(a)

(b)

(c)

Figure 11: The experimental bounds on the FV di-Higgs effective couplings
extracted from $B$-physics. (a) corresponds to the bounds from the measurement
of $R_{K^{+}}$, (b) and (c) correspond to the bounds obtained from the
measurement of $B_{d},B_{s}\rightarrow\mu^{+}\mu^{-}$ respectively. The solid
blue region is excluded.

##  V Matching to the SMEFT and the Scale of New Physics

Finally in this section, we show how our framework matches to the SMEFT, then
use the upper bounds on the FV Wilson coefficients to set lower limits on the
corresponding scale of NP. Working in the Warsaw basis Grzadkowski:2010es ,
there is only one class of operators at dimension-6 that contributes to the FV
di-Higgs couplings, which has the form $\psi^{2}H^{3}$. There are 3 operator
categories in $\psi^{2}H^{3}$

$\mathcal{L}_{\psi^{2}H^{3}}=\frac{C^{l}_{ij}}{\Lambda^{2}}(H^{\dagger}H)(\overline{l}_{i}He_{j})+\frac{C^{u}_{ij}}{\Lambda^{2}}(H^{\dagger}H)(\overline{q}_{i}\tilde{H}u_{j})+\frac{C^{d}_{ij}}{\Lambda^{2}}(H^{\dagger}H)(\overline{q}_{i}Hd_{j})+h.c.,$
(48)

which should be matched to the operators in eq. 6. The matching is identical
for all of the operators and is fairly straightforward: We simply plug the
Higgs doublet in eq. 48, then we match the $h^{2}$ term to eq. 6. setting
$C_{ij}^{\text{SMEFT}}=1$, we find

$C_{ij}=\frac{3v^{2}}{\Lambda^{2}}.$ (49)

Eq. 49 can be used to set a lower limit on the scale of NP $\Lambda$ from the
upper bounds on $C_{ij}$. We present these bounds in Tables 1 and 2. In the
lepton sector, we can see that lower bounds ranges between $\sim 1-10$ TeV. On
the other hand, the stronger bounds in the quark sector lead to much higher
scales $\Lambda$, ranging from a few TeV to $\sim 123$ TeV.

##  VI Conclusions

In this paper, we employed a completely model-independent bottom-up EFT to
investigate FV in the quark and the lepton interactions with the Higgs. In
this approach we dubbed the WSD, we did not resort to any power expansion, and
instead listed the most general FV interactions. This approach is a
generalization of the one introduced in Chang:2019vez ; Abu-Ajamieh:2020yqi ;
Abu-Ajamieh:2021egq ; Abu-Ajamieh:2022ppp ; Abu-Ajamieh:2021vnh  to the FV
case.

Unlike previous studies on FV in the Higgs sector which focused on FV Yukawa
couplings $Y_{ij}\neq\sqrt{2}m_{i}\delta_{ij}/v$. In this paper, we focused on
the next-to-minimal FV that arises from the di-Higgs effective couplings of
the form $h^{2}\overline{f}f$, and assumed that the Yukawa couplings are equal
to the SM predictions. To the best our knowledge, this is the first time
constraints are set on these types of couplings.

In the lepton sector, we investigated the bounds on the FV di-Higgs couplings
that arise form $l\rightarrow l_{1}l_{2}l_{3}$ decays, $l_{i}\rightarrow
l_{k}\gamma$ decays, muonium oscillations, the $g-2$ anomaly, LEP searches,
$\mu\rightarrow e$ conversion in nuclei, leptonic FV Higgs decays, and from
both flavor-conserving and FV $Z$ decays. We have set upper limits on both
individual effective couplings and products of the various couplings, and
found that these bounds in general vary between $\sim O(1)$ down to $\sim
O(10^{-5})$. In addition, we utilized the projections of some future
experiments, such as Belle II, the Mu2e experiment and the HL-LHC in order to
find future projections on some of these couplings, and found that these
bounds can be improved by roughly a factor ranging between a few and two
orders of magnitude. The bounds and future projections are summarized in Table
1 and Figures 9 and 10. On the other hand, bounds on the FV di-Higgs couplings
in the quark sector were obtained from meson oscillations and from $B$-physics
searches, and these range between $\sim O(10^{-2})$ down to $\sim O(10^{-5})$.
These bounds are summarized in Table 2 and Figure 11.

We have shown how our approach can be mapped to the SMEFT and have shown the
scale of NP that corresponds to the upper limits on the FV Wilson
coefficients. We saw that the scale of NP ranges between $\sim 1-10$ TeV in
the lepton sector, and between a few TeV up to $\sim 123$ TeV in the quark
sector. We believe that measuring the di-Higgs effective couplings, whether
flavor-conserving or FV, is of particular importance and should receive
adequate attention in the LHC searches and other low-energy experiments. The
proposed muon collider would be an interesting laboratory where these
couplings can be probed.

## Acknowledgments

The work of F.A. is supported by the C.V. Raman fellowship from CHEP at IISc.
S.K.V. thanks SERB Grant CRG/2021/007170 ”Tiny Effects from Heavy New Physics”
and Matrics grant MTR/2022/000255, ”Theoretical Aspects of Certain Physics
Beyond Standard Models” from the Department of Science and Technology,
Government of India.

##  Appendix A The decay width of $l\rightarrow l_{1}l_{2}l_{3}$

The matrix element of the decay shown in Figure 2 is given by

$\mathcal{M}_{123}=-i\frac{C_{ll_{1}}C_{l_{2}l_{3}}}{2v^{2}}\overline{u}(q_{1})\int\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{(k-p+q_{1})^{2}-M^{2}}\frac{1}{k^{2}-M^{2}}u(p)\overline{u}(q_{3})v(q_{2}),$
(50)

where $M$ is the mass of the particle in the loop. The loop is logarithmically
divergent and needs regularization. We use DR to perform the momentum integral

$\mathcal{M}_{123}=\frac{C_{ll_{1}}C_{l_{2}l_{3}}}{32\pi^{2}v^{2}}\overline{u}(q_{1})u(p)\overline{u}(q_{3})v(q_{2})\Gamma(2-\frac{d}{2})\Big{(}\frac{\mu^{2}}{M^{2}}\Big{)}^{2-\frac{d}{2}}\int_{0}^{1}dx\Big{[}1+x(x-1)(m^{2}+m_{1}^{2}-2p.q_{1})/M^{2}\Big{]}^{\frac{d}{2}-2},$
(51)

where $m,m_{1}$ are the masses of $l,l_{1}$ respectively and $\mu$ is the
renormalization scale. Before we perform the integral over the Feynman
parameter, we notice that in the limit $M\gg m,m_{1}$ applicable in our case,
we can drop the masses $m$ and $m_{1}$ from the integral. In addition, in the
rest frame of the decaying particle, $p.q_{1}=mE_{1}$, with $E_{1}$ being the
energy of $l_{1}$. Given that the upper limit of $E_{1}$ is $m/2$, we can drop
that term as well. Therefore, in the limit $M\gg m,m_{1}$ the integral in Eq.
(52) becomes trivial. Setting $d=4-2\epsilon$ and using the
$\overline{\text{MS}}$ scheme, the regularized matrix element reads

$\mathcal{M}_{123}=-\frac{C_{ll_{1}}C_{l_{2}l_{3}}}{32\pi^{2}v^{2}}\log{\Big{(}\frac{M^{2}}{m^{2}}\Big{)}}\overline{u}(q_{1})u(p)\overline{u},(q_{3})v(q_{2}),$
(52)

where have set $\mu^{2}=m^{2}$. Before we calculate the decay width, we point
out that depending on the decay, there could be either one or two Feynman
diagrams. For example, $\tau\rightarrow 3\mu$ obviously involves only one
Feynman diagram, i.e. $l=\tau$, and $l_{1},l_{2},l_{3}=\mu$. On the other
hand, a process like $\tau^{-}\rightarrow\mu^{+}\mu^{-}e^{-}$ involves two
diagram: the first with $l_{1}=e^{-}$, $l_{2}=\mu^{+}$, and $l_{3}=\mu^{-}$,
and the second with the $l_{1}$ and $l_{3}$ interchanged. The matrix elements
of the two diagrams should be added together, with the appropriate Fermi-Dirac
statistics taken into consideration. Here we show the decay width of processes
with only one Feynman diagram. Generalizing to processes with two Feynman
diagrams is straightforward.

Since the decays we are interested in are $\tau\rightarrow 3\mu$,
$\tau\rightarrow 3e$, $\tau\rightarrow\mu\mu e$, $\tau\rightarrow\mu ee$ and
$\mu\rightarrow 3e$, in all cases we have $m\gg m_{1},m_{2},m_{3}$. Thus, we
can treat the final states as massless. This simplifies the phase space
integral greatly and the final result reads

$\Gamma(l\rightarrow
l_{1}l_{2}l_{3})=\frac{m^{5}}{v^{4}}\Bigg{[}\frac{C_{ll_{1}}C_{l_{2}l_{3}}}{512\pi^{3}\sqrt{6\pi}}\log{\Big{(}\frac{M_{h}^{2}}{m_{l}^{2}}\Big{)}}\Bigg{]}^{2}.$
(53)

##  Appendix B Calculating the 2-loop diagram of $l_{i}\rightarrow
l_{k}\gamma$

Here we show the general calculation of the 2-loop diagram in Figure 3. This
diagram is the leading contribution to the decays $\tau\rightarrow\mu\gamma$,
$\tau\rightarrow e\gamma$ and $\mu\rightarrow e\gamma$. Notice that in each
case, the inner particle $j$ could be either $\tau$, $\mu$, or $e$, which
leads to different structures of the matrix element with the corresponding
effective couplings $C_{ij}$ and $C_{jk}$. We can write the matrix element as

$\mathcal{M}_{ijk}=\frac{eC_{ij}C_{jk}}{2v^{2}}\overline{u}(p-q,m_{k})I_{ijk}u(p,m_{i})\epsilon_{\mu}^{*}(q),$
(54)

where the two-loop momentum integral is given by

$I_{i,j,k}=\int\frac{d^{4}k_{1}}{(2\pi)^{4}}\frac{(\not{k}_{1}-\not{q}+m_{j})\gamma^{\mu}(\not{k}_{1}+m_{j})}{[(k_{1}-q)^{2}-m_{j}^{2}][k_{1}^{2}-m_{j}^{2}]}\int\frac{d^{4}k_{2}}{(2\pi)^{4}}\frac{1}{[(p-k_{1}-k_{2})^{2}-M^{2}][k_{2}^{2}-M^{2}]}.$
(55)

We can perform the integral over $k_{2}$ first, then combine the results with
the remaining integral over $k_{1}$, and finally perform the momentum integral
over $k_{1}$. Using DR, we find the following general form of the matrix
element

$\displaystyle\mathcal{M}_{ijk}$
$\displaystyle=\frac{eC_{ij}C_{jk}}{(4\pi)^{4}v^{2}}\Gamma(4-d)(-1)^{5-\frac{d}{2}}\Big{(}\frac{4\pi\mu^{2}}{M_{h}^{2}}\Big{)}^{4-d}\overline{u}(p-q,m_{k})u(p,m_{i})(p\cdot\epsilon^{*})$
(56)
$\displaystyle\times\int_{0}^{1}dx\int_{0}^{1}dy\int_{0}^{1}dzz^{1-\frac{d}{2}}(1-z)\frac{\Big{[}(a+b-1)(a\hskip
1.42262ptm_{i}+m_{j})-b(a\hskip
1.42262ptm_{k}+m_{j})\Big{]}}{\alpha^{\frac{d}{2}}\beta^{4-d}},$

where $M_{h}$ is the mass of the Higgs, $p^{\mu}$ the momentum of the initial
state lepton, and the functions $\alpha$, $\beta$, $a$ and $b$ are given by

$\displaystyle a$ $\displaystyle=\frac{xz(x-1)}{\alpha},$ (57) $\displaystyle
b$ $\displaystyle=\frac{y(z-1)}{\alpha},$ (58) $\displaystyle\alpha$
$\displaystyle=(x^{2}-x+1)z-1,$ (59) $\displaystyle\beta$
$\displaystyle=-z+\frac{xz(x-1)(y-1)(z-1)}{\alpha}\frac{m_{i}^{2}}{M_{h}^{2}}+(z-1)\frac{m_{j}^{2}}{M_{h}^{2}}-\frac{xyz(x-1)(z-1)}{\alpha}\frac{m_{k}^{2}}{M_{h}^{2}}.$
(60)

The integrals in Eq. (56) are badly divergent and care is needed to regularize
them. In addition, it is not possible to evaluate them exactly for any general
particles $i,j$ and $k$. Thus, was need to approximate them by assuming
$M_{h}\gg m_{\tau}\gg m_{\mu}\gg m_{e}$, and only keep the lepton with the
largest mass in each decay. Notice that in Eq. (B7), although $M\gg
m_{i,j,k}$, we need to keep the term with the largest lepton mass to keep the
integral IR finite. Therefore, evaluating Eq. (56) will depend on what the
particles $i$, $j$ and $k$ are. In order to set upper limits on the FV
couplings $C_{ij}$, we treat each case separately. For example, for the
process $\tau\rightarrow\mu$, we could have $j={\tau,\mu,e}$ running in the
loop. This furnished 9 distinct processes in total to consider. Here we show a
sample calculation, then quote the results for the rest of the process.

Consider the process $\tau\rightarrow e\gamma$ with $\mu$ in the loop. We
denote the corresponding matrix element by $\mathcal{M}_{\tau\mu e}$, with
$m_{i}=m_{\tau}$, $m_{j}=m_{\mu}$ and $m_{k}=m_{e}$. Dropping $m_{\mu},m_{e}$,
the integral in Eq. (56) simplifies to

$\mathcal{M}_{\tau\mu e}\simeq
m_{\tau}\int_{0}^{1}dx\int_{0}^{1}dy\int_{0}^{1}dzz^{1-\frac{d}{2}}(1-z)\frac{a(a+b-1)}{\alpha^{\frac{d}{2}}\beta^{4-d}},$
(61)

with

$\beta\simeq-z+\frac{xz(x-1)(y-1)(z-1)}{\alpha}\frac{m_{\tau}^{2}}{M^{2}}.$
(62)

The integral in Eq. (61) is still divergent. So, in order to regularize it, we
use the method described in Peskin:1995ev . First, we define the function

$f(z)\equiv(1-z)\frac{a(a+b-1)}{\alpha^{\frac{d}{2}}\beta^{4-d}}.$ (63)

then isolate the divergence by splitting the integral over $z$ as follows

$f(z)=\int_{0}^{1}dx\int_{0}^{1}dy\Bigg{[}\int_{0}^{1}dzz^{1-\frac{d}{2}}f(0)+\int_{0}^{1}dzz^{1-\frac{d}{2}}\Big{(}f(z)-f(0)\Big{)}\Bigg{]}=-\frac{1}{6}.$
(64)

Plugging Eq. (64) in Eq. (56), then setting $d=4-2\epsilon$ and using the
$\overline{\text{MS}}$ scheme, we arrive at final answer

$\mathcal{M}_{\tau\mu e}\simeq-\frac{eC_{\tau\mu}C_{\mu
e}m_{\tau}}{6(4\pi)^{4}v^{2}}\overline{u}_{e}(p-q)u_{\tau}(p)(p\cdot\epsilon^{*})\log{\Big{(}\frac{M_{h}^{2}}{m_{\mu}^{2}}\Big{)}},$
(65)

where we have set the renormalization scale $\mu^{2}=m_{\mu}^{2}$ in the
logarithm.

##  Appendix C $f_{i}\overline{f}_{j}\rightarrow f_{k}\overline{f}_{l}$
scattering

Here we show how to calculate the matrix element of the the scattering
$f_{i}\overline{f}_{j}\rightarrow f_{k}\overline{f}_{l}$, which will be used
to find the bounds from LEP, muonium-antimuonium oscillations and meson
oscillation. At 1-loop, the scattering proceeds through the s- and t-channels
as in Figure (3). The matrix element is given by

$\displaystyle i\mathcal{M}$
$\displaystyle=i\mathcal{M}_{s}+i\mathcal{M}_{t},$
$\displaystyle=\frac{C_{ij}C_{kl}}{4v^{2}}\overline{u}_{k}(k_{1})v_{l}(k_{2})\overline{v}_{j}(p_{2})u_{i}(p_{1})V(P_{s}^{2})-\frac{C_{ik}C_{jl}}{4v^{2}}\overline{u}_{k}(k_{1})u_{i}(p_{1})\overline{v}_{j}(p_{2})v_{l}(k_{2})V(P_{t}^{2}),$
(66)

where $P_{s}=p_{1}+p_{2}$, $P_{t}=p_{1}-k_{1}$, and $p_{1,2}$ ($k_{1,2}$) are
the initial (final) momenta. The loop integral is given by

$V(P^{2})=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{(k+P)^{2}-M^{2}}\frac{1}{k^{2}-M^{2}}.$
(67)

The integral in Eq. (67) is logarithmically divergent and needs
regularization. The suitable choice of regularization will depend on the type
of process at hand. In high energy scattering like in LEP, using a UV cutoff
is more appropriate. Evaluating the integral using a UV cutoff $\Lambda$, the
final result can be approximated by

$V(P^{2})\simeq\frac{i}{16\pi^{2}}\Bigg{(}1+\log{\Big{(}\frac{\Lambda^{2}}{M^{2}}\Big{)}}+\sqrt{1-\frac{4M^{2}}{P^{2}}}\log{\Bigg{[}\frac{\sqrt{1-4M^{2}/P^{2}}-1}{\sqrt{1-4M^{2}/P^{2}}+1}\Bigg{]}}\Bigg{)}.$
(68)

On the other hand, in the non-relativistic limit suitable for $M-\overline{M}$
and meson oscillation, it is more suitable to evaluate the integral using DR.
In the $\overline{\text{MS}}$ scheme, the integral evaluates to

$V(P^{2})\simeq\frac{i}{16\pi^{2}}\log{\Big{(}\frac{\mu^{2}}{M_{h}^{2}}\Big{)}},$
(69)

where $\mu$ is the renormalization scale. Notice that in the non-relativistic
limit $M^{2}\gg P^{2}$, Eq. (69) can be obtained from Eq. (68) by taking the
limit $P^{2}\rightarrow 0$ and then setting $\Lambda^{2}=e\mu^{2}$.

##  Appendix D Detailed calculation of $\mu\rightarrow e$ conversion in nuclei

The most general effective Lagrangian can be expressed as Kitano:2002mt

$\mathcal{L}_{\text{eff}}=c_{L}\frac{e}{8\pi^{2}}m_{\mu}(\overline{e}\sigma^{\mu\nu}P_{L}\mu)F_{\mu\nu}-\frac{1}{2}\sum_{q}\Big{[}g^{q}_{LS}(\overline{e}P_{R}\mu)(\overline{q}q)+g^{q}_{LP}(\overline{e}P_{R}\mu)(\overline{q}\gamma_{5}q)\\\
+g^{q}_{LV}(\overline{e}\gamma^{\mu}P_{L}\mu)(\overline{q}\gamma_{\mu}q)+g^{q}_{LA}(\overline{e}\gamma^{\mu}P_{L}\mu)(\overline{q}\gamma_{\mu}\gamma_{5}q)+\frac{1}{2}g^{q}_{LT}(\overline{e}\sigma^{\mu\nu}P_{R}\mu)(\overline{q}\sigma_{\mu\nu}q)\Big{]}+(L\leftrightarrow
R),$ (70)

where the sum is over all quarks. Here, the first term expresses the
contributions arising from the magnetic dipole operators as in the bottom
diagram of Figure 3. On the other hand, the terms inside the square brackets
refer to the scalar, pseudo-scalar, vector, pseudo-vector and tensor
contributions, respectively. As shown in Figure 6, only the scalar and tensor
contributions are non-vanishing. Furthermore, the tensor contribution is
expected to be small and the bounds are not expected to compete with those
from $l_{i}\rightarrow l_{k}\gamma$, therefore we neglect it as well.

The scalar contribution $g^{q}_{LS}$ and $g^{q}_{RS}$, are shown in the left
diagram of Figure 6. They can be calculated by integrating out the loop in the
non-relativistic limit and at vanishing momentum transfer, yielding

$g^{q}_{LS}=g^{q}_{RS}\equiv g^{q}_{S}=\frac{3\sqrt{2}C_{\mu
e}Y_{q}^{2}m_{N}}{64\pi^{2}vM_{h}^{2}}.$ (71)

where $Y_{q}$ is the quark Yukawa coupling and $m_{N}$ is the mass of the
nucleon. The $\mu\rightarrow e$ conversion rate receives contributions from
protons and neutrons and can be expressed as Kitano:2002mt

$\Gamma(\mu\rightarrow
e)=|\tilde{g}^{(p)}_{S}S^{(p)}+\tilde{g}^{(n)}_{S}S^{(n)}|^{2},$ (72)

where

$\tilde{g}^{(p)}_{S}=\sum_{q}g^{q}_{S}\frac{m_{p}}{m_{q}}f^{(q,p)},\hskip
56.9055pt\tilde{g}^{(n)}_{S}=\sum_{q}g^{q}_{S}\frac{m_{n}}{m_{q}}f^{(q,n)},$
(73)

where the nucleon matrix elements
$f^{(q,N)}\equiv\bra{N}m_{q}\overline{q}q\ket{N}/m_{N}$. These nucleon matrix
elements were calculated in Ellis:2008hf  but using an older value for the
nucleon sigma term $\Sigma_{\pi N}=64$ MeV. Using the updated value of $59.6$
MeV Gupta:2021ahb 999In Harnik:2012pb , the nucleon matrix elements where
calculated using the then latest value of $\Sigma_{\pi N}=55$ MeV, however,
there is an error in their equation A19. In particular, $f^{(u,n)}=0.018\neq
f^{(d,p)}$, and $f^{(d,n)}=0.043\neq f^{(u,p)}$. All other values were
correctly calculated for $\Sigma_{\pi N}=55$ MeV., the nucleon matrix elements
for the light quarks are given by

$\displaystyle f^{(u,p)}\simeq 0.022,\hskip 28.45274ptf^{(d,p)}\simeq
0.038,\hskip 28.45274ptf^{(s,p)}\simeq 0.342,$ (74) $\displaystyle
f^{(u,n)}\simeq 0.018,\hskip 28.45274ptf^{(d,n)}\simeq 0.049,\hskip
28.45274ptf^{(s,n)}\simeq 0.342,$ (75)

whereas the contribution for the heavy quarks is obtained from

$f^{(c,N)}=f^{(b,N)}=f^{(t,N)}=\frac{2}{27}\Big{(}1-\sum_{q=u,d,s}f^{(q,N)}\Big{)}\simeq
0.044,$ (76)

for both the neutron and proton. The coefficients $S^{(p)}$, $S^{(n)}$ are the
overlap integrals of the electron, muon and nuclear wavefunctions for the
proton and neutron respectively. They are tabulated for a variety of target
materials in Kitano:2002mt . According to SINDRUMII:2006dvw , gold provides
the strongest bound on the conversion rate

$\text{Br}^{\text{Au}}(\mu\rightarrow e)=\Bigg{[}\frac{\Gamma(\mu\rightarrow
e)}{\Gamma^{\mu}_{\text{Capture}}}\Bigg{]}_{\text{Au}}<7\times 10^{-13}\hskip
14.22636pt@\hskip 2.84526pt90\%\hskip 5.69054pt\text{C.L.},$ (77)

and we find from Kitano:2002mt  that
$\Gamma^{\text{Au}}_{\text{Capture}}=13.07\times 10^{6}\hskip
2.84526pt\text{s}^{-1}$. In addition, the overlap coefficients for gold are
given by $S^{(p)}=0.0614$ and $S^{(n)}=0.0918$ in units of $m_{\mu}^{5/2}$. On
the other hand, the Mu2e experiment is projected to improve the measurement of
the conversion rate by roughly 3 orders of magnitude through utilizing
aluminium as its stopping material. More specifically, the projected bound of
the Mu2e experiment is given by Kargiantoulakis:2019rjm

$\text{Br}^{\text{Al}}(\mu\rightarrow e)=\Bigg{[}\frac{\Gamma(\mu\rightarrow
e)}{\Gamma^{\mu}_{\text{Capture}}}\Bigg{]}_{\text{Al}}<10^{-16}\hskip
14.22636pt@\hskip 2.84526pt90\%\hskip 5.69054pt\text{C.L.},$ (78)

and we have $\Gamma^{\text{Al}}_{\text{Capture}}=0.7054\times 10^{6}\hskip
2.84526pt\text{s}^{-1}$, and the overlap coefficients for aluminium are given
by $S^{(p)}=0.0155$ and $S^{(n)}=0.0167$ in units of $m_{\mu}^{5/2}$.

## References

  * (1) G. D’Ambrosio, G. F. Giudice, G. Isidori and A. Strumia, Nucl. Phys. B 645 (2002), 155-187 doi:10.1016/S0550-3213(02)00836-2 [arXiv:hep-ph/0207036 [hep-ph]].
  * (2) A. Dery, A. Efrati, Y. Hochberg and Y. Nir, “What if $BR(h\to\mu\mu)$ $BR(h\to\tau\tau)\not=m_{\mu}^{2}/m_{\tau}^{2}$?,” JHEP 05, 039 (2013) hep-ph/1302.3229.
  * (3) R. Harnik, J. Kopp and J. Zupan, “Flavor Violating Higgs Decays,” JHEP 03, 026 (2013) hep-ph/1209.1397.
  * (4) J. de Blas, J. C. Criado, M. Perez-Victoria and J. Santiago, “Effective description of general extensions of the Standard Model: the complete tree-level dictionary,” JHEP 03 (2018), 109 hep-ph/1711.10391.
  * (5) F. del Aguila, M. Perez-Victoria and J. Santiago, “Effective description of quark mixing,” Phys. Lett. B 492 (2000), 98-106 hep-ph/hep-ph/0007160.
  * (6) B. Batell, S. Gori and L. T. Wang, “Higgs Couplings and Precision Electroweak Data,” JHEP 01 (2013), 139 hep-ph/1209.6382.
  * (7) C. Y. Chen, S. Dawson and E. Furlan, “Vectorlike fermions and Higgs effective field theory revisited,” Phys. Rev. D 96 (2017) no.1, 015006 hep-ph/1703.06134.
  * (8) L. Silvestrini and M. Valli, “Model-independent Bounds on the Standard Model Effective Theory from Flavour Physics,” Phys. Lett. B 799 (2019), 135062 hep-ph/1812.10913.
  * (9) S. Descotes-Genon, A. Falkowski, M. Fedele, M. González-Alonso and J. Virto, “The CKM parameters in the SMEFT,” JHEP 05 (2019), 172 hep-ph/1812.08163.
  * (10) J. Aebischer, J. Kumar, P. Stangl and D. M. Straub, “A Global Likelihood for Precision Constraints and Flavour Anomalies,” Eur. Phys. J. C 79 (2019) no.6, 509 hep-ph/1810.07698.
  * (11) A. Greljo and A. Palavrić, “Leading directions in the SMEFT,” hep-ph/2305.08898.
  * (12) A. Greljo, A. Palavrić and A. E. Thomsen, “Adding Flavor to the SMEFT,” JHEP 10 (2022), 010 hep-ph/2203.09561.
  * (13) S. Bruggisser, R. Schäfer, D. van Dyk and S. Westhoff, “The Flavor of UV Physics,” JHEP 05 (2021), 257 hep-ph/2101.07273.
  * (14) R. Aoude, T. Hurth, S. Renner and W. Shepherd, “The impact of flavour data on global fits of the MFV SMEFT,” JHEP 12 (2020), 113 hep-ph/2003.05432.
  * (15) T. Hurth, S. Renner and W. Shepherd, “Matching for FCNC effects in the flavour-symmetric SMEFT,” JHEP 06 (2019), 029 hep-ph/1903.00500.
  * (16) J. Aebischer and J. Kumar, “Flavour violating effects of Yukawa running in SMEFT,” JHEP 09 (2020), 187 hep-ph/2005.12283.
  * (17) D. A. Faroughy, G. Isidori, F. Wilsch and K. Yamamoto, “Flavour symmetries in the SMEFT,” JHEP 08 (2020), 166 hep-ph/2005.05366.
  * (18) F. Abu-Ajamieh and S. K. Vempati, “Can the Higgs Still Account for the g-2 Anomaly?,” hep-ph/2209.10898.
  * (19) S. Chang and M. A. Luty, “The Higgs Trilinear Coupling and the Scale of New Physics,” JHEP 03, 140 (2020) hep-ph/1902.05556.
  * (20) F. Abu-Ajamieh, S. Chang, M. Chen and M. A. Luty, “Higgs coupling measurements and the scale of new physics,” JHEP 07, 056 (2021) hep-ph/2009.11293.
  * (21) F. Abu-Ajamieh, “The scale of new physics from the Higgs couplings to $\gamma$$\gamma$ and $\gamma$Z,” JHEP 06, 091 (2022) hep-ph/2112.13529.
  * (22) F. Abu-Ajamieh, “The scale of new physics from the Higgs couplings to gg,” Phys. Lett. B 833, 137389 (2022) hep-ph/2203.07410.
  * (23) F. Abu-Ajamieh, “Model-independent Veltman condition, naturalness and the little hierarchy problem,” Chin. Phys. C 46, no.1, 013101 (2022) hep-ph/2101.06932.
  * (24) B. Grinstein and M. Trott, “A Higgs-Higgs bound state due to new physics at a TeV,” Phys. Rev. D 76, 073002 (2007) hep-ph/0704.1505.
  * (25) Z. N. Zhang, H. B. Zhang, J. L. Yang, S. M. Zhao and T. F. Feng, “Higgs boson decays with lepton flavor violation in the $B-L$ symmetric SSM,” Phys. Rev. D 103 (2021) no.11, 115015 hep-ph/2105.09799.
  * (26) A. Vicente, “Higgs lepton flavor violating decays in Two Higgs Doublet Models,” Front. in Phys. 7 (2019), 174 hep-ph/1908.07759.
  * (27) Y. Soreq, H. X. Zhu and J. Zupan, “Light quark Yukawa couplings from Higgs kinematics,” JHEP 12 (2016), 045 hep-ph/1606.09621.
  * (28) M. Buschmann, J. Kopp, J. Liu and X. P. Wang, “New Signatures of Flavor Violating Higgs Couplings,” JHEP 06 (2016), 149 hep-ph/1601.02616.
  * (29) M. Tanabashi et al. [Particle Data Group], “Review of Particle Physics,” Phys. Rev. D 98, no.3, 030001 (2018)
  * (30) T. Aushev, W. Bartel, A. Bondar, J. Brodzicka, T. E. Browder, P. Chang, Y. Chao, K. F. Chen, J. Dalseno and A. Drutskoy, et al. “Physics at Super B Factory,” hep-ex/1002.5012.
  * (31) L. Aggarwal et al. [Belle-II], “Snowmass White Paper: Belle II physics reach and plans for the next decade and beyond,” hep-ex/2207.06307.
  * (32) L. Calibbi and G. Signorelli, “Charged Lepton Flavour Violation: An Experimental and Theoretical Introduction,” Riv. Nuovo Cim. 41, no.2, 71-174 (2018) hep-ph/1709.00294.
  * (33) S. Banerjee, “Searches for Lepton Flavor Violation in Tau Decays at Belle II,” Universe 8, no.9, 480 (2022) hep-ex/2209.11639.
  * (34) J. Beringer et al. [Particle Data Group], “Review of Particle Physics (RPP),” Phys. Rev. D 86, 010001 (2012)
  * (35) L. Willmann, P. V. Schmidt, H. P. Wirtz, R. Abela, V. Baranov, J. Bagaturia, W. H. Bertl, R. Engfer, A. Grossmann and V. W. Hughes, et al. “New bounds from searching for muonium to anti-muonium conversion,” Phys. Rev. Lett. 82, 49-52 (1999) hep-ex/hep-ex/9807011.
  * (36) T. E. Clark and S. T. Love, “Muonium - anti-muonium oscillations and massive Majorana neutrinos,” Mod. Phys. Lett. A 19, 297-306 (2004) hep-ph/hep-ph/0307264.
  * (37) G. W. Bennett et al. [Muon g-2], “Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL,” Phys. Rev. D 73, 072003 (2006) hep-ex/hep-ex/0602035.
  * (38) B. Abi et al. [Muon g-2], “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm,” Phys. Rev. Lett. 126, no.14, 141801 (2021) hep-ex/2104.03281.
  * (39) T. Albahri et al. [Muon g-2], “Magnetic-field measurement and analysis for the Muon $g-2$ Experiment at Fermilab,” Phys. Rev. A 103, no.4, 042208 (2021) hep-ex/2104.03201.
  * (40) T. Albahri et al. [Muon g-2], “Measurement of the anomalous precession frequency of the muon in the Fermilab Muon $g-2$ Experiment,” Phys. Rev. D 103, no.7, 072002 (2021) hep-ex/2104.03247.
  * (41) 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. “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887, 1-166 (2020) hep-ph/2006.04822.
  * (42) 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. “Leading hadronic contribution to the muon magnetic moment from lattice QCD,” Nature 593, no.7857, 51-55 (2021) hep-lat/2002.12347.
  * (43) 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. “Window observable for the hadronic vacuum polarization contribution to the muon g-2 from lattice QCD,” Phys. Rev. D 106, no.11, 114502 (2022) hep-lat/2206.06582.
  * (44) C. Alexandrou, S. Bacchio, P. Dimopoulos, J. Finkenrath, R. Frezzotti, G. Gagliardi, M. Garofalo, K. Hadjiyiannakou, B. Kostrzewa and K. Jansen, et al. “Lattice calculation of the short and intermediate time-distance hadronic vacuum polarization contributions to the muon magnetic moment using twisted-mass fermions,” hep-lat/2206.15084.
  * (45) J. Alcaraz et al. [ALEPH, DELPHI, L3, OPAL and LEP Electroweak Working Group], “A Combination of preliminary electroweak measurements and constraints on the standard model,” hep-ex/hep-ex/0612034.
  * (46) R. Kitano, M. Koike and Y. Okada, “Detailed calculation of lepton flavor violating muon electron conversion rate for various nuclei,” Phys. Rev. D 66, 096002 (2002) [erratum: Phys. Rev. D 76, 059902 (2007)] hep-ph/hep-ph/0203110.
  * (47) W. H. Bertl et al. [SINDRUM II], “A Search for muon to electron conversion in muonic gold,” Eur. Phys. J. C 47, 337-346 (2006)
  * (48) M. Kargiantoulakis, “A Search for Charged Lepton Flavor Violation in the Mu2e Experiment,” Mod. Phys. Lett. A 35, no.19, 2030007 (2020) hep-ex/2003.12678.
  * (49) S. Banerjee, B. Bhattacherjee, M. Mitra and M. Spannowsky, “The Lepton Flavour Violating Higgs Decays at the HL-LHC and the ILC,” JHEP 07, 059 (2016) hep-ph/1603.05952.
  * (50) T. Behnke, J. E. Brau, B. Foster, J. Fuster, M. Harrison, J. M. Paterson, M. Peskin, M. Stanitzki, N. Walker and H. Yamamoto, “The International Linear Collider Technical Design Report - Volume 1: Executive Summary,” physics.acc-ph/1306.6327.
  * (51) M. Bona et al. [UTfit], “Model-independent constraints on $\Delta F=2$ operators and the scale of new physics,” JHEP 03, 049 (2008) hep-ph/0707.0636.
  * (52) D. London and J. Matias, “$B$ Flavour Anomalies: 2021 Theoretical Status Report,” Ann. Rev. Nucl. Part. Sci. 72, 37-68 (2022) hep-ph/2110.13270.
  * (53) [LHCb], “Measurement of lepton universality parameters in $B^{+}\to K^{+}\ell^{+}\ell^{-}$ and $B^{0}\to K^{*0}\ell^{+}\ell^{-}$ decays,” hep-ex/2212.09153.
  * (54) [CMS],“Measurement of the B${}^{0}_{\mathrm{S}}$$\to$$\mu^{+}\mu^{-}$ decay properties and search for the B0$\to$$\mu^{+}\mu^{-}$ decay in proton-proton collisions at $\sqrt{s}$ = 13 TeV,” hep-ex/2212.10311.
  * (55) B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, “Dimension-Six Terms in the Standard Model Lagrangian,” JHEP 10, 085 (2010) hep-ph/1008.4884.
  * (56) M. E. Peskin and D. V. Schroeder, “An Introduction to quantum field theory,” Addison-Wesley, 1995, ISBN 978-0-201-50397-5
  * (57) J. R. Ellis, K. A. Olive and C. Savage, “Hadronic Uncertainties in the Elastic Scattering of Supersymmetric Dark Matter,” Phys. Rev. D 77, 065026 (2008) hep-ph/0801.3656.
  * (58) R. Gupta, S. Park, M. Hoferichter, E. Mereghetti, B. Yoon and T. Bhattacharya, “Pion-Nucleon Sigma Term from Lattice QCD,” Phys. Rev. Lett. 127, no.24, 24 (2021) hep-lat/2105.12095.