compat=1.0.0 capbtabboxtable[][]

OU-HET-1160

A Large Muon EDM from Dark Matter

Kim Siang Khaw1,2, Yuichiro Nakai1,2, Ryosuke Sato3,
Yoshihiro Shigekami1,2 and Zhihao Zhang1,2
1Tsung-Dao Lee Institute, Shanghai Jiao Tong University,
520 Shengrong Road, Shanghai 201210, China
2School of Physics and Astronomy, Shanghai Jiao Tong University,
800 Dongchuan Road, Shanghai 200240, China
3Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan

We explore a model of dark matter (DM) that can explain the reported
discrepancy in the muon anomalous magnetic moment and predict a large electric
dipole moment (EDM) of the muon. The model contains a DM fermion and new
scalars whose exclusive interactions with the muon radiatively generate the
observed muon mass. Constraints from DM direct and indirect detection
experiments as well as collider searches are safely evaded. The model
parameter space that gives the observed DM abundance and explains the muon
$g-2$ anomaly leads to the muon EDM of $d_{\mu}\simeq(4\mathchar 45\relax
5)\times 10^{-22}\,e\,{\rm cm}$ that can be probed by the projected PSI muEDM
experiment. Another viable parameter space even achieves
$d_{\mu}=\mathcal{O}(10^{-21})\,e\,{\rm cm}$ reachable by the ongoing Fermilab
Muon $g-2$ experiment and the future J-PARC Muon $g-2$/EDM experiment.

##  1 Introduction

The near-future discovery of the muon electric dipole moment (EDM) is highly
expected by the reported discrepancy in the muon anomalous magnetic moment
$(g-2)_{\mu}$, which may indicate the existence of physics beyond the Standard
Model (SM) at or below the TeV scale [1, 2, 3, 4] (for a review, see ref.
[5]), because the same new physics contribution naturally has the imaginary
part which is relevant to the EDM. The current upper limit on the muon EDM is
$|d_{\mu}|<1.8\times 10^{-19}\,e\,{\rm cm}$ (95% C.L.) [6]. There is also a
study on indirect bounds on the muon EDM by measuring EDMs of heavy atoms and
molecules, which indicates $|d_{\mu}|<2\times 10^{-20}\,e\,{\rm cm}$ [7].
Moreover, the sensitivity to the muon EDM will be improved in the near future:
the ongoing Fermilab Muon $g-2$ experiments [8] and projected J-PARC Muon
$g-2$/EDM experiment [9] will explore the muon EDM at the level of
$10^{-21}\,e\,{\rm cm}$, while the Paul Scherrer Institute (PSI) muEDM
experiment [10, 11, 12] will reach the sensitivity of $6\times
10^{-23}\,e\,{\rm cm}$.

A fermion EDM $d_{f}$ is described by a dimension-five operator
$\mathcal{L}\supset-\frac{i}{2}d_{f}\bar{f}\sigma_{\mu\nu}\gamma^{5}fF^{\mu\nu}$
where $f$ is a Dirac fermion,
$\sigma^{\mu\nu}\equiv\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$ and $F^{\mu\nu}$
is the photon field strength. Since this operator requires a chirality flip
and left- and right-handed fermions carry different charges in the SM, we
actually need a Higgs field insertion which makes the EDM operator effectively
dimension-six. Therefore, a new physics contribution to a fermion EDM scales
as $v_{H}/M^{2}$ where $v_{H}$ and $M$ denote the Higgs vacuum expectation
value (VEV) and a new physics mass scale, respectively.

To estimate the expected size of the muon EDM, we can consider four classes of
new physics that generate the muon EDM as well as the anomalous magnetic
moment 111The similar classification has been presented for the case of the
electron EDM in ref. [13] (see also ref. [14]).:

  * •

Spurion approach. The chirality flip required to generate the EDM operator is
provided by the muon Yukawa coupling $y_{\mu}$ or some coupling proportional
to $y_{\mu}$. When the muon EDM $d_{\mu}$ is generated at the $k$-loop level,
we expect

$\displaystyle d_{\mu}\sim\delta_{\rm
CPV}\left(\frac{\lambda^{2}}{16\pi^{2}}\right)^{k}\frac{m_{\mu}}{M^{2}}\,.$ (
1.1)

Here, $\delta_{\rm CPV}$ and $\lambda$ represent the size of CP-violating
phases and couplings involved in the loop, and $m_{\mu}$ is the muon mass.
Models in this class have been discussed in refs. [15, 16, 17, 18, 19, 20,
21].

  * •

Flavor changing approach. If the muon is converted to the tau lepton by a
lepton flavor violating (LFV) interaction, the chirality flip can be provided
by the tau Yukawa coupling $y_{\tau}$. In this case, we find

$\displaystyle d_{\mu}\sim\delta_{\rm
CPV}\frac{y_{\mu\tau}^{2}}{\lambda^{2}}\left(\frac{\lambda^{2}}{16\pi^{2}}\right)^{k}\frac{m_{\tau}}{M^{2}}\,,$
( 1.2)

with a LFV coupling $y_{\mu\tau}$ and the tau lepton mass $m_{\tau}$. Refs.
[22, 23, 24, 25] have explored models in this class. Note that if the model
has a scalar leptoquark with appropriate charge assignment, the chirality flip
can be picked up from a quark Yukawa coupling, e.g., the top Yukawa coupling.
Moreover, there is an enhancement due to the color factor $N_{C}=3$. Refs.
[26, 27, 28, 29, 30, 31, 32] have explored such a possibility in the context
of the electron EDM and ref. [33] presented a general discussion for the case
of the muon $g-2$ and EDM. In addition, a model with extra vector-like leptons
also has a possibility to predict a large muon EDM [34, 35] due to the
chirality flip on a heavy lepton line.

  * •

Radiative stability approach. New physics that produces the muon EDM also
generates the muon mass by removing the attached photon. When we just assume
that such a contribution to the muon mass does not exceed the correct value,
the size of the muon EDM is expected to be

$\displaystyle d_{\mu}\sim\delta_{\rm CPV}\frac{m_{\mu}}{M^{2}}\,,$ ( 1.3)

because the same loop factor and coupling $\lambda$ are shared by the
generated muon mass and EDM.

  * •

Tuning approach. If the muon mass generated by new physics that produces the
muon EDM exceeds the correct value, a fine-tuning is required. This (unlikely
but logical) possibility allows us to obtain a very large muon EDM which is
bounded by

$\displaystyle d_{\mu}\sim\delta_{\rm
CPV}\lambda\left(\frac{\lambda^{2}}{16\pi^{2}}\right)^{k}\frac{v_{H}}{M^{2}}\lesssim\delta_{\rm
CPV}\frac{4\pi v_{H}}{M^{2}}\,,$ ( 1.4)

for $\lambda\lesssim 4\pi$.

Table 1 shows mass scales of new physics that produce the muon EDM at the
one/two-loop level probed by the projected PSI muEDM experiment [10, 11, 12]
(the ongoing Fermilab Muon $g-2$ experiment [8] and the future J-PARC Muon
$g-2$/EDM experiment [9]). Aside from the tuning approach, the table indicates
that the radiative stability approach generates the largest muon EDM and its
near-future measurements can probe mass scales larger than the TeV scale. The
present paper explores this fascinating possibility for the first time through
the study of a concrete model to realize the radiative stability approach
222In ref. [36], the author commented on the possibility of a large muon EDM
in the context of the radiative stability approach..

| 1-loop | 2-loop
---|---|---
Spurion | 300 GeV (75 GeV) | 16 GeV (4 GeV)
Flavor changing | 580 GeV (140 GeV) | 30 GeV (7 GeV)
Radiative stability | 5900 GeV (1400 GeV)
Tuning |  $1.0\times 10^{6}$ GeV ($2.5\times 10^{5}$ GeV)
Table 1: Mass scales of new physics that produce the muon EDM at the one/two-
loop level probed by the projected PSI muEDM experiment [10, 11, 12] (the
Fermilab Muon $g-2$ and J-PARC Muon $g-2$/EDM experiments [8, 9]) for each
approach presented in the main text. Here, we assume $\lambda\approx 0.65$ of
around the weak coupling constant and $y_{\mu\tau}\approx 0.3$ which is
roughly the maximum value allowed by the measurement of the branching ratio of
$h\to\mu\tau$ (see, e.g., ref. [23]). For the leptoquark model, the mass scale
is enhanced by $(y_{\mu t}/y_{\mu\tau})\sqrt{N_{C}m_{t}/m_{\tau}}\approx
57y_{\mu t}$ where $y_{\mu t}$ is the leptoquark coupling to the muon and the
top quark.

We consider a model of dark matter (DM) that can address the muon $g-2$
anomaly. A DM fermion and new scalars exclusively couple to the muon, which
leads to the radiative generation of the muon mass. The model contains a new
CP-violating phase and produces the muon EDM. We will find that the model
parameter space to give the observed DM abundance and explain the muon $g-2$
anomaly leads to the muon EDM of $d_{\mu}\simeq(4\mathchar 45\relax 5)\times
10^{-22}\,e\,{\rm cm}$ probed by the PSI muEDM experiment. Furthermore, it
will be shown that another viable parameter space even achieves
$d_{\mu}=\mathcal{O}(10^{-21})\,e\,{\rm cm}$ reached by the Fermilab Muon
$g-2$ and J-PARC Muon $g-2$/EDM experiments, which is consistent with the
estimate of Table 1.

The rest of the paper is organized as follows. Section 2 starts with the
description of our DM model and explores the mass spectrum. We then calculate
the radiatively generated muon mass and coupling to the Higgs boson and the
muon EDM as well as the anomalous magnetic moment. They are all induced at the
one-loop level. We also discuss deviations of the muon couplings to the Higgs
and $Z$ bosons from those of the SM. In section 3, phenomenology of DM in our
model is explored. Section 4 summarizes the independent parameters of the
model, and then presents our results to identify the parameter space that
gives the observed DM abundance and explains the muon $g-2$ anomaly and
indicate the size of the muon EDM. In section 5, we give conclusions and
discussions. Loop integrals and full one-loop expressions are summarized in
appendices.

##  2 Model description

Our DM model is based on models proposed in ref. [37], which radiatively
generate the muon mass and explain the muon $g-2$ anomaly. The model contains
a single Dirac fermion $\psi$ and two scalar fields $\phi,\eta$. Charge
assignments for the relevant particles are shown in Table 2.

| $L_{L}^{\mu}$ | $\mu_{R}$ | $H$ | $\psi_{L}$ | $\psi_{R}$ | $\phi$ | $\eta$
---|---|---|---|---|---|---|---
$SU(2)_{L}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$
$Y$ | $-\frac{1}{2}$ | $-1$ | $\frac{1}{2}$ | $Y_{\psi}$ | $Y_{\psi}$ | $Y_{\psi}+\frac{1}{2}$ | $Y_{\psi}+1$
$L_{\mu}$ | $-$ | $-$ | $+$ | $+$ | $+$ | $-$ | $-$
$X$ | $+$ | $+$ | $+$ | $-$ | $-$ | $-$ | $-$
$S_{a}$ | $+$ | $-$ | $+$ | $+$ | $+$ | $+$ | $-$
Table 2: Charge assignments for the relevant particles. $L_{L}^{\mu}$ and
$\mu_{R}$ represent the second generation of the left- and right-handed
leptons and $H$ is the SM Higgs field. The hypercharge of $\psi_{L,R}$ is
taken as $Y_{\psi}=0$ in the present paper. $L_{\mu}$ and $X$ are $Z_{2}$
symmetries associated with the muon number and the exotic particle number,
while $S_{a}$ is a softly broken $Z_{2}$ symmetry to forbid the tree-level
muon Yukawa coupling.

The present paper focuses on the case with $Y_{\psi}=0$ 333The model with
$Y_{\psi}=-1$ has a singlet real scalar $\eta$, and a CP phase appears in the
scalar sector. In this case, however, CP violating effects necessarily involve
the SM Higgs VEV, and therefore, the muon EDM is suppressed when the exotic
particle masses are set to be around TeV.. We introduce two $Z_{2}$ symmetries
$L_{\mu},X$ associated with the muon number and the exotic particle number,
respectively. The former symmetry 444The $L_{\mu}$ symmetry can be enhanced to
a global $U(1)_{L_{\mu}}$ symmetry when $\lambda^{\prime\prime}_{H\phi}$ in
Eq. (2.2) is set to be zero. This value is irrelevant to our current analysis.
Note that even if we do not have $U(1)_{L_{\mu}}$ symmetry, $B-3L_{e}$ number
and $B-3L_{\tau}$ number are conserved (for a review, see, e.g., Ref. [38])
and there is no washout of baryon asymmetry in the early universe. makes it
possible to avoid severe constraints from lepton flavor violating processes,
while the latter one stabilizes the lightest exotic particle which is
identified as DM. In addition, we assume a softly broken $Z_{2}$ symmetry
$S_{a}$ to forbid the tree-level muon Yukawa coupling. The charge assignments
lead to the following terms in the Lagrangian:

$\displaystyle\mathcal{L}$
$\displaystyle\supset\left(-y_{\phi}\overline{L_{L}^{\mu}}\phi^{\dagger}\psi_{R}-y_{\eta}\overline{\psi}_{L}\eta\mu_{R}-m_{D}\overline{\psi}_{L}\psi_{R}-\frac{m_{LL}}{2}\overline{\psi}_{L}\psi^{c}_{L}-\frac{m_{RR}}{2}\overline{\psi^{c}_{R}}\psi_{R}+{\rm
h.c.}\right)-V_{\rm scl}\,,$ ( 2.1) $\displaystyle V_{\rm scl}$
$\displaystyle=\sum_{s=H,\phi,\eta}\left[m_{s}^{2}s^{\dagger}s+\frac{\lambda_{s}}{2}(s^{\dagger}s)^{2}\right]+\lambda_{H\phi}(H^{\dagger}H)(\phi^{\dagger}\phi)+\lambda_{H\eta}(H^{\dagger}H)(\eta^{\dagger}\eta)+\lambda_{\phi\eta}(\phi^{\dagger}\phi)(\eta^{\dagger}\eta)$
$\displaystyle\hskip
11.99998pt+\lambda^{\prime}_{H\phi}(H^{\dagger}\phi)(\phi^{\dagger}H)+\left(aH\eta^{\dagger}\phi+\frac{\lambda^{\prime\prime}_{H\phi}}{2}(H^{\dagger}\phi)^{2}+{\rm
h.c.}\right)\,.$ ( 2.2)

Note that all couplings in $V_{\rm scl}$ and $y_{\phi,\eta}$ can be real and
positive by field redefinitions, while one phase of $m_{D}$, $m_{LL}$, and
$m_{RR}$ cannot be removed. In fact, a combination $m_{LL}m_{RR}/m_{D}^{2}$ is
independent of phase rotations, and we define a physical phase in the model as

$\displaystyle\theta_{\rm
phys}\equiv\frac{1}{2}\arg\left(\frac{m_{LL}m_{RR}}{m_{D}^{2}}\right)=\frac{1}{2}\Bigl{(}\theta_{L}+\theta_{R}-2\theta_{D}\Bigr{)}\,,$
( 2.3)

where $\theta_{D,L,R}$ denote phases of $m_{D}$, $m_{LL}$, and $m_{RR}$,
respectively. Since $\psi$ is singlet under the SM gauge symmetry, we can
define the left- and right-handed Majorana fermions as
$\psi_{L,R}^{M}\equiv\psi_{L,R}+(\psi_{L,R})^{c}$ with
$\psi_{L,R}=P_{L,R}\psi$ and $(\psi_{L,R})^{c}\equiv
i\gamma^{2}(\psi_{L,R})^{*}$. We assume that the exotic scalars $\phi,\eta$ do
not acquire nonzero VEVs. As a result, no mixing between $H$ and $\phi$/$\eta$
is induced, and hence we can parameterize the SM Higgs field $H$ as

$\displaystyle H=\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}\left(v_{H}+h^{0}+iG^{0}\right)\end{pmatrix}\,,$ ( 2.4)

where $v_{H}=246.22$ GeV is the SM Higgs VEV, $G^{+}$ and $G^{0}$ are Nambu-
Goldstone modes, and $h^{0}$ is the SM Higgs boson. Note that a minimization
condition leads to

$\displaystyle m_{H}^{2}$ $\displaystyle=-\frac{1}{2}\lambda_{H}v_{H}^{2}\,.$
( 2.5)

Below, we will present the mass spectrum of exotic particles and calculate the
radiatively generated muon mass and coupling to the Higgs boson and the muon
EDM as well as the anomalous magnetic moment. Deviations of the muon couplings
to the Higgs and $Z$ bosons from those of the SM will be also discussed. Note
that for the neutrino sector, we need a further extension to reproduce the
correct neutrino mixing angles, due to the muon number symmetry. We discuss
some possibilities of the extension in appendix A. We emphasize that such an
extension does not affect our numerical results.

###  2.1 Mass spectrum of exotic particles

From the Lagrangian (2.1), the mass matrix for $\psi_{L}$ and $\psi_{R}$ is

$\displaystyle-\frac{1}{2}\begin{pmatrix}\overline{\psi}_{L}&\overline{\psi^{c}_{R}}\end{pmatrix}\begin{pmatrix}m_{LL}&m_{D}\\\
m_{D}&m_{RR}\end{pmatrix}\begin{pmatrix}\psi_{L}^{c}\\\
\psi_{R}\end{pmatrix}+{\rm
h.c.}=-\frac{1}{2}\begin{pmatrix}\overline{\psi}_{L}&\overline{\psi^{c}_{R}}\end{pmatrix}\mathcal{M}_{\psi}\begin{pmatrix}\psi^{c}_{L}\\\
\psi_{R}\end{pmatrix}+{\rm h.c.}\,,$ ( 2.6)

where $\mathcal{M}_{\psi}$ is a complex symmetric matrix and diagonalized by a
unitary matrix $U_{\psi}$:

$\displaystyle\mathcal{M}_{\psi,{\rm
diag}}=U_{\psi}^{\dagger}\mathcal{M}_{\psi}U_{\psi}^{*}\,,\quad
U_{\psi}=\begin{pmatrix}c_{\alpha}&s_{\alpha}e^{-i\tau}\\\
-s_{\alpha}e^{i\tau}&c_{\alpha}\end{pmatrix}\,.$ ( 2.7)

Here, $c_{\alpha}\equiv\cos\alpha$ with mixing angle $\alpha$ and $\tau$ is
real. In our analysis, we take $m_{LL}$ and $m_{RR}$ to be real and positive,
while $m_{D}$ has a physical phase as $m_{D}=|m_{D}|e^{-i\theta_{\rm phys}}$.
We then obtain mass-squared eigenvalues of
$\mathcal{M}_{\psi}^{\dagger}\mathcal{M}_{\psi}$ as

$\displaystyle m_{\psi_{1}}^{2}$
$\displaystyle=\frac{1}{2}\Bigl{(}m_{LL}^{2}+m_{RR}^{2}+2\left|m_{D}\right|^{2}-\Delta
m_{\psi}^{2}\Bigr{)}\,,$ ( 2.8) $\displaystyle m_{\psi_{2}}^{2}$
$\displaystyle=\frac{1}{2}\Bigl{(}m_{LL}^{2}+m_{RR}^{2}+2\left|m_{D}\right|^{2}+\Delta
m_{\psi}^{2}\Bigr{)}\,,$ ( 2.9)

where $\Delta m_{\psi}^{2}\equiv m_{\psi_{2}}^{2}-m_{\psi_{1}}^{2}$ is given
by

$\displaystyle\Delta
m_{\psi}^{2}=\sqrt{\left(m_{LL}^{2}-m_{RR}^{2}\right)^{2}+4\left|m_{D}\right|^{2}\Bigl{|}m_{LL}e^{-i\theta_{\rm
phys}}+m_{RR}e^{i\theta_{\rm phys}}\Bigr{|}^{2}}\,,$ ( 2.10)

with $\theta_{\rm phys}$ defined in Eq. (2.3). The mixing angle $\alpha$ and
phase $\tau$ in $U_{\psi}$ are obtained as

$\displaystyle\sin 2\alpha$ $\displaystyle=\frac{2|m_{D}|}{\Delta
m_{\psi}^{2}}\Bigl{|}m_{LL}e^{-i\theta_{\rm phys}}+m_{RR}e^{i\theta_{\rm
phys}}\Bigr{|}\,,$ ( 2.11) $\displaystyle\tan\tau$
$\displaystyle=-\frac{m_{LL}-m_{RR}}{m_{LL}+m_{RR}}\tan\theta_{\rm phys}\,.$ (
2.12)

Due to the mass hierarchy, $m_{\psi_{2}}^{2}>m_{\psi_{1}}^{2}$, we can focus
on $0\leq\alpha\leq\pi/2$. Note that physical predictions are unchanged for
$\theta_{\rm phys}\to\theta_{\rm phys}+\pi$, and we focus on the range of
$-\pi/2<\theta_{\rm phys}\leq\pi/2$ in our analysis. $\psi_{L,R}$ can be
described in terms of mass eigenstates $\psi_{1,2}$ as

$\displaystyle\psi_{L}=\psi_{1}^{c}c_{\alpha}+\psi_{2}^{c}s_{\alpha}e^{-i\tau}\,,\quad\psi_{R}=-\psi_{1}s_{\alpha}e^{-i\tau}+\psi_{2}c_{\alpha}\,,$
( 2.13)

and the mass terms in Eq. (2.6) become

$\displaystyle-\frac{1}{2}\Bigl{[}m_{\psi_{1}}\overline{\psi_{1}^{M}}\psi_{1}^{M}+m_{\psi_{2}}\overline{\psi_{2}^{M}}\psi_{2}^{M}\Bigr{]}\,,$
( 2.14)

where $\psi_{1,2}^{M}\equiv\psi_{1,2}+\psi_{1,2}^{c}$ are Majorana fermions.

In order to analyze the mass spectrum for exotic scalar fields, we
parameterize them as

$\displaystyle\phi=\begin{pmatrix}\phi^{+}\\\
\frac{1}{\sqrt{2}}\left(\sigma_{\phi}+ia_{\phi}\right)\end{pmatrix}\,,\qquad\eta=\eta^{+}\,.$
( 2.15)

From Eq. (2.2), the mass-squared matrices for charged and neutral scalars (in
the basis of $(\phi^{+},\eta^{+})$ and $(\sigma_{\phi},a_{\phi})$,
respectively) are given by

$\displaystyle\mathcal{M}_{\pm}^{2}$
$\displaystyle=\begin{pmatrix}{\displaystyle
M_{\phi}^{2}}&{\displaystyle\frac{av_{H}}{\sqrt{2}}}\\\\[10.76385pt]
{\displaystyle\frac{av_{H}}{\sqrt{2}}}&{\displaystyle
M_{\eta}^{2}}\end{pmatrix}\,,$ ( 2.16) $\displaystyle\mathcal{M}_{0}^{2}$
$\displaystyle=\begin{pmatrix}{\displaystyle
M_{\phi}^{2}+\frac{\lambda^{+}_{H\phi}}{2}v_{H}^{2}}&0\\\\[8.61108pt]
0&{\displaystyle
M_{\phi}^{2}+\frac{\lambda^{-}_{H\phi}}{2}v_{H}^{2}}\end{pmatrix}\equiv\begin{pmatrix}m_{\sigma_{\phi}}^{2}&0\\\
0&m_{a_{\phi}}^{2}\end{pmatrix}\,,$ ( 2.17)

where $M_{\phi,\eta}^{2}\equiv
m_{\phi,\eta}^{2}+\frac{\lambda_{H\phi,H\eta}}{2}v_{H}^{2}$ and
$\lambda^{\pm}_{H\phi}\equiv\lambda^{\prime}_{H\phi}\pm\lambda^{\prime\prime}_{H\phi}$.
Note that since all quartic couplings are positive, $m_{a_{\phi}}^{2}$ is
always smaller than $m_{\sigma_{\phi}}^{2}$. Diagonalzation of
$\mathcal{M}_{\pm}^{2}$ can be done by an orthogonal matrix as

$\displaystyle\mathcal{M}_{\pm,{\rm
diag}}^{2}=U_{s}^{T}\mathcal{M}_{\pm}^{2}U_{s},\quad
U_{s}\equiv\begin{pmatrix}c_{\theta}&s_{\theta}\\\
-s_{\theta}&c_{\theta}\end{pmatrix}\,.$ ( 2.18)

The mass-squared eigenvalues and the mixing angle $\theta$ are

$\displaystyle m_{\varphi_{1}^{+}}^{2}$
$\displaystyle=\frac{1}{2}\left[M_{\phi}^{2}+M_{\eta}^{2}-\sqrt{(M_{\phi}^{2}-M_{\eta}^{2})^{2}+2a^{2}v_{H}^{2}}\right]\,,$
( 2.19) $\displaystyle m_{\varphi_{2}^{+}}^{2}$
$\displaystyle=\frac{1}{2}\left[M_{\phi}^{2}+M_{\eta}^{2}+\sqrt{(M_{\phi}^{2}-M_{\eta}^{2})^{2}+2a^{2}v_{H}^{2}}\right]\,,$
( 2.20) $\displaystyle\sin 2\theta$
$\displaystyle=\frac{\sqrt{2}av_{H}}{m_{\varphi_{2}^{+}}^{2}-m_{\varphi_{1}^{+}}^{2}}\,.$
( 2.21)

Then, $\phi^{\pm}$ and $\eta^{\pm}$ can be described in terms of mass
eigenstates $\varphi_{1,2}^{\pm}$ as

$\displaystyle\phi^{\pm}=\varphi_{1}^{\pm}c_{\theta}+\varphi_{2}^{\pm}s_{\theta}\,,\qquad\eta^{\pm}=-\varphi_{1}^{\pm}s_{\theta}+\varphi_{2}^{\pm}c_{\theta}\,.$
( 2.22)

Since the mass parameter $a$ can be set to be real and positive and
$m_{\varphi_{2}^{+}}^{2}-m_{\varphi_{1}^{+}}^{2}>0$, we can focus on
$0\leq\theta\leq\pi/2$.

###  2.2 Radiative mass and coupling of the muon

The mass and Yukawa coupling of the muon are induced by one-loop corrections.
When we move to the mass basis for exotic particles according to Eqs. (2.13)
and (2.22), the relevant terms are written as

$\displaystyle\mathcal{L}\supset\Bigl{(}-y_{L}^{ia}\bar{\mu}_{L}\varphi_{i}^{-}\psi_{a}-y_{R}^{ia}\overline{\psi_{a}^{c}}\varphi_{i}^{+}\mu_{R}+{\rm
h.c.}\Bigr{)}-\frac{A_{ij}}{\sqrt{2}}h^{0}\varphi_{i}^{-}\varphi_{j}^{+}\,,$ (
2.23)

where the explicit forms of $y_{L,R}^{ia}$ and $A_{ij}$ are summarized in
Table 3.

$(i,a)$ | $y_{L}^{ia}$ | $y_{R}^{ia}$ | $(i,j)$ | $A_{ij}$
---|---|---|---|---
$(1,1)$ | $-y_{\phi}c_{\theta}s_{\alpha}e^{-i\tau}$ | $-y_{\eta}s_{\theta}c_{\alpha}$ | $(1,1)$ | $-as_{2\theta}+\sqrt{2}v_{H}\left(\lambda_{H\phi}c_{\theta}^{2}+\lambda_{H\eta}s_{\theta}^{2}\right)$
$(1,2)$ | $y_{\phi}c_{\theta}c_{\alpha}$ | $-y_{\eta}s_{\theta}s_{\alpha}e^{i\tau}$ | $(1,2)$ | $ac_{2\theta}+\sqrt{2}v_{H}\left(\lambda_{H\phi}-\lambda_{H\eta}\right)s_{\theta}c_{\theta}$
$(2,1)$ | $-y_{\phi}s_{\theta}s_{\alpha}e^{-i\tau}$ | $y_{\eta}c_{\theta}c_{\alpha}$ | $(2,1)$ | $ac_{2\theta}+\sqrt{2}v_{H}\left(\lambda_{H\phi}-\lambda_{H\eta}\right)s_{\theta}c_{\theta}$
$(2,2)$ | $y_{\phi}s_{\theta}c_{\alpha}$ | $y_{\eta}c_{\theta}s_{\alpha}e^{i\tau}$ | $(2,2)$ | $as_{2\theta}+\sqrt{2}v_{H}\left(\lambda_{H\phi}s_{\theta}^{2}+\lambda_{H\eta}c_{\theta}^{2}\right)$
Table 3: Yukawa couplings for the muon and exotic particles and scalar
trilinear couplings in Eq. (2.23).

These couplings lead to the radiative mass and effective Yukawa coupling of
the muon at the one-loop level, through diagrams in Fig. 1:

$\displaystyle\mathcal{L}_{\rm eff}$ $\displaystyle\supset-m_{\mu}^{\rm
rad}\bar{\mu}_{L}\mu_{R}-\frac{y_{\mu}^{\rm
eff}}{\sqrt{2}}\bar{\mu}_{L}\mu_{R}h^{0}+{\rm h.c.}\,,$ ( 2.24) $\displaystyle
m_{\mu}^{\rm rad}$
$\displaystyle=\sum_{i,a}\frac{y_{L}^{ia}y_{R}^{ia}}{16\pi^{2}}m_{\psi_{a}}\left(\frac{m_{\psi_{a}}^{2}B_{0}(0,m_{\psi_{a}}^{2},m_{\psi_{a}}^{2})-m_{\varphi_{i}^{+}}^{2}B_{0}(0,m_{\varphi_{i}^{+}}^{2},m_{\varphi_{i}^{+}}^{2})}{m_{\varphi_{i}^{+}}^{2}-m_{\psi_{a}}^{2}}-1\right)$
$\displaystyle=\frac{y_{\phi}y_{\eta}}{16\pi^{2}}\frac{s_{2\theta}s_{2\alpha}}{4}\mathcal{F}(x_{1,1},x_{1,2},x_{2,1},x_{2,2})\,,$
( 2.25) $\displaystyle y_{\mu}^{\rm eff}(p_{h^{0}}^{2})$
$\displaystyle=-\sum_{i,j,a}\frac{y_{L}^{ia}y_{R}^{ja}A_{ij}}{16\pi^{2}}m_{\psi_{a}}C_{0}(m_{\mu}^{2},m_{\mu}^{2},p_{h^{0}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{j}^{+}}^{2})\,,$
( 2.26)

where $p_{h^{0}}$ is the four-momentum of the SM Higgs boson,
$m_{\psi_{a}}\equiv\sqrt{m_{\psi_{a}}^{2}}$, and $B_{0}(0,m^{2},m^{2})$ and
$C_{0}(m_{\mu}^{2},m_{\mu}^{2},p_{h^{0}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{j}^{+}}^{2})$
denote loop integrals for the self-energy and triangle type diagrams,
respectively, whose explicit forms are summarized in Appendix B, and
$\mathcal{F}(x_{1,1},x_{1,2},x_{2,1},x_{2,2})$ is defined as

$\displaystyle\mathcal{F}(x_{1,1},x_{1,2},x_{2,1},x_{2,2})$
$\displaystyle\equiv-
m_{\psi_{1}}e^{-i\tau}\left(\frac{x_{1,1}}{1-x_{1,1}}\mathop{\rm
ln}x_{1,1}-\frac{x_{2,1}}{1-x_{2,1}}\mathop{\rm ln}x_{2,1}\right)$
$\displaystyle\hskip
40.00006pt+m_{\psi_{2}}e^{i\tau}\left(\frac{x_{1,2}}{1-x_{1,2}}\mathop{\rm
ln}x_{1,2}-\frac{x_{2,2}}{1-x_{2,2}}\mathop{\rm ln}x_{2,2}\right)\,,$ ( 2.27)

with $x_{i,a}\equiv m_{\varphi_{i}^{+}}^{2}/m_{\psi_{a}}^{2}$. Eq. (2.26)
neglects sub-dominant contributions with a chirality flip on the muon line. We
show the full form for $y_{\mu}^{\rm eff}$ at the one-loop order in Appendix
C. To numerically evaluate the loop functions with a non-zero $p_{h^{0}}^{2}$,
we use LoopTools [39]. Note that due to the radiatively generated muon mass,
there is no standard relation between $m_{\mu}^{\rm rad}$ and $y_{\mu}^{\rm
eff}$, namely, $m_{\mu}^{\rm rad}\neq y_{\mu}^{\rm eff}v_{H}/\sqrt{2}$. Hence,
we need to check if the model satisfies a constraint from the measurement of
$h\to\mu^{+}\mu^{-}$. We discuss this constraint in Sec. 2.3.

$\mu_{R}$$\mu_{L}$$h^{0}$$\psi_{a}$$\varphi_{j}^{+}$$\varphi_{i}^{-}$$\mu_{R}$$\mu_{L}$$\gamma$$\psi_{a}$$\varphi_{i}^{+}$$\varphi_{i}^{-}$
Figure 1: Feynman diagrams for (a) the effective Yukawa coupling of the muon
and (b) dipole operators. The diagram for the radiative mass of the muon is
obtained by eliminating the external photon line from the diagram (b).

Since $m_{\mu}^{\rm rad}$ in Eq. (2.25) generally has a phase due to complex
couplings $y_{L,R}^{ia}$, we need to remove it by a chiral rotation of the
muon field as

$\displaystyle\mu\to e^{-i\theta_{\mu}\gamma_{5}/2}\mu\,,$ ( 2.28)

where $\theta_{\mu}$ is defined as $m_{\mu}^{\rm
rad}=m_{\mu}e^{i\theta_{\mu}}$. Here, $m_{\mu}$ is understood as the observed
muon mass and a real value, and $\theta_{\mu}$ can be obtained as
$\theta_{\mu}=\arg\bigl{[}\mathcal{F}(x_{1,1},x_{1,2},x_{2,1},x_{2,2}))\bigr{]}$.
This rotation affects dipole operators,

$\displaystyle\mathcal{L}_{\rm
dipole}=-\frac{e}{2}C_{T}(q^{2})\left(\bar{\mu}\sigma^{\alpha\beta}\mu\right)F_{\alpha\beta}-\frac{e}{2}C_{T^{\prime}}(q^{2})\left(\bar{\mu}i\sigma^{\alpha\beta}\gamma_{5}\mu\right)F_{\alpha\beta}\,,$
( 2.29)

where $q$ is the four-momentum for the photon. If there was no chiral
rotation, $C_{T}(0)$ and $C_{T^{\prime}}(0)$ would be the muon $g-2$ $a_{\mu}$
and the muon EDM $d_{\mu}$, respectively. After performing the chiral rotation
of Eq. (2.28), we can obtain the correct forms of $a_{\mu}$ and $d_{\mu}$ in
our model as

$\displaystyle\mathcal{L}_{\rm dipole}$
$\displaystyle=-\frac{e}{4m_{\mu}}a_{\mu}\left(\bar{\mu}\sigma^{\alpha\beta}\mu\right)F_{\alpha\beta}-\frac{i}{2}d_{\mu}\left(\bar{\mu}\sigma^{\alpha\beta}\gamma_{5}\mu\right)F_{\alpha\beta}\,,$
( 2.30) $\displaystyle a_{\mu}$
$\displaystyle=2m_{\mu}\left(C_{T}(0)\cos\theta_{\mu}+C_{T^{\prime}}(0)\sin\theta_{\mu}\right)\,,$
( 2.31) $\displaystyle d_{\mu}$
$\displaystyle=e\left(C_{T^{\prime}}(0)\cos\theta_{\mu}-C_{T}(0)\sin\theta_{\mu}\right)\,.$
( 2.32)

The leading contributions to $C_{T}(0)$ and $C_{T^{\prime}}(0)$ can be
estimated from the diagram (b) in Fig. 1 as

$\displaystyle C_{T}(0)$ $\displaystyle=\sum_{i,a}\frac{{\rm
Re}[y_{L}^{ia}y_{R}^{ia}]}{16\pi^{2}}m_{\psi_{a}}\Bigl{[}C_{0}(m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})+2C_{1}(m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\Bigr{]}\,,$
( 2.33) $\displaystyle C_{T^{\prime}}(0)$ $\displaystyle=\sum_{i,a}\frac{{\rm
Im}[y_{L}^{ia}y_{R}^{ia}]}{16\pi^{2}}m_{\psi_{a}}\Bigl{[}C_{0}(m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})+2C_{1}(m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\Bigr{]}\,,$
( 2.34)

where

$\displaystyle C_{0}(m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})$
$\displaystyle\equiv
C_{0}(m_{\mu}^{2},m_{\mu}^{2},0,m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\approx
C_{0}(0,0,0,m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\,,$
( 2.35) $\displaystyle C_{1}(m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})$
$\displaystyle\equiv
C_{1}(m_{\mu}^{2},0,m_{\mu}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\varphi_{i}^{+}}^{2})\approx
C_{1}(0,0,0,m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\varphi_{i}^{+}}^{2})\,,$
( 2.36)

are loop integrals for the triangle type diagram, and approximations in the
right hand sides are valid when $m_{\mu}^{2}\ll m_{\varphi_{i}^{+}}^{2}$,
$m_{\psi_{a}}^{2}$. In this case, we can obtain the following analytical forms
of $C_{0}$ and $C_{1}$:

$\displaystyle
C_{0}(0,0,0,m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})$
$\displaystyle=\frac{1}{m_{\psi_{a}}^{2}}\left[\frac{1}{1-x_{i,a}}+\frac{1}{(1-x_{i,a})^{2}}\mathop{\rm
ln}x_{i,a}\right]\,,$ ( 2.37) $\displaystyle
C_{1}(0,0,0,m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\varphi_{i}^{+}}^{2})$
$\displaystyle=-\frac{1}{m_{\psi_{a}}^{2}}\left[\frac{3-x_{i,a}}{4(1-x_{i,a})^{2}}+\frac{1}{2(1-x_{i,a})^{3}}\mathop{\rm
ln}x_{i,a}\right]\,.$ ( 2.38)

Here, $x_{i,a}$ is defined below Eq. (2.27). Since the leading contributions
to $C_{T}(0)$ and $C_{T^{\prime}}(0)$ are the same except for the overall
couplings, ${\rm Re}[y_{L}^{ia}y_{R}^{ia}]$ and ${\rm
Im}[y_{L}^{ia}y_{R}^{ia}]$, we can expect a sufficiently large $d_{\mu}$ to be
probed in near-future experiments when the muon $g-2$ is predicted to be
$\mathcal{O}(10^{-9})$. That is, when
$C_{T}(0)\cos\theta_{\mu}+C_{T^{\prime}}(0)\sin\theta_{\mu}\sim
C_{T^{\prime}}(0)\cos\theta_{\mu}-C_{T}(0)\sin\theta_{\mu}$ is satisfied, we
find

$\displaystyle d_{\mu}\sim\frac{e}{2m_{\mu}}\times a_{\mu}\simeq 2.34\times
10^{-22}\,e\,{\rm cm}\,,$ ( 2.39)

with $a_{\mu}\simeq 2.51\times 10^{-9}$. By using couplings in Table. 3 and
Eqs. (2.25), (2.37) and (2.38), $C_{T}(0)$ and $C_{T^{\prime}}(0)$ can be
rewritten as

$\displaystyle C_{T}(0)$
$\displaystyle=\frac{m_{\mu}{\cos\tau}}{\bigl{|}\mathcal{F}(x_{1,1},x_{1,2},x_{2,1},x_{2,2})\bigr{|}}\sum_{i,a=1}^{2}\frac{{(-1)^{i+a}}}{m_{\psi_{a}}}\frac{x_{i,a}^{2}-1-2x_{i,a}\mathop{\rm
ln}x_{i,a}}{2(1-x_{i,a})^{3}}\,,$ ( 2.40) $\displaystyle C_{T^{\prime}}(0)$
$\displaystyle=\frac{m_{\mu}{\sin\tau}}{\bigl{|}\mathcal{F}(x_{1,1},x_{1,2},x_{2,1},x_{2,2})\bigr{|}}\sum_{i,a=1}^{2}\frac{{(-1)^{i}}}{m_{\psi_{a}}}\frac{x_{i,a}^{2}-1-2x_{i,a}\mathop{\rm
ln}x_{i,a}}{2(1-x_{i,a})^{3}}\,.$ ( 2.41)

As $\bigl{|}\mathcal{F}\bigr{|}$ is proportional to $m_{\psi_{1,2}}$, their
scalings are consistent with the rough estimation given in Eq. (1.3). It is
notable that when we change $\theta_{\rm phys}\to-\theta_{\rm phys}$, signs of
$\sin\tau$ and $\sin\theta_{\mu}$ are flipped, the former of which leads to
$C_{T^{\prime}}(0)\to-C_{T^{\prime}}(0)$ through Eq. (2.41), and hence, this
change results in $d_{\mu}\to-d_{\mu}$ with $a_{\mu}$ unchanged. This fact
tells us that it is enough to focus on the range $0<\theta_{\rm phys}<\pi/2$,
because we are only interested in the prediction of $|d_{\mu}|$ here.
Furthermore, $\theta_{\rm phys}=0$ corresponds to a CP conserving limit which
gives $|d_{\mu}|=0$, while $\theta_{\rm phys}=\pi/2$ leads to
$\tau\approx\pi/2$ unless $m_{LL}=m_{RR}$ (see Eq. (2.12)), predicting
$\cos\theta_{\mu}\approx 0$, and hence, $|d_{\mu}|\propto
C_{T^{\prime}}(0)\cos\theta_{\mu}-C_{T}(0)\sin\theta_{\mu}\approx 0$.
Hereafter, we denote $d_{\mu}$ as its absolute value $|d_{\mu}|$ in our
analysis.

###  2.3 Muon coupling constraints

In our model, the muon Yukawa coupling to the Higgs boson is generated at the
one-loop level and does not follow the standard relation, $m_{\mu}^{\rm
rad}\neq y_{\mu}^{\rm eff}v_{H}/\sqrt{2}$. The ATLAS [40] and CMS [41]
experiments have searched for the Higgs boson decay $h\to\mu^{+}\mu^{-}$,
which lead to constraints on the $h$-$\mu$-$\mu$ coupling as

$\displaystyle\left|\kappa_{\mu}\right|$ $\displaystyle<1.47\quad({\rm
ATLAS})\,,$ ( 2.42) $\displaystyle 0.61<\left|\kappa_{\mu}\right|$
$\displaystyle<1.44\quad({\rm CMS})\,,$ ( 2.43)

where we use BR$(h\to\mu^{+}\mu^{-})_{\rm SM}\simeq 2.16\times 10^{-4}$ for
$m_{h}=125.25$ GeV [42], and $\kappa_{\mu}$ is defined by comparing the decay
width of $h\to\mu^{+}\mu^{-}$ to that of the SM,

$\displaystyle\Gamma_{h\to\mu^{+}\mu^{-}}^{\rm
SM}=\frac{m_{h}}{8\pi}\left(\frac{m_{\mu}}{v_{H}}\right)^{2}\left(1-\frac{4m_{\mu}^{2}}{m_{h}^{2}}\right)^{3/2}\,.$
( 2.44)

In our model, the width of $h\to\mu^{+}\mu^{-}$ is estimated as

$\displaystyle\Gamma_{h\to\mu^{+}\mu^{-}}=\frac{m_{h}}{16\pi}\sqrt{1-\frac{4m_{\mu}^{2}}{m_{h}^{2}}}\left[\left(1-\frac{4m_{\mu}^{2}}{m_{h}^{2}}\right)\left({\rm
Re}\>y_{\mu}^{\rm eff}\right)^{2}+\left({\rm Im}\>y_{\mu}^{\rm
eff}\right)^{2}\right]\,.$ ( 2.45)

Then, we find

$\displaystyle\left|\kappa_{\mu}\right|$
$\displaystyle=\frac{1}{\sqrt{2}}\frac{v_{H}}{m_{\mu}}\sqrt{\left({\rm
Re}\>y_{\mu}^{\rm
eff}\right)^{2}+\left(1-\frac{4m_{\mu}^{2}}{m_{h}^{2}}\right)^{-1}\left({\rm
Im}\>y_{\mu}^{\rm
eff}\right)^{2}}\approx\frac{1}{\sqrt{2}}\frac{v_{H}}{m_{\mu}}\left|y_{\mu}^{\rm
eff}\right|\,.$ ( 2.46)

Here, we have used $4m_{\mu}^{2}\ll m_{h}^{2}$.

Since exotic particles exclusively couple to the muon, the ratio between the
$Z\to e^{+}e^{-}$ and $Z\to\mu^{+}\mu^{-}$ decay widths may constrain our
parameter space. The current experimental status for this ratio is [43]

$\displaystyle\frac{\Gamma(Z\to\mu^{+}\mu^{-})}{\Gamma(Z\to
e^{+}e^{-})}=1.0009\pm 0.0028\,.$ ( 2.47)

The muon couplings to the $Z$ boson can be parameterized as

$\displaystyle\mathcal{L}_{Z}\supset\frac{g}{\cos\theta_{W}}\bar{\mu}\gamma^{\alpha}\Bigl{[}(g_{L}^{\mu}+\delta
g_{L}^{\mu})P_{L}+(g_{R}^{\mu}+\delta g_{R}^{\mu})P_{R}\Bigr{]}\mu
Z_{\alpha}\,,$ ( 2.48)

where $g$ denotes the $SU(2)_{L}$ gauge coupling, $\theta_{W}$ is the weak
mixing angle, and $g_{L}^{\mu}=-\frac{1}{2}+\sin^{2}\theta_{W}$, and
$g_{R}^{\mu}=\sin^{2}\theta_{W}$ are the muon couplings to the $Z$ boson in
the SM. In our model, new physics contributions $\delta g_{L,R}$ are induced
by the diagram (b) in Fig. 1 with replacing the photon to the $Z$ boson, and
their expressions are found in ref. [44]. The ratio in Eq. (2.47) is then
estimated as

$\displaystyle\frac{\Gamma(Z\to\mu^{+}\mu^{-})}{\Gamma(Z\to e^{+}e^{-})}\simeq
1+\frac{2g_{L}^{e}{\rm Re}\left(\delta g_{L}^{\mu}\right)+2g_{R}^{e}{\rm
Re}\left(\delta g_{R}^{\mu}\right)}{(g_{L}^{e})^{2}+(g_{R}^{e})^{2}}\equiv
1+\delta_{\mu\mu}\,,$ ( 2.49)

where $g_{L,R}^{e}=g_{L,R}^{\mu}$ are the electron couplings to the $Z$ boson
in the SM, and we assume that new physics contributions are smaller than those
of the SM, $\delta g_{L,R}^{\mu}\ll g_{L,R}^{\mu}$. Then, Eq. (2.47) indicates
that $|\delta_{\mu\mu}|$ must be less than $\mathcal{O}(10^{-3})$.

##  3 Dark matter

The candidate of DM in our model is the lightest Majorana fermion
$\psi_{1}^{M}$ or the lightest neutral scalar $a_{\phi}$, depending on their
masses. In the present paper, we focus on the case that $\psi_{1}^{M}$ is the
lightest exotic particle, and hence, gives the DM candidate. Hereafter, we
denote $\psi_{1}^{M}$ as $\psi_{1}$ for simplicity. For the case that
$a_{\phi}$ is the DM candidate, there is no direct correlation with the muon
EDM, because the mass $m_{a_{\phi}}$ does not contribute to the muon EDM at
the one-loop level.

The main annihilation mode of the DM fermion is
$\psi_{1}\psi_{1}\to\mu\bar{\mu}$ through the $t$-channel exchange of
$\varphi_{i}^{\pm}$, as shown in the left of Fig. 2.

$\psi_{1}$$\mu$$\psi_{1}$$\bar{\mu}$$\varphi_{i}^{\pm}$$\psi_{1}$$\mu$$\varphi_{1}^{+}$$V$$\varphi_{i}^{-}$$\varphi_{1}^{-}$$\mu$$\varphi_{1}^{+}$$\bar{\mu}$$\psi_{a}$
Figure 2: Feynman diagrams for DM self-annihilation (left), DM-charged scalar
coannihilation (center) and charged scalar self-annihilation (right). $V$
denotes the $\gamma/Z$ boson. Other diagrams can be obtained by changing final
states with an appropriate mediator.

In the expansion of the thermally averaged cross section by the DM velocity
$v$, $\langle\sigma
v\rangle_{\mu\bar{\mu}}=a_{\mu\bar{\mu}}+b_{\mu\bar{\mu}}v^{2}+\mathcal{O}(v^{4})$,
$s$-wave and $p$-wave contributions are given by

$\displaystyle a_{\mu\bar{\mu}}$ $\displaystyle=\frac{1}{16\pi
m_{\psi_{1}}^{2}}\left[\left|\sum_{i=1,2}\frac{y_{L}^{i1}y_{R}^{i1}}{1+x_{i,1}}\right|^{2}+\mathcal{O}\left(\frac{m_{\mu}^{2}}{m_{\psi_{1}}^{2}}\right)\right]\,,$
( 3.1) $\displaystyle b_{\mu\bar{\mu}}$ $\displaystyle=\frac{1}{48\pi
m_{\psi_{1}}^{2}}\left(\sum_{i=1,2}\frac{|y_{L}^{i1}|^{2}}{(1+x_{i,1})^{2}}\sqrt{1+x_{i,1}^{2}}\right)^{2}+(L\to
R)\,,$ ( 3.2)

where the second term in Eq. (3.1) is suppressed by $m_{\mu}/m_{\psi_{1}}$.
Thus, the $s$-wave contribution dominates the total DM annihilation cross
section in our focused parameter space. Note that for the annihilation mode
$\psi_{1}\psi_{1}\to\nu_{\mu}\bar{\nu}_{\mu}$, the $s$-wave contribution is
suppressed by a tiny neutrino mass, because there is no right-handed coupling
$y_{R}^{i1}$ for the neutrino. The other annihilation cross sections, such as
$\psi_{1}\psi_{1}\to\gamma\gamma$ and $\psi_{1}\psi_{1}\to\mu\bar{\mu}\gamma$,
are several orders of magnitude smaller than that of
$\psi_{1}\psi_{1}\to\mu\bar{\mu}$.

In the thermal freeze-out scenario, the number density of DM is calculated by
the Boltzmann equation,

$\displaystyle\frac{dn_{\psi_{1}}(t)}{dt}+3H(t)n_{\psi_{1}}(t)=-\langle\sigma
v\rangle_{\rm eff}\Bigl{[}n_{\psi_{1}}(t)^{2}-n_{\psi_{1}}^{\rm
eq}(t)^{2}\Bigr{]}\,.$ ( 3.3)

Here, $H(t)$ denotes the Hubble rate, and $n_{\psi_{1}}$ is the number density
of $\psi_{1}$, while $n_{\psi_{1}}^{\rm eq}$ is that in equilibrium. The
effective annihilation cross section $\langle\sigma v\rangle_{\rm eff}$ is
estimated by summing all possible annihilation modes, i.e.,
$\psi_{1}\psi_{1}\to\ell\bar{\ell},VV^{\prime},\ell\bar{\ell}V$
($\ell=\mu,\nu;V,V^{\prime}=\gamma,Z,W$). However, when the DM and charged
scalar masses are almost degenerate, coannihilation processes should be taken
into account for solving the Boltzmann equation. In this case, we have [45]

$\displaystyle(\sigma v)_{\rm eff}$
$\displaystyle=\frac{1}{(g_{\psi_{1}}+\bar{g}_{\varphi_{1}^{+}})^{2}}\Bigl{[}g_{\psi_{1}}^{2}(\sigma
v)_{\psi_{1}\psi_{1}}+g_{\psi_{1}}\bar{g}_{\varphi_{1}^{+}}(\sigma
v)_{\psi_{1}\varphi_{1}^{+}}+\bar{g}_{\varphi_{1}^{+}}^{2}(\sigma
v)_{\varphi_{1}^{+}\varphi_{1}^{-}}\Bigr{]}\,,$ ( 3.4)
$\displaystyle\text{with
}\bar{g}_{\varphi_{1}^{+}}=g_{\varphi_{1}^{+}}\left(\frac{m_{\varphi_{1}^{+}}}{m_{\psi_{1}}}\right)^{3/2}e^{-(m_{\varphi_{1}^{+}}-m_{\psi_{1}})/T}\,,$
( 3.5)

where $g_{\psi_{1}}=2$ and $g_{\varphi_{1}^{+}}=2$ are internal degrees of
freedom for $\psi_{1}$ and $\varphi_{1}^{+}$, respectively, $T$ is the
temperature, and $(\sigma v)_{XY}$ denotes the (co)annihilation cross section
whose initial state is $XY$. The corresponding diagrams are shown in Fig. 2.
The second term in Eq. (3.4) is suppressed by the exponential factor in Eq.
(3.5) and the third term is more suppressed due to the squared exponential
factor when $m_{\varphi_{1}^{+}}\gg m_{\psi_{1}}$. As $m_{\varphi_{1}^{+}}$
decreases and is close to $m_{\psi_{1}}$, the second term gives a non-
negligible contribution to $(\sigma v)_{\rm eff}$ [46, 47]. The resultant DM
relic density is given by

$\displaystyle\Omega
h^{2}=\frac{m_{\psi_{1}}n_{\psi_{1}}(t_{0})}{\rho_{c}}h^{2}\,,$ ( 3.6)

where $\rho_{c}$ is the critical density of the Universe and
$n_{\psi_{1}}(t_{0})$ is the today’s number density of $\psi_{1}$ obtained by
solving the Boltzmann equation (3.3). To calculate the DM relic density
$\Omega h^{2}$ including appropriate coannihilation processes, we use
micrOMEGAs_5.2.13 [48, 49].

Although our DM particle $\psi_{1}$ does not couple to the SM quarks and
gluons, the DM-nucleon scattering is induced by contact and non-contact type
interactions. In our model, relevant interactions for the scattering are

$\displaystyle\mathcal{L}_{\rm eff}\supset
a_{\psi_{1}}\overline{\psi_{1}}\gamma^{\mu}\gamma^{5}\psi_{1}\partial^{\nu}F_{\mu\nu}+C_{S,p}\overline{\psi_{1}}\psi_{1}\bar{p}p+C_{S,n}\overline{\psi_{1}}\psi_{1}\bar{n}n\,.$
( 3.7)

Here, $p$ and $n$ represent the proton and the neutron, and the effective
coefficients are estimated as

$\displaystyle a_{\psi_{1}}$
$\displaystyle=-\frac{e}{16\pi^{2}m_{\psi_{1}}^{2}}\sum_{i=1,2}\Bigl{(}|y_{L}^{i1}|^{2}+|y_{R}^{i1}|^{2}\Bigr{)}\hat{a}_{\psi_{1}}(x_{i,1},\epsilon_{\mu})\,,$
( 3.8) $\displaystyle C_{S,N}$
$\displaystyle=-m_{N}\sum_{q}C_{S,q}\frac{f_{Tq}^{(N)}}{m_{q}}=-\frac{y_{\psi_{1}}^{\rm
eff}m_{N}}{\sqrt{2}m_{h}^{2}v_{H}}\sum_{q}f_{Tq}^{(N)}\qquad(N=p,n)\,,$ ( 3.9)

where $x_{i,1}$ is defined below Eq. (2.27), $\epsilon_{\mu}\equiv
m_{\mu}^{2}/m_{\psi_{1}}^{2}$, $\hat{a}_{\psi_{1}}(x,y)$ is the loop function
for the anapole operator, which is given by

$\displaystyle\hat{a}_{\psi_{1}}(x,y)=\frac{1}{12}\left[\frac{3}{2}\mathop{\rm
ln}\left(\frac{y}{x}\right)+\frac{3x-3y+1}{\Delta(x,y)^{1/2}}\tanh^{-1}\left(\frac{\Delta(x,y)^{1/2}}{x+y-1}\right)\right]\,,$
( 3.10)

with $\Delta(x,y)=x^{2}+(y-1)^{2}-2x(y+1)$, $C_{S,q}$ denotes the effective
coupling of an operator $\overline{\psi_{1}}\psi_{1}\bar{q}q$ with the SM
quark $q$, and $f_{Tq}^{(N)}$ is related to the quark mass contribution to the
nucleon mass, whose value can be found in refs. [50, 51, 52, 53, 54, 55, 56,
57, 58, 59]. $y_{\psi_{1}}^{\rm eff}$ is the effective Yukawa coupling of
$\psi_{1}$, and it can be obtained by the replacement of
$m_{\psi_{1}}\leftrightarrow m_{\mu}$ and $y_{L}\leftrightarrow y_{R}$ in the
expression of $y_{\mu}^{\rm eff}$ given in Eq. (C.1). This effective Yukawa
coupling increases when $m_{\psi_{1}}$ becomes large, because it is
proportional to $m_{\psi_{1}}$ like $y_{\mu}^{\rm eff}$ (see Eq. (2.26)). The
function $\hat{a}_{\psi_{1}}(x,y)$ is enhanced when $x\to 1$ with $y\neq 0$.
Therefore, the limit of $m_{\psi_{1}}\simeq m_{\varphi_{1}^{+}}$ leads to a
large contribution to the cross section from $a_{\psi_{1}}$. Note that for the
Majorana DM model, there are other contributions through the $Z$-penguin which
lead to effective interactions such as
$(\overline{\psi_{1}}\gamma^{\mu}\gamma^{5}\psi_{1})(\bar{q}\gamma_{\mu}q)$
and
$(\overline{\psi_{1}}\gamma^{\mu}\gamma^{5}\psi_{1})(\bar{q}\gamma_{\mu}\gamma^{5}q)$.
However, these contributions are suppressed by the lepton mass (and the DM
velocity for the former interaction), and we neglect their effects in our
analysis. Using the effective couplings in Eq. (3.7), the differential cross
section with respect to the recoil energy $E_{R}$ is estimated as

$\displaystyle\frac{d\sigma}{dE_{R}}$
$\displaystyle=\left\\{\frac{2m_{N}f_{A}^{2}}{\pi v^{2}}+\frac{4\alpha
Z^{2}}{v^{2}}a_{\psi_{1}}^{2}\left[2m_{N}v^{2}-\frac{(m_{N}+m_{\psi_{1}})^{2}}{m_{\psi_{1}}^{2}}E_{R}\right]\right\\}|F_{\rm
Helm}(E_{R})|^{2}$ $\displaystyle\hskip
11.99998pt+\frac{8m_{N}^{2}\mu_{A}^{2}}{\pi
v^{2}}E_{R}a_{\psi_{1}}^{2}\frac{J_{A}+1}{3J_{A}}|F_{\rm spin}(E_{R})|^{2}\,,$
( 3.11)

where $v$ is the DM velocity, $\alpha$ is the fine structure constant,
$f_{A}=ZC_{S,p}+(A-Z)C_{S,n}$ with an atomic number $Z$ and a mass number $A$,
and $m_{N}$, $\mu_{A}$ and $J_{A}$ are the mass, magnetic moment and spin of
the nucleus, respectively. $F_{\rm Helm}(E_{R})$ and $F_{\rm spin}(E_{R})$
denote form factors found in refs. [60, 61]. It can be seen from Eq. (3.11)
that the anapole contribution is suppressed by the DM velocity $v$ or the
recoil energy $E_{R}$. On the other hand, there is no suppression for
contributions from the contact-type interactions. It is notable that in our
model, $C_{S,N}$ in $f_{A}$ is enhanced by $y_{\psi_{1}}^{\rm eff}$ due to the
absence of the tree-level muon Yukawa coupling. Recently, the LUX-ZEPLIN (LZ)
experiment has reported their first results for spin-independent (SI) and
spin-dependent (SD) DM-nucleon scattering cross sections [62]. The upper limit
on the SI cross section has been improved, compared with previous results from
the XENON1T [63, 64] and PandaX-4T [65, 66] experiments. The corresponding
cross sections in our model are given by [50, 67, 68]

$\displaystyle\sigma_{\rm SI}^{\rm scalar}$
$\displaystyle=\frac{4\mu_{N}^{2}}{\pi}\frac{|y_{\psi_{1}}^{\rm
eff}|^{2}}{2m_{h}^{4}}\frac{m_{N}^{2}}{v_{H}^{2}}\left(\sum_{q}f_{Tq}^{(N)}\right)^{2}\,,$
( 3.12) $\displaystyle\sigma_{\rm SI}^{\rm anapole}$
$\displaystyle=8\alpha|a_{\psi_{1}}|^{2}\mu_{N}^{2}v^{2}\,.$ ( 3.13)

Here, $\mu_{N}$ is the reduced mass for $m_{\psi_{1}}$ and the nucleon mass
$\simeq 0.939$ GeV.

##  4 Numerical analysis

In this section, we first summarize the independent parameters in our model.
Then, the parameter space that gives the correct DM relic density and explains
the muon $g-2$ anomaly is identified and the size of the muon EDM is indicated
in that region as well as more general parameter regions. We take account of
muon coupling constraints presented in section 2.3 and also discuss
constraints from DM direct and indirect detection experiments as well as
collider searches.

###  4.1 Independent parameters

The Lagrangian of our model contains 18 parameters,

$\begin{split}&y_{\phi},y_{\eta},|m_{D}|,|m_{LL}|,|m_{RR}|,\theta_{\rm
phys},\\\
&m_{H}^{2},m_{\phi}^{2},m_{\eta}^{2},a,\lambda_{H},\lambda_{\phi},\lambda_{\eta},\lambda_{H\phi},\lambda_{H\eta},\lambda_{\phi\eta},\lambda^{\prime}_{H\phi},\lambda^{\prime\prime}_{H\phi}\,.\end{split}$
( 4.1)

Note that some of them are irrelevant to our analysis on the calculation of
the muon $g-2$, the muon EDM, the radiative mass, and the effective Yukawa
coupling of the muon. The Higgs mass-squared parameter $m_{H}^{2}$ is fixed by
the minimization condition in Eq. (2.5), and $\lambda_{H}$ should be
determined so that the SM Higgs mass, $m_{h}=125.25$ GeV, is correctly
reproduced. The quartic couplings $\lambda_{\phi}$, $\lambda_{\eta}$ and
$\lambda_{\phi\eta}$ are irrelevant to the mass spectrum of exotic particles,
although these values should be consistent with perturbative unitarity bounds
(commented below) and also chosen to avoid an unstable minimum of the scalar
potential. Moreover, $y_{\phi}y_{\eta}$ can be fixed by using Eq. (2.25), but
we need to check that values of the couplings do not exceed $\sqrt{4\pi}$. As
a result, the relevant (and independent) input parameters for the analysis can
be read as

$\displaystyle y_{\phi},|m_{D}|,|m_{LL}|,|m_{RR}|,\theta_{\rm
phys},m_{\phi}^{2},m_{\eta}^{2},a,\lambda_{H\phi},\lambda^{\prime}_{H\phi},\lambda^{\prime\prime}_{H\phi},\lambda_{H\eta}\,.$
( 4.2)

Note that $\lambda^{\prime}_{H\phi}$ and $\lambda^{\prime\prime}_{H\phi}$ are
relevant only to the masses of heavy neutral scalars, $m_{\sigma_{\phi}}^{2}$
and $m_{a_{\phi}}^{2}$ (see Eq. (2.17)), and irrelevant to our following
analysis as long as the DM candidate of the model is $\psi_{1}$. Furthermore,
we discuss our results by using $M_{\phi}^{2}$ and $M_{\eta}^{2}$ instead of
$m_{\phi}^{2},m_{\eta}^{2},\lambda_{H\phi}$ and $\lambda_{H\eta}$ (see below
Eq. (2.17)).

We here comment on perturbative unitarity bounds [69, 70], which are related
to $2\to 2$ scattering processes of scalar particles. At the tree level, it is
clear that quartic couplings are related to their amplitudes. In addition,
trilinear couplings also contribute to them through $s$-, $t$\- and
$u$-channel processes if scalar particles are not so heavy. There are studies
on the bounds, e.g., for models extended by singlet scalars [71, 72, 73] and
doublet scalars [74, 75, 76, 77, 78, 79, 80]. Since our model is a hybrid
extension with one singlet and one doublet scalars, there are lots of
scattering processes like $hh\to\sigma_{\phi}\sigma_{\phi}$,
$h\sigma_{\phi}\to\varphi_{i}^{+}\varphi_{j}^{-}$ and
$\varphi_{i}^{+}\varphi_{j}^{-}\to a_{\phi}a_{\phi}$. To obtain perturbative
unitarity bounds in our model, we use the SARAH/SPheno framework [81, 82, 83,
84, 85, 86, 87]. The details of the calculation for general scalar couplings
can be found in ref. [88].

###  4.2 Results

Figure 3: The predictions of $(g-2)_{\mu}$ and the muon EDM $d_{\mu}$ in terms
of $m_{\psi_{1}}$ and $m_{\varphi_{1}^{+}}$. We fix the model parameters as
follows: $y_{\phi}=1.2$, $m_{D}=700\,{\rm GeV}$, $m_{RR}=1000\,{\rm GeV}$,
$\theta_{\rm phys}=1.0$, $M_{\phi}^{2}=(1000\,{\rm GeV})^{2}$, $a=900\,{\rm
GeV}$. $m_{LL}$ and $M_{\eta}^{2}$ are changed as $320\,{\rm GeV}\leq
m_{LL}\leq 1200\,{\rm GeV}$ and $(540\,{\rm GeV})^{2}\leq
M_{\eta}^{2}\leq(1000\,{\rm GeV})^{2}$, respectively. Green and yellow shaded
regions denote $1\sigma$ and $2\sigma$ bands for $(g-2)_{\mu}$, and black
lines correspond to contours for $d_{\mu}$ in $10^{-23}\,e\,{\rm cm}$ unit.
The future prospects of the Fermilab Muon $g-2$ [8] and J-PARC Muon $g-2$/EDM
[9] experiments are shown as the orange shaded region. The red band is the
correct DM relic density, $\Omega h^{2}=0.120\pm 0.001$, and contours for
$\Omega h^{2}=0.5$ and $1.0$ are shown in blue and turquoise lines,
respectively. The dot-dashed magenta lines correspond to $\kappa_{\mu}=0.996$
and 1.004. The dashed gray line gives $m_{\psi_{1}}=m_{\varphi_{1}^{+}}$.

Fig. 3 shows the predictions of $(g-2)_{\mu}$ and the muon EDM $d_{\mu}$ in
our model. Here, we fix the relevant parameters as 555If we change $m_{RR}$
instead of $m_{LL}$ and set $m_{LL}=1000$ GeV, the prediction for $d_{\mu}$ is
totally the same as that of Fig. 3, because the mass eigenvalues of
$\psi_{1,2}$ and the mixing angle $\alpha$ are symmetric under
$m_{LL}\leftrightarrow m_{RR}$. Although $\tan\tau$ changes its sign, it is
irrelevant to the absolute value of $d_{\mu}$.

$\begin{split}&y_{\phi}=1.2,\quad M_{\phi}^{2}=(1000\,{\rm GeV})^{2},\quad
a=900\,{\rm GeV},\\\ &m_{D}=700\,{\rm GeV},\quad m_{RR}=1000\,{\rm
GeV},\quad\theta_{\rm phys}=1.0,\end{split}$ ( 4.3)

and change $m_{LL}$ and $M_{\eta}^{2}$ as $320\,{\rm GeV}\leq m_{LL}\leq
1200\,{\rm GeV}$ and $(540\,{\rm GeV})^{2}\leq M_{\eta}^{2}\leq(1000\,{\rm
GeV})^{2}$, respectively. For this parameter choice, the lightest particle
among $X$-odd particles is either $\psi_{1}$ or $\varphi_{1}^{+}$. The other
parameters which are not shown in Eq. (4.2), i.e. scalar quartic couplings,
does not affect the analysis here and are taken to be moderate values to
satisfy perturbative unitarity bounds. Note that when the scalar trilinear
coupling $a$ becomes large, some of the quartic couplings should be
$\mathcal{O}(1)$ to avoid the instability of the vacuum to give the correct
electroweak symmetry breaking. We numerically check that the SM vacuum is
stable when all of quartic couplings are within the range of $0.2$-$0.5$ with
the parameter choice shown above. These values of quartic couplings also
satisfy perturbative unitarity bounds, which is checked by the SARAH/SPheno
framework. In addition, since we have an additional $SU(2)_{L}$ doublet
scalar, there is a new contribution to the $T$-parameter [89, 90]. We have
calculated the contribution by following refs. [91, 92] and found that our
parameter choice leads to $\Delta T\sim 0.002$, which satisfies the current
constraint [93].

The current discrepancy of $(g-2)_{\mu}$ is [1, 2, 3, 4]

$\displaystyle\Delta a_{\mu}=(2.51\pm 0.59)\times 10^{-9}\,,$ ( 4.4)

whose $1\sigma$ and $2\sigma$ bands are shown as green and yellow shaded
regions in Fig. 3. Note that a lighter $m_{\psi_{1}}$ predicts a larger
$\Delta a_{\mu}$ due to its dependence, $\Delta a_{\mu}\sim
1/m_{\psi_{1}}^{2}$. In the figure, black lines correspond to contours for
$d_{\mu}$ in $10^{-23}\,e\,{\rm cm}$ unit. The future prospect of the muon
EDM, which is reported as $\mathcal{O}(10^{-21})\,e\,{\rm cm}$ at the Fermilab
Muon $g-2$ experiment [8] and the J-PARC Muon $g-2$/EDM experiment [9], is
shown as the orange shaded region. The red band shows the parameter space
where the correct DM relic density, $\Omega h^{2}=0.120\pm 0.001$ [94], is
obtained. Outside of this band, the relic density changes rapidly, as one can
see from the blue and turquoise contours which correspond to $\Omega
h^{2}=0.5$ and $1.0$, respectively. Note that in the whole parameter space of
Fig. 3, a new physics contribution to the ratio between the decay widths of
$Z\to e^{+}e^{-}$ and $Z\to\mu^{+}\mu^{-}$ in Eq. (2.49) is sufficiently
small, and we obtain $|\delta_{\mu\mu}|\lesssim 4\times 10^{-4}$ which is
consistent with the current data (2.47).

For the case of $m_{\psi_{1}}>m_{\varphi_{1}^{+}}$ (the region below the
dashed gray line), the $\psi_{1}\to\varphi_{1}^{\pm}+\mu^{\mp}$ decay occurs
at the tree level, and therefore, $\psi_{1}$ cannot be a DM candidate 666Since
$m_{\psi_{1}}<m_{\sigma_{\phi},a_{\phi}}$ in the current input parameters (see
Eqs. (2.17) and (4.3)), the DM cannot decay into $\nu_{\mu}+\sigma_{\phi}$,
$\nu_{\mu}+a_{\phi}$ in the plotted region of Fig. 3.. Without any interaction
to break the exotic number symmetry $X$, $\varphi_{1}^{+}$ is a stable exotic
particle, which may be cosmologically dangerous. However, we can consider, for
example, an interaction with the right-handed electron,
$\overline{L_{L}^{\mu}}\phi^{\dagger}e_{R}$, to make $\varphi_{1}^{+}$ decay
into $\nu_{\mu}+e^{+}$ 777This lepton flavor violating (LFV) interaction does
not induce LFV processes such as $\mu\to e\gamma$ because of the muon number
symmetry $L_{\mu}$. However, the interaction with a sizable coupling may be
constrained by the muonium-antimuonium oscillation [95, 96, 97, 98, 99]
although we do not need a large coupling for our purpose.. Interestingly, the
parameter region predicts a large $d_{\mu}$ due to a small value of
$m_{\varphi_{1}^{+}}$ and may be also explored by Higgs coupling measurements
at future collider experiments. Ref. [100] summarizes future sensitivities for
the measurements of the SM Higgs couplings. In particular, the Future Circular
Collider (FCC) may be able to measure $\kappa_{\mu}$ with relative precision
of $\sim 0.4\%$ whose contours are shown as dot-dashed magenta lines in Fig. 3
888Here, we just assume that the central value of $\kappa_{\mu}$ is $1$ at
future collider experiments..

In the whole parameter region shown in the figure, the muon EDM is predicted
to be larger than the future sensitivity of the PSI muEDM experiment, $6\times
10^{-23}\,e\,{\rm cm}$ [10, 11, 12]. The $2\sigma$ discrepancy of
$(g-2)_{\mu}$ can be explained for $560\,{\rm GeV}<m_{\psi_{1}}<780\,{\rm
GeV}$, while only the region of $m_{\psi_{1}}\simeq m_{\varphi_{1}^{+}}$ is
favored for the correct DM relic density. With the current parameter choice,
the coannihilation process plays an important role in obtaining the correct
relic density.

For $m_{\psi_{1}}\simeq 850\mathchar 45\relax 860$ GeV, the DM sector
contribution to the muon EDM accidentally disappears. This behavior can be
understood as follows. In this region, $m_{LL}\simeq 1000$ GeV which means
$m_{LL}\simeq m_{RR}$ for our current setup. Eq. (2.12) tells us that
$\theta_{\rm phys}=0$ or $|m_{LL}|=|m_{RR}|$ can lead to $\tan\tau=0$, which
makes $y_{L}^{ia}$ and $y_{R}^{ia}$ real. Therefore, there is no contribution
to the muon EDM for $|m_{LL}|=|m_{RR}|$, even when the physical phase has a
non-zero value, $\theta_{\rm phys}\neq 0$.

We now comment on constraints from DM searches at colliders and DM direct and
indirect detection experiments.

  * 1.

Collider searches
At the Large Hadron Collider (LHC), we expect a pair production of exotic
charged scalars $\varphi_{1}$ decaying into muons and DM fermions:

$\displaystyle pp\to\varphi_{1}\varphi_{1}\to\mu\mu+\psi_{1}\psi_{1}\,,$ (
4.5)

whose signal is two muons plus a large missing energy. The signal is similar
to that of a pair production of sleptons decaying into leptons and a missing
energy. Then, the ATLAS [101, 102] and CMS [103] experiments put a lower bound
on the DM mass. However, it is less than 500 GeV [101], which is outside the
plot range of Fig. 3. Ref. [104] has performed the numerical analysis to
obtain a bound on the DM mass for the similar model, and it was found to be
200-300 GeV, depending on the size of the mixing angle $s_{\theta}$ in Eq.
(2.18) and the mass of $\varphi_{1}^{+}$. Ref. [102] has investigated the case
with $m_{\psi_{1}}\sim m_{\varphi_{1}^{+}}$ and put a lower bound on
$m_{\varphi_{1}^{+}}$ less than 250 GeV for
$m_{\varphi_{1}^{+}}-m_{\psi_{1}}=\mathcal{O}(10)$ GeV.

  * 2.

Indirect detection
As mentioned in Sec. 3, the annihilation cross section of our DM $\psi_{1}$ is
dominated by $\psi_{1}\psi_{1}\to\mu\bar{\mu}$. For the parameter region in
Fig. 3, we obtain the prediction of $\langle\sigma
v\rangle_{\mu\bar{\mu}}\simeq\mathcal{O}(10^{-27})$-$\mathcal{O}(10^{-28})\,{\rm
cm}^{3}/{\rm s}$. Ref. [105] has studied a constraint on the annihilation
cross section of a Majorana DM, whose annihilation modes are
$\psi_{1}\psi_{1}\to\ell\bar{\ell}\gamma$ and
$\psi_{1}\psi_{1}\to\gamma\gamma$. The combination of the thermally averaged
cross sections, $\langle\sigma v\rangle_{\mu\bar{\mu}\gamma}+2\langle\sigma
v\rangle_{\gamma\gamma}$, is constrained to be less than
$10^{-26}$-$10^{-27}\,{\rm cm}^{3}/{\rm s}$, depending on the DM mass. The
cross sections of these annihilation processes, however, are several orders of
magnitude smaller than that of $\langle\sigma v\rangle_{\mu\bar{\mu}}$ in our
model, and therefore, ref. [105] does not put a constraint on the parameter
region shown in Fig. 3. It is notable that due to the gauge invariance, we
should consider the $\psi_{1}\psi_{1}\to\ell\bar{\ell}\gamma$ process together
with the $\psi_{1}\psi_{1}\to\ell\bar{\ell}Z$ process. Such processes may be
able to be explored by the PAMELA anti-proton search [106], and there are
studies for Majorana DM models [107, 108], although they indicate that it is
difficult to observe a Majorana DM at current and future telescopes.

  * 3.

Direct detection
Our Majorana DM scattering with the nucleon is induced by interactions
presented in Eq. (3.7). For the current parameter set, we obtain the SI DM-
nucleon scattering cross section of $\mathcal{O}(10^{-47}\mathchar 45\relax
10^{-50})\,{\rm cm}^{2}$. This is smaller than the current limit from the LZ
experiment [62], which is $(1.5\mathchar 45\relax 2.4)\times 10^{-46}\,{\rm
cm}^{2}$ for the DM mass range in Fig. 3. The LZ experiment [62] also put
constraints on the SD DM-proton and DM-neutron scattering cross sections, but
both are weaker than that of the SI cross section, and therefore, no region of
the parameter space is excluded by direct detection experiments. With the
future sensitivity of the LZ experiment, the upper limit on the SI cross
section will be improved by one order of magnitude [109, 110], which is still
not sufficient to explore our parameter space. By the future sensitivity of
PandaX-4T with 5.6 tonne$\cdot$year exposure [111], we may be able to explore
the parameter space in Fig. 3. Their current limit on the SI DM-nucleon
scattering cross section [65] can be read as $(2.6\mathchar 45\relax
4.2)\times 10^{-46}\,{\rm cm}^{2}$ for the DM mass range of $550\,{\rm
GeV}\leq m_{\psi_{1}}\leq 860\,{\rm GeV}$, and hence, if the future limit is
improved by a few orders of magnitude, a heavy DM mass region with
$m_{\psi_{1}}\simeq m_{\varphi_{1}^{+}}$ will be explored at the PandaX-4T
experiment.

Finally, let us discuss how our new physics contributions to $\Delta a_{\mu}$
and $d_{\mu}$ depend on input parameter choices. First of all, a different
choice of $\theta_{\rm phys}$ can change our predictions for $\Delta a_{\mu}$
and $d_{\mu}$ shown in Fig. 3. It is expected that $d_{\mu}$ is maximized by
choosing $\theta_{\rm phys}\sim\pi/4$, because $\theta_{\rm phys}\to 0$ or
$\theta_{\rm phys}\to\pi/2$ leads to $d_{\mu}\approx 0$. On the other hand,
the contribution to $\Delta a_{\mu}$ does not have such a clear dependence on
$\theta_{\rm phys}$. Actually, both observables strongly depend on the input
parameter set of $(m_{D},m_{LL},m_{RR},\theta_{\rm phys})$. For example, if
$(m_{D},m_{LL},m_{RR})=(700\,{\rm GeV},200\,{\rm GeV},1000\,{\rm GeV})$ are
chosen, $\Delta a_{\mu}<0$ for $0<\theta_{\rm phys}\lesssim 0.42$ and the
$2\sigma$ deviation can be explained for $\pi/6<\theta_{\rm phys}<\pi/4$, and
$d_{\mu}$ is maximized around $\theta_{\rm phys}\simeq\pi/6$. Instead, if we
choose $(m_{D},m_{LL},m_{RR})=(500\,{\rm GeV},990\,{\rm GeV},1000\,{\rm
GeV})$, $\Delta a_{\mu}<0$ is predicted in almost all range of $\theta_{\rm
phys}$ and the $2\sigma$ deviation can be explained only around $\theta_{\rm
phys}\simeq 1.45$, and the peak of $d_{\mu}$ appears at $\theta_{\rm
phys}\simeq 1.35$. In any case, our prediction of the muon EDM is
$d_{\mu}>\mathcal{O}(10^{-22})\,e\,{\rm cm}$. These observables also depend on
the values of $M_{\phi}^{2},M_{\eta}^{2}$. As one can see in Fig. 3, the
predictions of $\Delta a_{\mu}$ and $d_{\mu}$ become small as
$m_{\varphi_{1}^{+}}$ increases. In contrast, a large $a$ enhances the
contributions by a few %. This is because $a$ is related to the difference
between $m_{\varphi_{1}^{+}}^{2}$ and $m_{\varphi_{2}^{+}}^{2}$ (see Eqs.
(2.19) and (2.20)), and hence, a larger $a$ leads to a slightly smaller
$m_{\varphi_{1}^{+}}^{2}$.

##  5 Conclusion

In the present paper, we have investigated a prediction for the muon EDM
obtained in a model of DM. As shown in Table 1, the radiative stability
approach has a clear advantage to enhance the muon EDM, and we focused on a
model in which the muon mass is generated radiatively. With appropriate
discrete symmetries, exotic particles, $\psi$, $\phi$ and $\eta$, have
couplings to the muon (and also to the SM Higgs doublet). In this model, one
of the complex phases in the couplings cannot be removed by any field
redefinition and provides a physical CP phase, which leads to a new
contribution to the muon EDM. $\psi_{L,R}$ are singlet under the SM gauge
groups and the lightest mode gives a candidate of the Majorana fermion DM.

We found that even when the DM mass is heavier than the current collider
bound, $m_{\psi_{1}}>500$ GeV, the model predicts a muon EDM larger than
$10^{-22}\,e\,{\rm cm}$ which can be tested at the PSI muEDM experiment. In
the parameter space where the discrepancy of the muon $g-2$ and the correct DM
relic density are explained at the same time, the model predicts
$d_{\mu}\simeq(4\mathchar 45\relax 5)\times 10^{-22}\,e\,{\rm cm}$. For the
case of $m_{\psi_{1}}>m_{\varphi_{1}^{+}}$, the muon EDM can be even larger,
$d_{\mu}\simeq(7\mathchar 45\relax 8)\times 10^{-22}\,e\,{\rm cm}$, due to a
small value of $m_{\varphi_{1}^{+}}$, although $\psi_{1}$ does not give a DM
candidate. Furthermore, once we forget a new physics explanation for the muon
$g-2$ discrepancy as well as the DM relic density, the muon EDM can be larger
than the future sensitivities of the ongoing Fermilab Muon $g-2$ and projected
J-PARC Muon $g-2$/EDM experiments.

One of the most promising approaches to probe our DM model is a future muon
collider (see e.g. ref. [112] and references therein) because a muon collider
is expected to have a new particle mass reach higher than that of the LHC and
also our DM fermion directly couples to the muon. It would be interesting to
explore the phenomenology of our DM model to generate the radiative muon mass,
the muon $g-2$ and the muon EDM at a muon collider, which is left for future
study.

## Acknowledgements

KSK is supported by Natural Science Foundation of China (NSFC) under grant No.
12050410233. YN is supported by NSFC under grant No. 12150610465.

##  Appendix A Neutrino sector

###  A.1 A scalar triplet extension

In our model, due to the muon number symmetry, we need a further extension for
obtaining the correct neutrino mixing angles. One of the simplest way to
reproduce them is to introduce a scalar triplet, as discussed in appendix A of
ref. [113]. At first, we can write down dimension-five operators which are
related to lepton doublets as

$\displaystyle-\frac{c_{HH}^{ab}}{M_{HH}}\left(\overline{L_{L}^{a}}\widetilde{H}\right)\left(\widetilde{H}^{T}(L_{L}^{c})^{b}\right)+{\rm
h.c.}\,,$ ( A.1)

where $\widetilde{H}=i\sigma_{2}H^{*}$, and $a,b=e,\mu,\tau$ are indices for
lepton species. Note that the exotic number symmetry forbids terms of
$\left(\overline{L_{L}^{a}}\widetilde{H}\right)\left(\widetilde{\phi}^{T}(L_{L}^{c})^{b}\right)$
and
$\left(\overline{L_{L}^{a}}\widetilde{\phi}\right)\left(\widetilde{H}^{T}(L_{L}^{c})^{b}\right)$,
and the term of
$\left(\overline{L_{L}^{a}}\widetilde{\phi}\right)\left(\widetilde{\phi}^{T}(L_{L}^{c})^{b}\right)$
is irrelevant to the discussion on the neutrino mixing angles, because $\phi$
does not acquire a nonzero VEV. The coefficient $c_{HH}$ is a $3\times 3$
matrix, but symmetries of the model makes it have the form,

$\displaystyle c_{HH}=\begin{pmatrix}c_{HH}^{ee}&0&c_{HH}^{e\tau}\\\
0&c_{HH}^{\mu\mu}&0\\\ c_{HH}^{\tau e}&0&c_{HH}^{\tau\tau}\end{pmatrix}\,.$ (
A.2)

Therefore, the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix obtained from our
model cannot be consistent with the experimental result at this stage.
However, once we introduce a $SU(2)_{L}$ triplet scalar $\Delta$ to the model,
which has the $U(1)_{Y}$ charge of $-1$ and odd muon number, we can write
additional terms of

$\displaystyle-
c_{\Delta}^{ab}\overline{L_{L}^{a}}\Delta(L_{L}^{c})^{b}\,,\quad
c_{\Delta}=\begin{pmatrix}0&c_{\Delta}^{e\mu}&0\\\ c_{\Delta}^{\mu
e}&0&c_{\Delta}^{\mu\tau}\\\ 0&c_{\Delta}^{\tau\mu}&0\\\ \end{pmatrix}\,.$ (
A.3)

When $\Delta$ acquires a nonzero VEV, $v_{\Delta}\neq 0$, we can reproduce all
elements for the neutrino mass matrix as

$\displaystyle
m_{\nu}=\begin{pmatrix}(m_{\nu})^{ee}&(m_{\nu})^{e\mu}&(m_{\nu})^{e\tau}\\\
(m_{\nu})^{\mu e}&(m_{\nu})^{\mu\mu}&(m_{\nu})^{\mu\tau}\\\ (m_{\nu})^{\tau
e}&(m_{\nu})^{\tau\mu}&(m_{\nu})^{\tau\tau}\end{pmatrix}\,.$ ( A.4)

Here, each $(m_{\nu})^{ab}$ is estimated by

$\displaystyle(m_{\nu})^{ab}=\begin{cases}{\displaystyle\frac{c_{HH}^{ab}}{2}\frac{v_{H}^{2}}{M_{HH}}}&\text{for
$(a,b)=(e,e),(e,\tau),(\mu,\mu),(\tau,e),(\tau,\tau)$}\\\\[10.33327pt]
{\displaystyle c_{\Delta}^{ab}v_{\Delta}}&\text{for
$(a,b)=(e,\mu),(\mu,e),(\mu,\tau),(\tau,\mu)$}\end{cases}\,.$ ( A.5)

Then, if $v_{\Delta}$ is $\mathcal{O}(v_{H}^{2}/M_{HH})$ with $c_{HH}^{ab}\sim
c_{\Delta}^{ab}$, all elements of Eq. (A.4) have the similar order, and hence,
the large mixing angles for the PMNS matrix can be obtained. For the neutrino
masses of $\mathcal{O}({\rm eV})$, the mass scale for $M_{HH}$ is required to
be $\mathcal{O}(10^{13{\rm-}14})$ GeV.

###  A.2 Right-handed neutrinos

Another possibility to reproduce the correct PMNS matrix is to introduce three
generations of the right-handed neutrinos (RHNs), denoted as
$N_{R}^{e,\mu,\tau}$. Similar to the charged lepton sector, only $N_{R}^{\mu}$
is odd under the muon number $L_{\mu}$, and we have additional Dirac Yukawa
couplings and Majorana mass terms for neutrinos as
($\ell^{(^{\prime})}=e,\tau$)

$\displaystyle\mathcal{L}_{N}=-y_{\nu}^{\ell\ell^{\prime}}\overline{L_{L}^{\ell}}\widetilde{H}N_{R}^{\ell^{\prime}}-y_{\nu}^{\mu\mu}\overline{L_{L}^{\mu}}\widetilde{H}N_{R}^{\mu}-m_{N}^{\ell\ell^{\prime}}\overline{N_{R}^{\ell\,c}}N_{R}^{\ell^{\prime}}-m_{N}^{\mu\mu}\overline{N_{R}^{\mu\,c}}N_{R}^{\mu}-m_{N}^{\ell\mu}\overline{N_{R}^{\ell\,c}}N_{R}^{\mu}+{\rm
h.c.}\,,$ ( A.6)

where the last term breaks the $L_{\mu}$ symmetry softly, which is required
for the correct neutrino mixing angles. This can be understood
diagrammatically, as shown in Fig. 4. In addition to this diagram, the mixing
between $\nu_{L}^{e}$ and $\nu_{L}^{\tau}$ is also induced by the same diagram
with changing $\mu\to e$ or $\tau$. The mass matrix of $\nu_{L}^{e,\mu,\tau}$
can be obtained as

$\displaystyle
m_{\nu_{L}}=-\frac{v_{H}^{2}}{2}\begin{pmatrix}y_{\nu}^{ee}&0&y_{\nu}^{e\tau}\\\
0&y_{\nu}^{\mu\mu}&0\\\ y_{\nu}^{\tau e}&0&y_{\nu}^{\tau\tau}\\\
\end{pmatrix}\begin{pmatrix}m_{N}^{ee}&m_{N}^{e\mu}&m_{N}^{e\tau}\\\
m_{N}^{\mu e}&m_{N}^{\mu\mu}&m_{N}^{\mu\tau}\\\ m_{N}^{\tau
e}&m_{N}^{\tau\mu}&m_{N}^{\tau\tau}\\\
\end{pmatrix}^{-1}\begin{pmatrix}y_{\nu}^{ee}&0&y_{\nu}^{\tau e}\\\
0&y_{\nu}^{\mu\mu}&0\\\ y_{\nu}^{e\tau}&0&y_{\nu}^{\tau\tau}\\\
\end{pmatrix}\,.$ ( A.7)

Therefore, if the mass scales of $m_{N}^{\ell\ell^{\prime}}$ are similar with
each other, we can obtain a full $3\times 3$ matrix for light neutrino states,
which can be consistent with experimental results on the PMNS matrix. Similar
to the previous method, the mass scale of $m_{N}^{\ell\ell^{\prime}}$ is
required to be $\mathcal{O}(10^{13{\rm-}14})$ GeV for the neutrino masses of
$\mathcal{O}({\rm eV})$ if $y_{\nu}^{\ell\ell^{\prime}}\sim\mathcal{O}(1)$.

$\nu_{L}^{\mu}$$\nu_{L}^{e,\tau}$$\langle H\rangle$$\langle
H\rangle$$\times$$N_{R}^{\mu}$$N_{R}^{e,\tau}$ Figure 4: A Feynman diagram for
the neutrino mixing between $\nu_{L}^{\mu}$ and $\nu_{L}^{e,\tau}$ where
“$\times$” indicates the soft breaking mixing between $N_{R}^{\mu}$ and
$N_{R}^{e,\tau}$.

###  A.3 Comment on LFV

For the second example, however, we have LFV processes due to the soft
$L_{\mu}$ breaking terms. To see this, we focus on the first two generations
of leptons, namely, the electron-muon system. From Eq. (A.6) and the SM Yukawa
interactions for the electron, we have one-loop contributions to off-diagonal
elements of Yukawa couplings of charged leptons, as shown in Fig. 5.

$e_{R}$$L_{L}^{\mu}$$H$$H$$L_{L}^{e}$$\nu_{R}^{\alpha}$ Figure 5: One-loop
contributions to off-diagonal elements of Yukawa couplings of charged leptons.
Here, $\nu_{R}^{\alpha}$ is the mass eigenstate of RHNs.

Due to the soft $L_{\mu}$ breaking, the mass eigenstate $\nu_{R}^{\alpha}$ can
be written by one mixing angle $\theta_{N}$ as

$\displaystyle\begin{pmatrix}\nu_{R}^{1}\\\
\nu_{R}^{2}\end{pmatrix}=\begin{pmatrix}\cos\theta_{N}&\sin\theta_{N}\\\
-\sin\theta_{N}&\cos\theta_{N}\end{pmatrix}\begin{pmatrix}N_{R}^{e}\\\
N_{R}^{\mu}\end{pmatrix}\,.$ ( A.8)

After integrating out the RHNs, we obtain the off-diagonal element,

$\displaystyle\mathcal{L}_{\rm eff}\supset-\delta y_{\mu
e}\overline{L_{L}^{\mu}}He_{R}+{\rm h.c.}\,,$ ( A.9)

where $\delta y_{\mu e}$ is obtained by Fig. 5 and roughly estimated as

$\displaystyle\delta y_{\mu
e}\simeq\frac{y_{e}}{16\pi^{2}}y_{\nu}^{ee}y_{\nu}^{\mu\mu}\sin\theta_{N}\cos\theta_{N}\mathop{\rm
ln}\frac{M_{\nu_{R}^{1}}^{2}}{M_{\nu_{R}^{2}}^{2}}\,,$ ( A.10)

with $y_{e}$ being the SM Yukawa coupling of the term
$\overline{L_{L}^{e}}He_{R}$. This off-diagonal element can be removed by
field redefinition of left-handed lepton doublets,

$\displaystyle\ell_{L}^{e}=L_{L}^{e}\cos\theta_{e\mu}+L_{L}^{\mu}\sin\theta_{e\mu}\,,\leavevmode\nobreak\
\ell_{L}^{\mu}=-L_{L}^{e}\sin\theta_{e\mu}+L_{L}^{\mu}\cos\theta_{e\mu}\,,$ (
A.11)

which lead to an electron coupling with exotic particles,

$\displaystyle y_{\phi}\overline{L_{L}^{\mu}}\phi^{\dagger}\psi_{R}\to
y_{\phi}^{e}\overline{\ell_{L}^{e}}\phi^{\dagger}\psi_{R}+y_{\phi}^{\mu}\overline{\ell_{L}^{\mu}}\phi^{\dagger}\psi_{R}\,,$
( A.12)

where we define $y_{\phi}^{e}\equiv y_{\phi}\sin\theta_{e\mu}$ and
$y_{\phi}^{\mu}\equiv y_{\phi}\cos\theta_{e\mu}$. Then, we have a one-loop
contribution to $\mu\to e\gamma$ by replacing $\mu_{L}\to e_{L}$ in the
diagram (b) of Fig. 1. Its branching ratio can be calculated as [114]

$\displaystyle{\rm BR}(\mu\to
e\gamma)=\frac{3(4\pi)^{3}\alpha}{2G_{F}^{2}m_{\mu}^{2}}\Bigl{(}\left|C_{T}^{e\mu}\right|^{2}+\left|C_{T^{\prime}}^{e\mu}\right|^{2}\Bigr{)}\,{\rm
BR}(\mu\to e\nu_{\mu}\bar{\nu}_{e})\,,$ ( A.13)

where $G_{F}$ is the Fermi constant, BR$(\mu\to
e\nu_{\mu}\bar{\nu}_{e})\approx 1$, and $C_{T}^{e\mu}$ and
$C_{T^{\prime}}^{e\mu}$ are coefficients of dipole operators, whose
definitions are

$\displaystyle\mathcal{L}_{\rm
LFV}\supset-\frac{e}{2}C_{T}^{e\mu}\left(\bar{e}\sigma^{\alpha\beta}\mu\right)F_{\alpha\beta}-\frac{e}{2}C_{T^{\prime}}^{e\mu}\left(\bar{e}i\sigma^{\alpha\beta}\gamma_{5}\mu\right)F_{\alpha\beta}\,.$
( A.14)

From Eqs. (2.28), (2.33) and (2.34), leading contributions to $C_{T}^{e\mu}$
and $C_{T^{\prime}}^{e\mu}$ can be easily estimated by replacing
$y_{L}^{ia}\to y_{L}^{ia}\sin\theta_{e\mu}$. Similar to this replacement, our
predictions of $a_{\mu}$ and $d_{\mu}$ in Eqs. (2.31) and (2.32) are changed
as

$\displaystyle a_{\mu}$
$\displaystyle=2m_{\mu}\left(C_{T}(0)\cos\theta_{\mu}+C_{T^{\prime}}(0)\sin\theta_{\mu}\right)\cos\theta_{e\mu}\,,$
( A.15) $\displaystyle d_{\mu}$
$\displaystyle=e\left(C_{T^{\prime}}(0)\cos\theta_{\mu}-C_{T}(0)\sin\theta_{\mu}\right)\cos\theta_{e\mu}\,.$
( A.16)

From these facts, we have the following relation between $C_{T}^{e\mu}$,
$C_{T^{\prime}}^{e\mu}$ and $a_{\mu}$, $d_{\mu}$ as

$\displaystyle|C_{T}^{e\mu}|^{2}+|C_{T^{\prime}}^{e\mu}|^{2}$
$\displaystyle=\left(\left|C_{T}(0)c_{\theta_{\mu}/2}+C_{T^{\prime}}(0)s_{\theta_{\mu}/2}\right|^{2}+\left|C_{T^{\prime}}(0)c_{\theta_{\mu}/2}-C_{T}(0)s_{\theta_{\mu}/2}\right|^{2}\right)\sin^{2}\theta_{e\mu}$
$\displaystyle=\left(|C_{T}(0)|^{2}+|C_{T^{\prime}}(0)|^{2}\right)\sin^{2}\theta_{e\mu}$
$\displaystyle=\left(\left|\frac{a_{\mu}}{2m_{\mu}}\right|^{2}+\left|\frac{d_{\mu}}{e}\right|^{2}\right)\tan^{2}\theta_{e\mu}\,,$
( A.17)

and hence, the branching ratio in Eq. (A.13) can be expressed by $a_{\mu}$ and
$d_{\mu}$ as

$\displaystyle{\rm BR}(\mu\to
e\gamma)=\frac{3(4\pi)^{3}\alpha}{2G_{F}^{2}m_{\mu}^{2}}\left(\left|\frac{a_{\mu}}{2m_{\mu}}\right|^{2}+\left|\frac{d_{\mu}}{e}\right|^{2}\right)\tan^{2}\theta_{e\mu}\,.$
( A.18)

By assuming $a_{\mu}=2.51\times 10^{-9}$ and $d_{\mu}=4.5\times
10^{-22}\,e\,{\rm cm}$ as we found in our model, the upper limit on
$\tan\theta_{e\mu}$ can be obtained as

$\displaystyle\tan\theta_{e\mu}\lesssim 6.7\times 10^{-6}\,,$ ( A.19)

where we have used the current upper bound on the branching ratio, BR$(\mu\to
e\gamma)<4.2\times 10^{-13}$ [115]. Then, the mixing angle in Eq. (A.11)
should be tiny, which means a small off-diagonal element compared to $y_{e}$,
$\delta y_{\mu e}\ll y_{e}$. To realize the constraint in Eq. (A.19), we
roughly need $\delta y_{\mu e}/y_{e}\sim 6\times 10^{-6}$. Assuming
$\mathop{\rm ln}(M_{\nu_{R}^{1}}^{2}/M_{\nu_{R}^{2}}^{2})\sim\mathcal{O}(1)$
and $\sin\theta_{N}\sim\cos\theta_{N}\sim 1/\sqrt{2}$ due to the similar order
for all $m_{N}^{\ell\ell^{\prime}}$ to obtain large mixing angles for
neutrinos, we obtain $y_{\nu}^{ee}y_{\nu}^{\mu\mu}\sim 2\times 10^{-3}$, and
the scale of RHNs will be $\mathcal{O}(10^{10})$ GeV for the light neutrino
masses of $\mathcal{O}({\rm eV})$.

##  Appendix B Loop integrals

In the appendix, we summarize loop integrals relevant to our analysis. Note
that we use dimensional regularization, and
$\Delta_{\epsilon}\equiv\frac{2}{\epsilon}-\gamma_{E}+\mathop{\rm ln}4\pi$
with the Euler constant $\gamma_{E}$ and $\epsilon=4-D$ diverges when $D\to
4$.

  * – Self-energy integral –

$\displaystyle B_{0}(p^{2},m_{0}^{2},m_{1}^{2})$
$\displaystyle=\Delta_{\epsilon}-\int_{0}^{1}\\!dx\mathop{\rm
ln}\left[\frac{-x(1-x)p^{2}+xm_{1}^{2}+(1-x)m_{0}^{2}}{\mu^{2}}\right]\,.$ (
B.1)

When $p^{2}=0$ and $m_{0}^{2}=m_{1}^{2}\equiv m^{2}$, we can obtain the
following simple form:

$\displaystyle B_{0}(0,m^{2},m^{2})$
$\displaystyle=\Delta_{\epsilon}-\mathop{\rm ln}\frac{m^{2}}{\mu^{2}}\,.$ (
B.2)

Note that $B_{0}$ has a divergent part which should be cancelled for physical
predictions.

  * – Triangle integrals –

$\displaystyle
C_{0}(p_{1}^{2},p_{12}^{2},p_{2}^{2},m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=-\int_{0}^{1}\\!dx_{1}\int_{0}^{1-x_{1}}\\!dx_{2}\frac{1}{C_{3}}\,,$
( B.3) $\displaystyle
C_{1}(p_{1}^{2},p_{12}^{2},p_{2}^{2},m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=\int_{0}^{1}\\!dx_{1}\int_{0}^{1-x_{1}}\\!dx_{2}\frac{x_{1}}{C_{3}}\,,$
( B.4) $\displaystyle
C_{ij}(p_{1}^{2},p_{12}^{2},p_{2}^{2},m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=-\int_{0}^{1}\\!dx_{1}\int_{0}^{1-x_{1}}\\!dx_{2}\frac{x_{i}x_{j}}{C_{3}}\,,$
( B.5)

where
$C_{3}=x_{1}(x_{1}-1)p_{1}^{2}+x_{2}(x_{2}-1)p_{2}^{2}+x_{1}x_{2}(p_{1}^{2}+p_{2}^{2}-p_{12}^{2})+x_{1}m_{1}^{2}+x_{2}m_{2}^{2}+(1-x_{1}-x_{2})m_{0}^{2}$
with $p_{12}^{2}=(p_{1}-p_{2})^{2}$. When all of $p_{1}^{2}$, $p_{2}^{2}$ and
$p_{12}^{2}$ are $0$, these integrals are simplified as follows. First, the
$C_{0}$ function is (here we use $r_{i}\equiv m_{i}^{2}/m_{0}^{2}$):

$\displaystyle C_{0}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=\frac{1}{m_{0}^{2}}\frac{r_{1}\mathop{\rm
ln}r_{1}-r_{1}r_{2}\mathop{\rm ln}r_{1}-r_{2}\mathop{\rm
ln}r_{2}+r_{1}r_{2}\mathop{\rm ln}r_{2}}{(1-r_{1})(1-r_{2})(r_{1}-r_{2})}\,,$
( B.6) $\displaystyle C_{0}(0,0,0,m_{0}^{2},m_{1}^{2},m_{1}^{2})$
$\displaystyle=\frac{1}{m_{0}^{2}}\frac{1-r_{1}+\mathop{\rm
ln}r_{1}}{(1-r_{1})^{2}}\,,$ ( B.7) $\displaystyle
C_{0}(0,0,0,m_{0}^{2},m_{0}^{2},m_{0}^{2})$
$\displaystyle=-\frac{1}{2m_{0}^{2}}\,.$ ( B.8)

Note that $C_{0}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})$ is symmetric under
exchanging two of $m_{i}^{2}$, e.g.,
$C_{0}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})=C_{0}(0,0,0,m_{1}^{2},m_{0}^{2},m_{2}^{2})=C_{0}(0,0,0,m_{0}^{2},m_{2}^{2},m_{1}^{2})$.

Second, the $C_{1}$ function is simplified as

$\displaystyle C_{1}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=-\frac{1}{2m_{0}^{2}}\left[\frac{r_{1}}{(1-r_{1})(r_{1}-r_{2})}+\frac{r_{1}(r_{1}-2r_{2}+r_{1}r_{2})}{(1-r_{1})^{2}(r_{1}-r_{2})^{2}}\mathop{\rm
ln}r_{1}\right.$ $\displaystyle\hskip
60.00009pt\left.+\frac{r_{2}^{2}}{(1-r_{2})(r_{1}-r_{2})^{2}}\mathop{\rm
ln}r_{2}\right]\,,$ ( B.9) $\displaystyle
C_{1}(0,0,0,m_{0}^{2},m_{1}^{2},m_{1}^{2})$
$\displaystyle=-\frac{1}{4m_{0}^{2}}\frac{r_{1}^{2}-4r_{1}+3+2\mathop{\rm
ln}r_{1}}{(1-r_{1})^{3}}\,,$ ( B.10) $\displaystyle
C_{1}(0,0,0,m_{0}^{2},m_{0}^{2},m_{0}^{2})$
$\displaystyle=\frac{1}{6m_{0}^{2}}\,.$ ( B.11)

Note that there is a symmetric property only for $m_{0}^{2}\leftrightarrow
m_{2}^{2}$,
$C_{1}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})=C_{1}(0,0,0,m_{2}^{2},m_{1}^{2},m_{0}^{2})$.

Third, the $C_{11}$ function is

$\displaystyle C_{11}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=-\frac{1}{3m_{0}^{2}}\left[\frac{r_{1}(r_{1}^{2}-3r_{1}+5r_{2}-3r_{1}r_{2})}{2(1-r_{1})^{2}(r_{1}-r_{2})^{2}}\right.$
$\displaystyle\hskip
40.00006pt-\frac{r_{1}(r_{1}^{2}r_{2}^{2}+r_{1}^{2}r_{2}+r_{1}^{2}-3r_{1}r_{2}^{2}-3r_{1}r_{2}+3r_{2}^{2})}{(1-r_{1})^{3}(r_{1}-r_{2})^{3}}\mathop{\rm
ln}r_{1}$ $\displaystyle\hskip
40.00006pt\left.+\frac{r_{2}^{3}}{(1-r_{2})(r_{1}-r_{2})^{3}}\mathop{\rm
ln}r_{2}\right]\,,$ ( B.12) $\displaystyle
C_{11}(0,0,0,m_{0}^{2},m_{1}^{2},m_{1}^{2})$
$\displaystyle=-\frac{1}{18m_{0}^{2}}\frac{2r_{1}^{3}-9r_{1}^{2}+18r_{1}-11-6\mathop{\rm
ln}r_{1}}{(1-r_{1})^{4}}\,,$ ( B.13) $\displaystyle
C_{11}(0,0,0,m_{0}^{2},m_{0}^{2},m_{0}^{2})$
$\displaystyle=-\frac{1}{12m_{0}^{2}}\,.$ ( B.14)

Finally, the $C_{12}$ function is

$\displaystyle C_{12}(0,0,0,m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=-\frac{1}{6m_{0}^{2}}\left[\frac{r_{1}^{2}r_{2}-r_{1}^{2}+r_{1}r_{2}^{2}-r_{2}^{2}}{(1-r_{1})(1-r_{2})(r_{1}-r_{2})^{2}}\right.$
$\displaystyle\hskip
60.00009pt-\frac{r_{1}^{2}(2r_{1}r_{2}+r_{1}-3r_{2})}{(1-r_{1})^{2}(r_{1}-r_{2})^{3}}\mathop{\rm
ln}r_{1}$ $\displaystyle\hskip
60.00009pt\left.+\frac{r_{2}^{2}(2r_{1}r_{2}-3r_{1}+r_{2})}{(1-r_{2})^{2}(r_{1}-r_{2})^{3}}\mathop{\rm
ln}r_{2}\right]\,,$ ( B.15) $\displaystyle
C_{12}(0,0,0,m_{0}^{2},m_{1}^{2},m_{1}^{2})$
$\displaystyle=-\frac{1}{36m_{0}^{2}}\frac{2r_{1}^{3}-9r_{1}^{2}+18r_{1}-11-6\mathop{\rm
ln}r_{1}}{(1-r_{1})^{4}}\,,$ ( B.16) $\displaystyle
C_{12}(0,0,0,m_{0}^{2},m_{0}^{2},m_{0}^{2})$
$\displaystyle=-\frac{1}{24m_{0}^{2}}\,.$ ( B.17)

##  Appendix C Full forms for $y_{\mu}^{\rm eff}$, $C_{T}$ and
$C_{T^{\prime}}$

We here present full expressions for $y_{\mu}^{\rm eff}$, $C_{T}$ and
$C_{T^{\prime}}$ at the one-loop order. In their calculations, FeynCalc [116,
117, 118] is used, and to deal with $\gamma_{5}$ in dimensional
regularization, we adopt the ’t Hooft-Veltman-Breitenlohner-Maison (BMHV)
prescription [119, 120]. The results are

$\displaystyle y_{\mu}^{\rm eff}(p_{h^{0}}^{2})$
$\displaystyle=\sum_{i,j,a}\Biggl{\\{}\Biggl{.}-\frac{y_{L}^{ia}y_{R}^{ja}A_{ij}}{16\pi^{2}}m_{\psi_{a}}C_{0}(m_{\mu}^{2},m_{\mu}^{2},p_{h^{0}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{j}^{+}}^{2})$
$\displaystyle\hskip 40.00006pt+\frac{A_{ij}}{16\pi^{2}}m_{\mu}^{\rm
rad}\left[y_{R}^{ia\,*}y_{R}^{ja}C_{1}(m_{\mu}^{2},p_{h^{0}}^{2},m_{\mu}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\varphi_{j}^{+}}^{2})\right.$
$\displaystyle\hskip
100.00015pt\left.+y_{L}^{ia}y_{L}^{ja\,*}C_{2}(m_{\mu}^{2},p_{h^{0}}^{2},m_{\mu}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2},m_{\varphi_{j}^{+}}^{2})\right]\Biggl{.}\Biggr{\\}}\,,$
( C.1) $\displaystyle C_{T}(q^{2})$
$\displaystyle=\sum_{i,a}\Biggl{\\{}\Biggr{.}\frac{{\rm
Re}[y_{L}^{ia}y_{R}^{ia}]}{16\pi^{2}}m_{\psi_{a}}\Bigl{[}\Bigr{.}Q_{S}\left(C_{0}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})+2C_{1}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\right)$
$\displaystyle\hskip
140.00021pt-2Y_{\psi}C_{1}(q^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2})\Bigl{.}\Bigr{]}$
$\displaystyle\hskip
40.00006pt-\frac{\left|y_{L}^{ia}\right|^{2}+\left|y_{R}^{ia}\right|^{2}}{16\pi^{2}}{\rm
Re}[m_{\mu}^{\rm rad}]\Bigl{[}\Bigr{.}Q_{S}C_{\rm
sub}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})$ $\displaystyle\hskip
180.00027pt+Y_{\psi}C_{\rm
sub}(q^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2})\Bigl{.}\Bigr{]}\Biggl{.}\Biggr{\\}}\,,$
( C.2) $\displaystyle C_{T^{\prime}}(q^{2})$
$\displaystyle=\sum_{i,a}\Biggl{\\{}\Biggr{.}\frac{{\rm
Im}[y_{L}^{ia}y_{R}^{ia}]}{16\pi^{2}}m_{\psi_{a}}\Bigl{[}\Bigr{.}Q_{S}\left(C_{0}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})+2C_{1}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\right)$
$\displaystyle\hskip
140.00021pt-2Y_{\psi}C_{1}(q^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2})\Bigl{.}\Bigr{]}$
$\displaystyle\hskip
40.00006pt-\frac{\left|y_{L}^{ia}\right|^{2}+\left|y_{R}^{ia}\right|^{2}}{16\pi^{2}}{\rm
Im}[m_{\mu}^{\rm rad}]\Bigl{[}\Bigr{.}Q_{S}C_{\rm
sub}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})$ $\displaystyle\hskip
180.00027pt+Y_{\psi}C_{\rm
sub}(q^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2})\Bigl{.}\Bigr{]}\Biggl{.}\Biggr{\\}}\,,$
( C.3)

where $Q_{S}=1+Y_{\psi}$, and we have defined

$\displaystyle C_{0}(q^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})$
$\displaystyle\equiv
C_{0}(m_{\mu}^{2},m_{\mu}^{2},q^{2},m_{\varphi_{i}^{+}}^{2},m_{\psi_{a}}^{2},m_{\varphi_{i}^{+}}^{2})\,,$
( C.4) $\displaystyle C_{N}(q^{2},m_{A}^{2},m_{B}^{2})$ $\displaystyle\equiv
C_{N}(m_{\mu}^{2},q^{2},m_{\mu}^{2},m_{A}^{2},m_{B}^{2},m_{B}^{2})\qquad(N=1,11,12)\,,$
( C.5)

and for sub-leading contributions, we have also defined

$\displaystyle C_{\rm sub}(q^{2},m_{A}^{2},m_{B}^{2})$ $\displaystyle\equiv
C_{1}(q^{2},m_{A}^{2},m_{B}^{2})+C_{11}(q^{2},m_{A}^{2},m_{B}^{2})+C_{12}(q^{2},m_{A}^{2},m_{B}^{2})\,.$
( C.6)

## References

  * [1] Muon g-2 Collaboration, G. W. Bennett et al., “Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL,” Phys. Rev. D 73 (2006) 072003, arXiv:hep-ex/0602035.
  * [2] A. Keshavarzi, D. Nomura, and T. Teubner, “Muon $g-2$ and $\alpha(M_{Z}^{2})$: a new data-based analysis,” Phys. Rev. D 97 no. 11, (2018) 114025, arXiv:1802.02995 [hep-ph].
  * [3] T. Aoyama et al., “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887 (2020) 1–166, arXiv:2006.04822 [hep-ph].
  * [4] Muon g-2 Collaboration, B. Abi et al., “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm,” Phys. Rev. Lett. 126 no. 14, (2021) 141801, arXiv:2104.03281 [hep-ex].
  * [5] A. Keshavarzi, K. S. Khaw, and T. Yoshioka, “Muon $g-2$: A review,” Nucl. Phys. B 975 (2022) 115675, arXiv:2106.06723 [hep-ex].
  * [6] Muon (g-2) Collaboration, G. W. Bennett et al., “An Improved Limit on the Muon Electric Dipole Moment,” Phys. Rev. D 80 (2009) 052008, arXiv:0811.1207 [hep-ex].
  * [7] Y. Ema, T. Gao, and M. Pospelov, “Improved Indirect Limits on Muon Electric Dipole Moment,” Phys. Rev. Lett. 128 no. 13, (2022) 131803, arXiv:2108.05398 [hep-ph].
  * [8] Muon g-2 Collaboration, R. Chislett, “The muon EDM in the g-2 experiment at Fermilab,” EPJ Web Conf. 118 (2016) 01005.
  * [9] M. Abe et al., “A New Approach for Measuring the Muon Anomalous Magnetic Moment and Electric Dipole Moment,” PTEP 2019 no. 5, (2019) 053C02, arXiv:1901.03047 [physics.ins-det].
  * [10] A. Adelmann et al., “Search for a muon EDM using the frozen-spin technique,” arXiv:2102.08838 [hep-ex].
  * [11] M. Sakurai et al., “muEDM: Towards a search for the muon electric dipole moment at PSI using the frozen-spin technique,” in 24th International Symposium on Spin Physics.  1, 2022.  arXiv:2201.06561 [hep-ex].
  * [12] muon EDM initiative Collaboration, K. S. Khaw et al., “Search for the muon electric dipole moment using frozen-spin technique at PSI,” PoS NuFact2021 (2022) 136, arXiv:2201.08729 [hep-ex].
  * [13] C. Cesarotti, Q. Lu, Y. Nakai, A. Parikh, and M. Reece, “Interpreting the Electron EDM Constraint,” JHEP 05 (2019) 059, arXiv:1810.07736 [hep-ph].
  * [14] Y. Nakai and M. Reece, “Electric Dipole Moments in Natural Supersymmetry,” JHEP 08 (2017) 031, arXiv:1612.08090 [hep-ph].
  * [15] T. Ibrahim and P. Nath, “The Neutron and the electron electric dipole moment in N=1 supergravity unification,” Phys. Rev. D 57 (1998) 478–488, arXiv:hep-ph/9708456. [Erratum: Phys.Rev.D 58, 019901 (1998), Erratum: Phys.Rev.D 60, 079903 (1999), Erratum: Phys.Rev.D 60, 119901 (1999)].
  * [16] T.-F. Feng, T. Huang, X.-Q. Li, X.-M. Zhang, and S.-M. Zhao, “Lepton dipole moments and rare decays in the CP violating MSSM with nonuniversal soft supersymmetry breaking,” Phys. Rev. D 68 (2003) 016004, arXiv:hep-ph/0305290.
  * [17] T.-F. Feng, X.-Q. Li, L. Lin, J. Maalampi, and H.-S. Song, “The Two-loop supersymmetric corrections to lepton anomalous magnetic and electric dipole moments,” Phys. Rev. D 73 (2006) 116001, arXiv:hep-ph/0604171.
  * [18] T.-F. Feng, L. Sun, and X.-Y. Yang, “Electroweak and supersymmetric two-loop corrections to lepton anomalous magnetic and electric dipole moments,” Nucl. Phys. B 800 (2008) 221–252, arXiv:0805.1122 [hep-ph].
  * [19] S.-M. Zhao, T.-F. Feng, X.-J. Zhan, H.-B. Zhang, and B. Yan, “The study of lepton EDM in CP violating BLMSSM,” JHEP 07 (2015) 124, arXiv:1411.4210 [hep-ph].
  * [20] L.-H. Su, D. He, X.-X. Dong, T.-F. Feng, and S.-M. Zhao, “Study of lepton EDMs in the U(1)X SSM *,” Chin. Phys. C 46 no. 9, (2022) 093103, arXiv:2201.00517 [hep-ph].
  * [21] Y. Nakai, R. Sato, and Y. Shigekami, “Muon electric dipole moment as a probe of flavor-diagonal CP violation,” Phys. Lett. B 831 (2022) 137194, arXiv:2204.03183 [hep-ph].
  * [22] G. Hiller, K. Huitu, T. Ruppell, and J. Laamanen, “A Large Muon Electric Dipole Moment from Flavor?,” Phys. Rev. D 82 (2010) 093015, arXiv:1008.5091 [hep-ph].
  * [23] Y. Omura, E. Senaha, and K. Tobe, “$\tau$\- and $\mu$-physics in a general two Higgs doublet model with $\mu-\tau$ flavor violation,” Phys. Rev. D 94 no. 5, (2016) 055019, arXiv:1511.08880 [hep-ph].
  * [24] Y. Abe, T. Toma, and K. Tsumura, “A $\mu$-$\tau$-philic scalar doublet under $Z_{n}$ flavor symmetry,” JHEP 06 (2019) 142, arXiv:1904.10908 [hep-ph].
  * [25] W.-S. Hou, G. Kumar, and S. Teunissen, “Charged lepton EDM with extra Yukawa couplings,” JHEP 01 (2022) 092, arXiv:2109.08936 [hep-ph].
  * [26] K.-m. Cheung, “Muon anomalous magnetic moment and leptoquark solutions,” Phys. Rev. D 64 (2001) 033001, arXiv:hep-ph/0102238.
  * [27] J. M. Arnold, B. Fornal, and M. B. Wise, “Phenomenology of scalar leptoquarks,” Phys. Rev. D 88 (2013) 035009, arXiv:1304.6119 [hep-ph].
  * [28] I. Doršner, S. Fajfer, A. Greljo, J. F. Kamenik, and N. Košnik, “Physics of leptoquarks in precision experiments and at particle colliders,” Phys. Rept. 641 (2016) 1–68, arXiv:1603.04993 [hep-ph].
  * [29] W. Dekens, J. de Vries, M. Jung, and K. K. Vos, “The phenomenology of electric dipole moments in models of scalar leptoquarks,” JHEP 01 (2019) 069, arXiv:1809.09114 [hep-ph].
  * [30] W. Altmannshofer, S. Gori, H. H. Patel, S. Profumo, and D. Tuckler, “Electric dipole moments in a leptoquark scenario for the $B$-physics anomalies,” JHEP 05 (2020) 069, arXiv:2002.01400 [hep-ph].
  * [31] K. S. Babu, P. S. B. Dev, S. Jana, and A. Thapa, “Unified framework for $B$-anomalies, muon $g-2$ and neutrino masses,” JHEP 03 (2021) 179, arXiv:2009.01771 [hep-ph].
  * [32] A. Crivellin and M. Hoferichter, “Consequences of chirally enhanced explanations of (g-2)μ for h $\to$ $\mu$$\mu$ and Z $\to$ $\mu$$\mu$,” JHEP 07 (2021) 135, arXiv:2104.03202 [hep-ph]. [Erratum: JHEP 10, 030 (2022)].
  * [33] A. Crivellin, M. Hoferichter, and P. Schmidt-Wellenburg, “Combined explanations of $(g-2)_{\mu,e}$ and implications for a large muon EDM,” Phys. Rev. D 98 no. 11, (2018) 113002, arXiv:1807.11484 [hep-ph].
  * [34] G. Hiller, C. Hormigos-Feliu, D. F. Litim, and T. Steudtner, “Model Building from Asymptotic Safety with Higgs and Flavor Portals,” Phys. Rev. D 102 no. 9, (2020) 095023, arXiv:2008.08606 [hep-ph].
  * [35] K. Hamaguchi, N. Nagata, G. Osaki, and S.-Y. Tseng, “Probing New Physics in the Vector-like Lepton Model by Lepton Electric Dipole Moments,” arXiv:2211.16800 [hep-ph].
  * [36] W. Yin, “Radiative lepton mass and muon g-2 with suppressed lepton flavor and CP violations,” JHEP 08 (2021) 043, arXiv:2103.14234 [hep-ph].
  * [37] M. J. Baker, P. Cox, and R. R. Volkas, “Radiative muon mass models and $(g-2)_{\mu}$,” JHEP 05 (2021) 174, arXiv:2103.13401 [hep-ph].
  * [38] W. Buchmuller, R. D. Peccei, and T. Yanagida, “Leptogenesis as the origin of matter,” Ann. Rev. Nucl. Part. Sci. 55 (2005) 311–355, arXiv:hep-ph/0502169.
  * [39] T. Hahn and M. Perez-Victoria, “Automatized one loop calculations in four-dimensions and D-dimensions,” Comput. Phys. Commun. 118 (1999) 153–165, arXiv:hep-ph/9807565.
  * [40] ATLAS Collaboration, G. Aad et al., “A search for the dimuon decay of the Standard Model Higgs boson with the ATLAS detector,” Phys. Lett. B 812 (2021) 135980, arXiv:2007.07830 [hep-ex].
  * [41] CMS Collaboration, A. M. Sirunyan et al., “Evidence for Higgs boson decay to a pair of muons,” JHEP 01 (2021) 148, arXiv:2009.04363 [hep-ex].
  * [42] LHC Higgs Cross Section Working Group Collaboration, D. de Florian et al., “Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector,” arXiv:1610.07922 [hep-ph].
  * [43] ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group, SLD Heavy Flavour Group Collaboration, S. Schael et al., “Precision electroweak measurements on the $Z$ resonance,” Phys. Rept. 427 (2006) 257–454, arXiv:hep-ex/0509008.
  * [44] M. J. Baker, P. Cox, and R. R. Volkas, “Has the Origin of the Third-Family Fermion Masses been Determined?,” JHEP 04 (2021) 151, arXiv:2012.10458 [hep-ph].
  * [45] K. Griest and D. Seckel, “Three exceptions in the calculation of relic abundances,” Phys. Rev. D 43 (1991) 3191–3203.
  * [46] P. Gondolo and G. Gelmini, “Cosmic abundances of stable particles: Improved analysis,” Nucl. Phys. B 360 (1991) 145–179.
  * [47] J. Edsjo and P. Gondolo, “Neutralino relic density including coannihilations,” Phys. Rev. D 56 (1997) 1879–1894, arXiv:hep-ph/9704361.
  * [48] G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, “micrOMEGAs5.0 : Freeze-in,” Comput. Phys. Commun. 231 (2018) 173–186, arXiv:1801.03509 [hep-ph].
  * [49] G. Belanger, A. Mjallal, and A. Pukhov, “Recasting direct detection limits within micrOMEGAs and implication for non-standard Dark Matter scenarios,” Eur. Phys. J. C 81 no. 3, (2021) 239, arXiv:2003.08621 [hep-ph].
  * [50] E. Del Nobile, “The Theory of Direct Dark Matter Detection: A Guide to Computations,” arXiv:2104.12785 [hep-ph].
  * [51] J. R. Ellis, A. Ferstl, and K. A. Olive, “Reevaluation of the elastic scattering of supersymmetric dark matter,” Phys. Lett. B 481 (2000) 304–314, arXiv:hep-ph/0001005.
  * [52] P. Gondolo, J. Edsjo, P. Ullio, L. Bergstrom, M. Schelke, and E. A. Baltz, “DarkSUSY: Computing supersymmetric dark matter properties numerically,” JCAP 07 (2004) 008, arXiv:astro-ph/0406204.
  * [53] J. R. Ellis, K. A. Olive, and C. Savage, “Hadronic Uncertainties in the Elastic Scattering of Supersymmetric Dark Matter,” Phys. Rev. D 77 (2008) 065026, arXiv:0801.3656 [hep-ph].
  * [54] G. Belanger, F. Boudjema, A. Pukhov, and A. Semenov, “Dark matter direct detection rate in a generic model with micrOMEGAs 2.2,” Comput. Phys. Commun. 180 (2009) 747–767, arXiv:0803.2360 [hep-ph].
  * [55] H.-Y. Cheng and C.-W. Chiang, “Revisiting Scalar and Pseudoscalar Couplings with Nucleons,” JHEP 07 (2012) 009, arXiv:1202.1292 [hep-ph].
  * [56] G. Belanger, F. Boudjema, A. Pukhov, and A. Semenov, “micrOMEGAs$\\_$3: A program for calculating dark matter observables,” Comput. Phys. Commun. 185 (2014) 960–985, arXiv:1305.0237 [hep-ph].
  * [57] A. Crivellin, M. Hoferichter, and M. Procura, “Accurate evaluation of hadronic uncertainties in spin-independent WIMP-nucleon scattering: Disentangling two- and three-flavor effects,” Phys. Rev. D 89 (2014) 054021, arXiv:1312.4951 [hep-ph].
  * [58] M. Hoferichter, J. Ruiz de Elvira, B. Kubis, and U.-G. Meißner, “High-Precision Determination of the Pion-Nucleon $\sigma$ Term from Roy-Steiner Equations,” Phys. Rev. Lett. 115 (2015) 092301, arXiv:1506.04142 [hep-ph].
  * [59] J. Ellis, N. Nagata, and K. A. Olive, “Uncertainties in WIMP Dark Matter Scattering Revisited,” Eur. Phys. J. C 78 no. 7, (2018) 569, arXiv:1805.09795 [hep-ph].
  * [60] R. H. Helm, “Inelastic and Elastic Scattering of 187-Mev Electrons from Selected Even-Even Nuclei,” Phys. Rev. 104 (1956) 1466–1475.
  * [61] J. D. Lewin and P. F. Smith, “Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil,” Astropart. Phys. 6 (1996) 87–112.
  * [62] LZ Collaboration, J. Aalbers et al., “First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment,” arXiv:2207.03764 [hep-ex].
  * [63] XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 no. 11, (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * [64] XENON Collaboration, E. Aprile et al., “Constraining the spin-dependent WIMP-nucleon cross sections with XENON1T,” Phys. Rev. Lett. 122 no. 14, (2019) 141301, arXiv:1902.03234 [astro-ph.CO].
  * [65] PandaX-4T Collaboration, Y. Meng et al., “Dark Matter Search Results from the PandaX-4T Commissioning Run,” Phys. Rev. Lett. 127 no. 26, (2021) 261802, arXiv:2107.13438 [hep-ex].
  * [66] PandaX Collaboration, J. Liu, “The first results of PandaX-4T,” Int. J. Mod. Phys. D 31 no. 05, (2022) 2230007.
  * [67] M. Pospelov and T. ter Veldhuis, “Direct and indirect limits on the electromagnetic form-factors of WIMPs,” Phys. Lett. B 480 (2000) 181–186, arXiv:hep-ph/0003010.
  * [68] C. M. Ho and R. J. Scherrer, “Anapole Dark Matter,” Phys. Lett. B 722 (2013) 341–346, arXiv:1211.0503 [hep-ph].
  * [69] B. W. Lee, C. Quigg, and H. B. Thacker, “The Strength of Weak Interactions at Very High-Energies and the Higgs Boson Mass,” Phys. Rev. Lett. 38 (1977) 883–885.
  * [70] B. W. Lee, C. Quigg, and H. B. Thacker, “Weak Interactions at Very High-Energies: The Role of the Higgs Boson Mass,” Phys. Rev. D 16 (1977) 1519.
  * [71] G. Cynolter, E. Lendvai, and G. Pocsik, “Note on unitarity constraints in a model for a singlet scalar dark matter candidate,” Acta Phys. Polon. B 36 (2005) 827–832, arXiv:hep-ph/0410102.
  * [72] S. K. Kang and J. Park, “Unitarity Constraints in the standard model with a singlet scalar field,” JHEP 04 (2015) 009, arXiv:1306.6713 [hep-ph].
  * [73] R. Costa, A. P. Morais, M. O. P. Sampaio, and R. Santos, “Two-loop stability of a complex singlet extended Standard Model,” Phys. Rev. D 92 (2015) 025024, arXiv:1411.4048 [hep-ph].
  * [74] R. Casalbuoni, D. Dominici, R. Gatto, and C. Giunti, “Strong Interacting Two Doublet and Doublet Singlet Higgs Models,” Phys. Lett. B 178 (1986) 235.
  * [75] R. Casalbuoni, D. Dominici, F. Feruglio, and R. Gatto, “Testing the Standard Model in Terms of a Possible Strong Scalar Sector,” Phys. Lett. B 200 (1988) 495–500.
  * [76] J. Maalampi, J. Sirkka, and I. Vilja, “Tree level unitarity and triviality bounds for two Higgs models,” Phys. Lett. B 265 (1991) 371–376.
  * [77] S. Kanemura, T. Kubota, and E. Takasugi, “Lee-Quigg-Thacker bounds for Higgs boson masses in a two doublet model,” Phys. Lett. B 313 (1993) 155–160, arXiv:hep-ph/9303263.
  * [78] I. F. Ginzburg and I. P. Ivanov, “Tree level unitarity constraints in the 2HDM with CP violation,” arXiv:hep-ph/0312374.
  * [79] A. G. Akeroyd, A. Arhrib, and E.-M. Naimi, “Note on tree level unitarity in the general two Higgs doublet model,” Phys. Lett. B 490 (2000) 119–124, arXiv:hep-ph/0006035.
  * [80] J. Horejsi and M. Kladiva, “Tree-unitarity bounds for THDM Higgs masses revisited,” Eur. Phys. J. C 46 (2006) 81–91, arXiv:hep-ph/0510154.
  * [81] W. Porod, “SPheno, a program for calculating supersymmetric spectra, SUSY particle decays and SUSY particle production at e+ e- colliders,” Comput. Phys. Commun. 153 (2003) 275–315, arXiv:hep-ph/0301101.
  * [82] W. Porod and F. Staub, “SPheno 3.1: Extensions including flavour, CP-phases and models beyond the MSSM,” Comput. Phys. Commun. 183 (2012) 2458–2469, arXiv:1104.1573 [hep-ph].
  * [83] F. Staub, “SARAH,” arXiv:0806.0538 [hep-ph].
  * [84] F. Staub, “From Superpotential to Model Files for FeynArts and CalcHep/CompHep,” Comput. Phys. Commun. 181 (2010) 1077–1086, arXiv:0909.2863 [hep-ph].
  * [85] F. Staub, “Automatic Calculation of supersymmetric Renormalization Group Equations and Self Energies,” Comput. Phys. Commun. 182 (2011) 808–833, arXiv:1002.0840 [hep-ph].
  * [86] F. Staub, “SARAH 3.2: Dirac Gauginos, UFO output, and more,” Comput. Phys. Commun. 184 (2013) 1792–1809, arXiv:1207.0906 [hep-ph].
  * [87] F. Staub, “SARAH 4 : A tool for (not only SUSY) model builders,” Comput. Phys. Commun. 185 (2014) 1773–1790, arXiv:1309.7223 [hep-ph].
  * [88] M. D. Goodsell and F. Staub, “Unitarity constraints on general scalar couplings with SARAH,” Eur. Phys. J. C 78 no. 8, (2018) 649, arXiv:1805.07306 [hep-ph].
  * [89] M. E. Peskin and T. Takeuchi, “A New constraint on a strongly interacting Higgs sector,” Phys. Rev. Lett. 65 (1990) 964–967.
  * [90] M. E. Peskin and T. Takeuchi, “Estimation of oblique electroweak corrections,” Phys. Rev. D 46 (1992) 381–409.
  * [91] W. Grimus, L. Lavoura, O. M. Ogreid, and P. Osland, “A Precision constraint on multi-Higgs-doublet models,” J. Phys. G 35 (2008) 075001, arXiv:0711.4022 [hep-ph].
  * [92] W. Grimus, L. Lavoura, O. M. Ogreid, and P. Osland, “The Oblique parameters in multi-Higgs-doublet models,” Nucl. Phys. B 801 (2008) 81–96, arXiv:0802.4353 [hep-ph].
  * [93] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [94] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [95] B. Pontecorvo, “Mesonium and anti-mesonium,” Sov. Phys. JETP 6 (1957) 429.
  * [96] G. Feinberg and S. Weinberg, “Conversion of Muonium into Antimuonium,” Phys. Rev. 123 (1961) 1439–1443.
  * [97] B. W. Lee and R. E. Shrock, “Natural Suppression of Symmetry Violation in Gauge Theories: Muon - Lepton and Electron Lepton Number Nonconservation,” Phys. Rev. D 16 (1977) 1444.
  * [98] B. W. Lee, S. Pakvasa, R. E. Shrock, and H. Sugawara, “Muon and Electron Number Nonconservation in a V-A Gauge Model,” Phys. Rev. Lett. 38 (1977) 937. [Erratum: Phys.Rev.Lett. 38, 1230 (1977)].
  * [99] L. Willmann et al., “New bounds from searching for muonium to anti-muonium conversion,” Phys. Rev. Lett. 82 (1999) 49–52, arXiv:hep-ex/9807011.
  * [100] J. de Blas et al., “Higgs Boson Studies at Future Particle Colliders,” JHEP 01 (2020) 139, arXiv:1905.03764 [hep-ph].
  * [101] ATLAS Collaboration, G. Aad et al., “Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions using the ATLAS detector,” Eur. Phys. J. C 80 no. 2, (2020) 123, arXiv:1908.08215 [hep-ex].
  * [102] ATLAS Collaboration, G. Aad et al., “Searches for electroweak production of supersymmetric particles with compressed mass spectra in $\sqrt{s}=$ 13 TeV $pp$ collisions with the ATLAS detector,” Phys. Rev. D 101 no. 5, (2020) 052005, arXiv:1911.12606 [hep-ex].
  * [103] CMS Collaboration, A. M. Sirunyan et al., “Search for supersymmetric partners of electrons and muons in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 790 (2019) 140–166, arXiv:1806.05264 [hep-ex].
  * [104] J. Kawamura, S. Okawa, and Y. Omura, “Current status and muon $g-2$ explanation of lepton portal dark matter,” JHEP 08 (2020) 042, arXiv:2002.12534 [hep-ph].
  * [105] M. Garny, A. Ibarra, M. Pato, and S. Vogl, “Internal bremsstrahlung signatures in light of direct dark matter searches,” JCAP 12 (2013) 046, arXiv:1306.6342 [hep-ph].
  * [106] PAMELA Collaboration, O. Adriani et al., “PAMELA results on the cosmic-ray antiproton flux from 60 MeV to 180 GeV in kinetic energy,” Phys. Rev. Lett. 105 (2010) 121101, arXiv:1007.0821 [astro-ph.HE].
  * [107] M. Garny, A. Ibarra, and S. Vogl, “Antiproton constraints on dark matter annihilations from internal electroweak bremsstrahlung,” JCAP 07 (2011) 028, arXiv:1105.5367 [hep-ph].
  * [108] M. Garny, A. Ibarra, and S. Vogl, “Dark matter annihilations into two light fermions and one gauge boson: General analysis and antiproton constraints,” JCAP 04 (2012) 033, arXiv:1112.5155 [hep-ph].
  * [109] B. J. Mount et al., “LUX-ZEPLIN (LZ) Technical Design Report,” arXiv:1703.09144 [physics.ins-det].
  * [110] LZ Collaboration, D. S. Akerib et al., “The LUX-ZEPLIN (LZ) Experiment,” Nucl. Instrum. Meth. A 953 (2020) 163047, arXiv:1910.09124 [physics.ins-det].
  * [111] PandaX Collaboration, H. Zhang et al., “Dark matter direct search sensitivity of the PandaX-4T experiment,” Sci. China Phys. Mech. Astron. 62 no. 3, (2019) 31011, arXiv:1806.02229 [physics.ins-det].
  * [112] H. Al Ali et al., “The muon Smasher’s guide,” Rept. Prog. Phys. 85 no. 8, (2022) 084201, arXiv:2103.14043 [hep-ph].
  * [113] T. Abe, R. Sato, and K. Yagyu, “Muon specific two-Higgs-doublet model,” JHEP 07 (2017) 012, arXiv:1705.01469 [hep-ph].
  * [114] M. Lindner, M. Platscher, and F. S. Queiroz, “A Call for New Physics : The Muon Anomalous Magnetic Moment and Lepton Flavor Violation,” Phys. Rept. 731 (2018) 1–82, arXiv:1610.06587 [hep-ph].
  * [115] MEG Collaboration, A. M. Baldini et al., “Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment,” Eur. Phys. J. C 76 no. 8, (2016) 434, arXiv:1605.05081 [hep-ex].
  * [116] R. Mertig, M. Bohm, and A. Denner, “FEYN CALC: Computer algebraic calculation of Feynman amplitudes,” Comput. Phys. Commun. 64 (1991) 345–359.
  * [117] V. Shtabovenko, R. Mertig, and F. Orellana, “New Developments in FeynCalc 9.0,” Comput. Phys. Commun. 207 (2016) 432–444, arXiv:1601.01167 [hep-ph].
  * [118] V. Shtabovenko, R. Mertig, and F. Orellana, “FeynCalc 9.3: New features and improvements,” Comput. Phys. Commun. 256 (2020) 107478, arXiv:2001.04407 [hep-ph].
  * [119] G. ’t Hooft and M. J. G. Veltman, “Regularization and Renormalization of Gauge Fields,” Nucl. Phys. B 44 (1972) 189–213.
  * [120] P. Breitenlohner and D. Maison, “Dimensional Renormalization and the Action Principle,” Commun. Math. Phys. 52 (1977) 11–38.