compat=1.1.1 capbtabboxtable[][]

# Probing CP Violation in Dark Sector through
the Electron Electric Dipole Moment

Jia Liu<EMAIL_ADDRESS>School of Physics and State Key Laboratory of
Nuclear Physics and Technology,
Peking University, Beijing 100871, China  Center for High Energy Physics,
Peking University, Beijing 100871, China  Yuichiro Nakai<EMAIL_ADDRESS>Tsung-Dao Lee Institute, Shanghai Jiao Tong University,
520 Shengrong Road, Shanghai 201210, China  School of Physics and Astronomy,
Shanghai Jiao Tong University,
800 Dongchuan Road, Shanghai 200240, China
Yoshihiro Shigekami<EMAIL_ADDRESS>Tsung-Dao Lee Institute, Shanghai
Jiao Tong University,
520 Shengrong Road, Shanghai 201210, China  School of Physics and Astronomy,
Shanghai Jiao Tong University,
800 Dongchuan Road, Shanghai 200240, China  Muyuan Song<EMAIL_ADDRESS>Center for High Energy Physics, Peking University, Beijing 100871, China
School of Physics and State Key Laboratory of Nuclear Physics and Technology,
Peking University, Beijing 100871, China

###### Abstract

The Two Higgs Doublet Model (2HDM) stands as a promising framework for
exploring physics beyond the Standard Model (SM). Within this context, we
explore the possibility that the two Higgs doublets may serve as a window into
CP-violating dark sectors, neutral under the SM gauge groups. Specifically,
our focus is on investigating the electric dipole moment (EDM) of the
electron, generated solely by CP violation in the dark sector. We present a
general formula for the electron EDM, without specifying the structure of the
dark sectors, and discuss the current constraints on various dark sector
models. It is noteworthy that even in the case of a CP-conserving 2HDM, the
resulting electron EDM is capable of reaching the current experimental limit,
with CP violation arising exclusively from the dark sectors. Furthermore, we
introduce a heavy dark sector (HDS) approximation for the analytic calculation
of the EDM, assuming that the dark sector particles are much heavier than the
physical states in the 2HDM. This approximation yields simplified analytic
results that are consistent with the full numerical calculations.

###### Contents

  1. 1 Introduction
  2. 2 The framework
    1. 2.1 Two Higgs Doublet Model
    2. 2.2 2HDM portal to the dark sector
  3. 3 EDMs from general dark sector
    1. 3.1 Current limits on EDMs
    2. 3.2 The CP-violating diagrams
    3. 3.3 The four-fermion diagram
    4. 3.4 The Barr-Zee diagrams
      1. 3.4.1 The heavy dark sector limit
  4. 4 Numerical analysis for benchmark models
    1. 4.1 Scalar dark sector model
      1. 4.1.1 Numerical results
      2. 4.1.2 Discussion on the HDS approximation
    2. 4.2 Fermion dark sector model
      1. 4.2.1 Numerical results
      2. 4.2.2 Discussion on the HDS approximation
  5. 5 Conclusion
  6. A The 2HDM
  7. B Potential analysis for the fermion dark sector model
  8. C Rephasing-invariant CP-violating phases
    1. C.1 CPV phases in the scalar dark sector model
    2. C.2 CPV phases in the fermion dark sector model
  9. D Details of calculations
    1. D.1 General dark sector
    2. D.2 Heavy dark sector
    3. D.3 Decomposed expressions for the electron EDM

##  1 Introduction

Identifying the structure of the Higgs sector to spontaneously break the
electroweak symmetry is one of the most important missions in modern particle
physics. Despite that the discovery of the Higgs boson at the Large Hadron
Collider (LHC) [1, 2] largely advances our understanding of the Higgs sector,
the current experimental data still allow a lot of freedom for its structure.
Among others, two Higgs doublet models (2HDMs) [3] are classified into the
simplest and viable extensions of the Standard Model (SM) Higgs sector (see
$e.g.$ ref. [4] for a review). 2HDMs are motivated by various phenomenological
reasons. One of them is supersymmetry, which requires two Higgs doublets to
give masses to both up-type and down-type fermions [5, 6, 7, 8]. Another
motivation is dark matter (DM), which can be provided by a stable neutral
scalar in some versions of 2HDMs, such as the inert doublet model [9, 10, 11,
12, 13, 14]. Therefore, the paradigm of 2HDMs has been considered one of the
top candidates of physics beyond SM. In the present paper, we point out a new
interesting feature of 2HDMs: two Higgs doublets become a window with a good
view of CP-violating dark sectors that are neutral under the SM gauge groups
but interact with SM particles through two Higgs doublet portal couplings.

Exploration of dark sector models has received a lot of attention recently. In
particular, if a dark sector contains a new source of CP violation, it may be
able to address the matter-antimatter asymmetry of the Universe. In fact,
refs. [15, 16] have demonstrated a connection between CP violation in a dark
sector and electroweak baryogenesis, which is responsible for the observed
baryon asymmetry during the electroweak phase transition (EWPT). While new
light particles in a dark sector are often probed by high-intensity
accelerator-based experiments, a CP-violating dark sector may generate
electric dipole moments (EDMs) of nucleons, atoms, and molecules that are
precisely measured. Among others, the measurement of the electron EDM is
fascinating. The recent results of the ACME [17] and JILA [18, 19] experiments
have reached $\mathcal{O}(10^{-29}-10^{-30})\,e\rm cm$. Since the electron EDM
is immune to QCD effects and its precise calculation is possible, it seems to
provide a sensitive probe of a CP-violating dark sector. The authors of ref.
[20] have initiated to study of this issue by considering various dark sectors
feebly interacting with the SM sector. They found that the effect of dark
sector CP violation appears only from higher loops and is too tiny to be
observed unless a neutrino portal coupling contains a nonzero complex phase.
Ref. [21] has extended the Higgs sector by a complex singlet coupling to a
dark sector fermion. The authors showed that such a complex singlet extension
does not generate the electron EDM up to at least the two-loop level.

In the present paper, we explore the effect of dark sector CP violation within
the framework of 2HDMs and calculate molecule EDMs that constrain the electron
EDM.111 Ref. [22] has discussed a similar setup in a different context and
studied a dark sector contribution to anomalous gauge couplings. A new CP-
violating phase, which can be $\mathcal{O}(1)$, exists in an interaction
involving dark sector particles, while the ordinary 2HDM sector does not
contain CP violation. We find that the produced electron EDM reaches the
current experimental limit in some parameter space. The 2HDM portal enhances
the sensitivity of the electron EDM to a CP-violating dark sector.

The rest of the paper is organized as follows. In section 2, we introduce our
setup. Section 3 then presents a general formula for the electron EDM without
specifying the structure of a dark sector. To discuss the current constraints
numerically, we specify the dark sector structure and consider two cases of a
dark complex scalar and a dark vector-like fermion in section 4. Section 5 is
devoted to conclusions and discussions. Some calculational details are
summarized in appendices.

##  2 The framework

Let us start with an overview of the features of our CP-conserving 2HDMs and
CP-violating dark sectors under consideration in the present paper. The 2HDM
is a well-studied case of physics beyond the SM, involving an extra Higgs
doublet scalar field. Here, our objective is to provide a brief description of
the relevant physical modes and interactions. More details are summarized in
appendix A. While the current section does not specify a concrete dark sector
model, dark sector CP violation will be encoded through a (radiatively
generated) mixing term between CP-even and CP-odd scalars of the 2HDM.

###  2.1 Two Higgs Doublet Model

The potential of two Higgs doublets $H_{1,2}$ can be expressed as

$\displaystyle V_{H}$
$\displaystyle=\sum_{i=1,2}\Bigl{[}m_{i}^{2}H_{i}^{\dagger}H_{i}+\frac{\lambda_{i}}{2}\left(H_{i}^{\dagger}H_{i}\right)^{2}\Bigr{]}+\lambda_{3}\left(H_{1}^{\dagger}H_{1}\right)\left(H_{2}^{\dagger}H_{2}\right)$
$\displaystyle\hskip
11.99998pt+\lambda_{4}\left(H_{1}^{\dagger}H_{2}\right)\left(H_{2}^{\dagger}H_{1}\right)+\left[-m_{12}^{2}H_{1}^{\dagger}H_{2}+\frac{\lambda_{5}}{2}\left(H_{1}^{\dagger}H_{2}\right)^{2}+{\rm
h.c.}\right]\,,$ ( 2.1)

where $m_{12}^{2}$ and $\lambda_{5}$ are generally complex, while the other
parameters are real. By field redefinition of either $H_{1}$ or $H_{2}$, we
can set $m_{12}^{2}$ or $\lambda_{5}$ to be real, while
$\arg[\lambda_{5}^{*}(m_{12}^{2})^{2}]$ remains invariant. This phase is
severely constrained by the electron EDM measurement, as previously studied in
refs. [23, 24, 25, 26, 27, 28, 29, 30, 31]. Since our focus is on dark sector
CP violation and how it is transferred into the 2HDM sector, we assume that
the pure 2HDM sector is CP-conserving,
$\arg[\lambda_{5}^{*}(m_{12}^{2})^{2}]=0$, and both $\lambda_{5}$ and
$m_{12}^{2}$ are real. The two $SU(2)_{L}$ doublets $H_{1,2}$ are parametrized
as

$\displaystyle H_{i}=\begin{pmatrix}H_{i}^{+}\\\
\left(v_{i}+H_{i}^{0}+iA_{i}^{0}\right)/\sqrt{2}\end{pmatrix}\leavevmode\nobreak\
\leavevmode\nobreak\ (i=1,2)\,.$ ( 2.2)

Here, $v_{1}$ and $v_{2}$ are nonzero vacuum expectation values (VEVs) for the
neutral components of $H_{1}$ and $H_{2}$, respectively. They are related by
$v_{1}^{2}+v_{2}^{2}=v_{\rm SM}^{2}=(246.22\,\rm GeV)^{2}$. We define their
ratio as

$\displaystyle\tan\beta\equiv\frac{v_{2}}{v_{1}}\,.$ ( 2.3)

Upon diagonalizing the mass matrices derived from the scalar potential in Eq.
(2.1), we obtain the mass eigenstates for the neutral scalars:

$\displaystyle H_{1}^{0}$ $\displaystyle=H^{0}c_{\alpha}-h^{0}s_{\alpha}\,,$
$\displaystyle H_{2}^{0}=H^{0}s_{\alpha}+h^{0}c_{\alpha}\,,$ ( 2.4)
$\displaystyle A_{1}^{0}$ $\displaystyle=G^{0}c_{\beta}-A^{0}s_{\beta}\,,$
$\displaystyle A_{2}^{0}=G^{0}s_{\beta}+A^{0}c_{\beta}\,,$ ( 2.5)

with mixing angle $\alpha$ whose expression is provided in Eq. (A.11). Here,
$G^{0}$ represents the Nambu-Goldstone (NG) field, while $h^{0}$ and $H^{0}$
denote the two CP-even scalars, with the former being the SM Higgs, and
$A^{0}$ represents the CP-odd heavy scalar.

Type | Yukawa terms
---|---
Type-I ($d_{R}:+$, $e_{R}:+$) | $-Y_{u}^{(2)}\overline{Q_{L}}\widetilde{H_{2}}u_{R}-Y_{d}^{(2)}\overline{Q_{L}}H_{2}d_{R}-Y_{e}^{(2)}\overline{L_{L}}H_{2}e_{R}+{\rm h.c.}$
Type-II ($d_{R}:-$, $e_{R}:-$) | $-Y_{u}^{(2)}\overline{Q_{L}}\widetilde{H_{2}}u_{R}-Y_{d}^{(1)}\overline{Q_{L}}H_{1}d_{R}-Y_{e}^{(1)}\overline{L_{L}}H_{1}e_{R}+{\rm h.c.}$
Type-X ($d_{R}:+$, $e_{R}:-$) | $-Y_{u}^{(2)}\overline{Q_{L}}\widetilde{H_{2}}u_{R}-Y_{d}^{(2)}\overline{Q_{L}}H_{2}d_{R}-Y_{e}^{(1)}\overline{L_{L}}H_{1}e_{R}+{\rm h.c.}$
Type-Y ($d_{R}:-$, $e_{R}:+$) | $-Y_{u}^{(2)}\overline{Q_{L}}\widetilde{H_{2}}u_{R}-Y_{d}^{(1)}\overline{Q_{L}}H_{1}d_{R}-Y_{e}^{(2)}\overline{L_{L}}H_{2}e_{R}+{\rm h.c.}$
Table 1: The four types of Yukawa couplings in the 2HDM, showing the $Z_{2}$
charges of $d_{R}$ and $e_{R}$ along with the type names in the first column
and the allowed Yukawa couplings in the second column. Note that $Q_{L}$ and
$L_{L}$ are assumed to be even under the $Z_{2}$.

The Yukawa couplings for SM fermions in the 2HDM are represented by the
following Lagrangian term:

$\displaystyle\mathcal{L}_{\rm
Yuk}=-Y_{u}^{(i)}\overline{Q_{L}}\widetilde{H}_{i}u_{R}-Y_{d}^{(i)}\overline{Q_{L}}H_{i}d_{R}-Y_{e}^{(i)}\overline{L_{L}}H_{i}e_{R}+{\rm
h.c.}\,,$ ( 2.6)

where $\widetilde{H}_{i}=i\sigma_{2}H_{i}^{*}$, $Q_{L}$ and $L_{L}$ denote the
left-handed $SU(2)_{L}$ doublet quark and lepton, and $u_{R}$, $d_{R}$ and
$e_{R}$ are the right-handed $SU(2)_{L}$ singlet up-type quark, down-type
quark and charged lepton. In general, there are Yukawa couplings associated
with both $H_{1}$ and $H_{2}$ for each SM fermion. However, including both
Yukawa couplings easily leads to dangerous flavor-changing neutral current
(FCNC) processes. To avoid such a situation, one can assume different symmetry
properties for $H_{1}$ and $H_{2}$ [32]. For example, $H_{1}$ behaves as odd
while $H_{2}$ as even under a discrete $Z_{2}$ symmetry.222 Under the $Z_{2}$
symmetry, the scalar potential term $m_{12}^{2}H_{1}^{\dagger}H_{2}$ in Eq.
(2.1) can be considered as a soft $Z_{2}$ breaking term. Without loss of
generality, one can select a scenario in which $\overline{Q_{L}}u_{R}$ couples
exclusively to $H_{2}$. Then, there are four types of Yukawa coupling
arrangements, namely Type-I, II, X, and Y, depending on the charges of $d_{R}$
and $e_{R}$, which are summarized in Table 1. Due to the mixings in Eqs. (2.4)
and (2.5), the Yukawa interactions between the SM fermions and the physical
neutral scalars are given by

$\displaystyle\mathcal{L}_{\rm Yuk}\supset-\sum_{f=u,d,e}\frac{m_{f}}{v_{\rm
SM}}\left(\xi_{h}^{f}\bar{f}fh^{0}+\xi_{H}^{f}\bar{f}fH^{0}-\xi_{A}^{f}\bar{f}i\gamma^{5}fA^{0}\right)\,,$
( 2.7)

where $m_{f}$ is the SM fermion mass, and the coefficients $\xi_{h,H,A}^{f}$
are summarized in Table 2. It is notable that the lepton couplings to the
heavy scalars are significantly enhanced by a large $\tan\beta$ for Type-II
and Type-X 2HDMs, leading to a larger effect on the electron EDM.

| $\xi_{h}^{u}$ | $\xi_{H}^{u}$ | $\xi_{A}^{u}$ | $\xi_{h}^{d}$ | $\xi_{H}^{d}$ | $\xi_{A}^{d}$ | $\xi_{h}^{e}$ | $\xi_{H}^{e}$ | $\xi_{A}^{e}$
---|---|---|---|---|---|---|---|---|---
Type-I | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $\cot\beta$ | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $-\cot\beta$ | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $-\cot\beta$
Type-II | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $\cot\beta$ | $-s_{\alpha}/c_{\beta}$ | $c_{\alpha}/c_{\beta}$ | $\tan\beta$ | $-s_{\alpha}/c_{\beta}$ | $c_{\alpha}/c_{\beta}$ | $\tan\beta$
Type-X | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $\cot\beta$ | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $-\cot\beta$ | $-s_{\alpha}/c_{\beta}$ | $c_{\alpha}/c_{\beta}$ | $\tan\beta$
Type-Y | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $\cot\beta$ | $-s_{\alpha}/c_{\beta}$ | $c_{\alpha}/c_{\beta}$ | $\tan\beta$ | $c_{\alpha}/s_{\beta}$ | $s_{\alpha}/s_{\beta}$ | $-\cot\beta$
Table 2: The coefficients between SM fermions and physical neutral scalars.

Due to the additional Higgs doublet, the $W$ boson couplings to CP-even
scalars are modified as

$\displaystyle\mathcal{L}\supset\frac{2m_{W}^{2}}{v_{\rm
SM}}s_{\beta-\alpha}W_{\mu}^{+}W^{-\,\mu}h^{0}+\frac{2m_{W}^{2}}{v_{\rm
SM}}c_{\beta-\alpha}W_{\mu}^{+}W^{-\,\mu}H^{0}\,.$ ( 2.8)

Note that if the model does not mix between CP-even and CP-odd scalars, there
is no $W$ boson coupling to the CP-odd scalar, $A^{0}$ [33].

###  2.2 2HDM portal to the dark sector

To couple the 2HDM to a dark sector, which is neutral under the SM gauge
symmetry, we make the following assumptions:

  * (i)

The dark sector incorporates physical CP phase(s), and it does not contain any
new scalar field that acquires a nonzero VEV.

  * (ii)

New scalar potential terms that couple $H_{1,2}$ to the dark sector are
expressed as

$\displaystyle
V_{H\mathcal{O}}=\lambda_{ij}H_{i}^{\dagger}H_{j}\mathcal{O}_{ij}^{D}+{\rm
h.c.}\,,$ ( 2.9)

where $\mathcal{O}^{D}_{ij}$ represents operators composed of dark sector
fields.

The assumption (i) implies the absence of tree-level mass mixing between CP-
even and CP-odd scalars in $H_{1,2}$. Loop corrections through the
interactions (2.9) generate effective mixing terms for those scalars:

$\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman[large] \vertex(a1)
{${h^{0},H^{0}}$}; \vertex[blob] [right=2.8cm of a1] (b1) {};
\vertex[right=2.5cm of b1] (c1) {$A^{0}$}; \diagram[medium] { (a1) -- [scalar,
momentum=$k$] (b1) -- [scalar] (c1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}},$ ( 2.10)

where the gray shaded blob represents loops of dark sector particles. Nonzero
mass mixing between CP-even ($h^{0},\leavevmode\nobreak\ H^{0}$) and CP-odd
($A^{0}$) scalars in the 2HDM is induced by dark sector CP violation.

##  3 EDMs from general dark sector

We now summarize the current experimental limits on molecule EDMs which
constrain the electron EDM, and derive analytical expressions for dark sector
contributions to the EDMs without specifying the detailed structure of the
dark sector. As mentioned in the previous section, the mixing between CP-even
and CP-odd scalars is induced via loops of CP-violating dark sector particles.

###  3.1 Current limits on EDMs

Measuring EDMs of leptons, nucleons, atoms, and molecules stands as one of the
most sensitive probes of new physics with CP violation, and numerous
experiments set upper limits on their magnitudes [17, 18, 34, 35], which put
stringent constraints on models incorporating CP violation, such as the
complex 2HDM and supersymmetry (see e.g. refs. [36, 37]). For an in-depth
understanding of EDMs and their measurement techniques, along with theoretical
aspects, one can refer to comprehensive reviews such as refs. [38, 39, 40].

One intriguing approach in EDM experiments is to utilize paramagnetic atoms or
molecules. A particularly stringent constraint on the electron EDM from this
type of experiment was reported by the ACME collaboration before 2022 [17].
They conducted a measurement on the EDM of the paramagnetic ThO molecule,
establishing an upper bound of $|d_{\rm ThO}|<1.1\times 10^{-29}\,e\,{\rm cm}$
at $90\%$ confidence level (CL). For paramagnetic molecules, the EDM
measurement offers sensitivity not only to the electron EDM $d_{e}$ but also
to the CP-odd electron-nucleon interaction $C_{S}$. The relevant dimension-
five and dimension-six operators are described by the following Lagrangian:

$\displaystyle\mathcal{L}_{{\rm
EDM},eN}\supset-\frac{i}{2}d_{e}\bar{e}\sigma^{\mu\nu}\gamma^{5}eF_{\mu\nu}-C_{S}\bar{e}i\gamma^{5}e\bar{N}N\,,$
( 3.1)

where $\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$, and $N$
represents the nucleon. These CP-violating interactions contribute to the EDM
of the paramagnetic molecule, and the bound from the ACME collaboration is
expressed as [17]

$\displaystyle\left|d_{e}+k_{\rm ThO}C_{S}\right|<1.1\times 10^{-29}\,e\,{\rm
cm}\,,$ ( 3.2)

with $k_{\rm ThO}\simeq 1.8\times 10^{-15}\,{\rm GeV}^{2}e\,{\rm cm}$, where
we have used $\mathcal{E}_{\rm eff}^{\rm ThO}\approx 78\,{\rm GVcm}^{-1}$ for
the effective electric field and $W_{S}^{\rm ThO}=-282\,{\rm kHz}\times h$ for
the molecule specific constant.

Towards the end of 2022, the JILA experiments presented a new upper bound on
the electron EDM, utilizing the HfF+ molecule ion [18, 19]. The result is
expressed as $|d_{{\rm HfF}^{+}}|<4.1\times 10^{-30}\,e\,{\rm cm}$ at 90% CL,
which can be interpreted as

$\displaystyle\left|d_{e}+k_{{\rm HfF}^{+}}C_{S}\right|<4.1\times
10^{-30}\,e\,{\rm cm}\,,$ ( 3.3)

with $k_{{\rm HfF}^{+}}\simeq 1.1\times 10^{-15}\,{\rm GeV}^{2}e\,{\rm cm}$,
where we have used $\mathcal{E}_{\rm eff}^{{\rm HfF}^{+}}\approx 23\,{\rm
GVcm}^{-1}$ for the effective electric field and $W_{S}^{{\rm
HfF}^{+}}=-51\,{\rm kHz}\times h$ for the molecule specific constant.

In addition to EDMs of paramagnetic molecules, there are experiments to
measure EDMs of diamagnetic atoms, such as the mercury EDM [41] and the
neutron EDM [35]. These experiments, in principle, can impose constraints on
the 2HDM or supersymmetry, and the details can be found in refs. [42, 43].
However, upon careful examination, we have verified that the electron EDM
constraint in the JILA experiment is considerably stronger than the bounds
stemming from the other EDM experiments. Consequently, we focus on the
constraint obtained by the JILA experiment in the following discussion.

###  3.2 The CP-violating diagrams

Let us now consider dark sector contributions to the CP-violating Lagrangian
in Eq. (3.1). The relevant contributions are shown in Fig. 1. The left diagram
leads to the four-fermion operator $C_{S}\bar{e}i\gamma^{5}e\bar{N}N$, while
the right corresponds to the Barr-Zee (BZ) diagram that results in the
electron EDM operator $d_{e}\bar{e}\sigma^{\mu\nu}\gamma^{5}eF_{\mu\nu}$ [44].
We can also consider one-loop diagrams induced by loops of neutral scalars,
but when compared with the diagrams $(a)$, $(b)$ in Fig. 1, these one-loop
diagrams are suppressed by the cubic power of the electron mass $m_{e}^{3}$:
one arises from a chirality flip on the electron line and the others stem from
electron-neutral scalar couplings. Therefore, for the current study, we can
neglect these one-loop diagrams.

$f$$f$$f^{\prime}$$f^{\prime}$$(a)$$\phi_{j}$$\eta_{a}$$e$$e$$\gamma$$(b)$$G$${f,W}$
Figure 1: The two diagrams relevant to the CP-violating Lagrangian in Eq.
(3.1). Left: The four-fermion operators mediated by the scalars. Right: The
Barr-Zee diagram contributing to the electron EDM, with $G$ denoting the
photon or $Z$ boson. In both diagrams, the scalar lines represent the CP-even
scalars $\phi_{j}$ or CP-odd scalars $\eta_{a}$, and the gray shaded blob
symbolizes the effective mixing between $\phi_{j}$ and $\eta_{a}$ arising from
the CP-violating dark sector contribution, which corresponds to $F_{\rm
DS}^{\phi_{j}\eta_{a}}(k^{2},M^{2})$ in Eq. (3.6).

We keep our calculation as general as possible by incorporating a general CP-
violating dark sector contribution into the two diagrams of Fig. 1 as
propagator corrections. Moreover, we consider a general scalar model with
$N_{E}$ CP-even and $N_{O}$ CP-odd neutral scalars contributing to the
diagrams, denoted as $\phi_{j}$ and $\eta_{a}$, respectively. For example, the
2HDM can be identified with $(N_{E},N_{O})=(2,2)$ and
$(\phi_{1},\phi_{2},\eta_{1},\eta_{2})=(h^{0},H^{0},G^{0},A^{0})$. Note that
we assume $\eta_{1}$ always represents the NG field $G^{0}$, and thus, the sum
of the index $a$ should start from $2$ for $\eta_{a}$ scalars. $\phi_{j}$ and
$\eta_{a}$ are considered as mass eigenstates, and hence, masses for these
scalars are described as

$\displaystyle V_{\rm
scl}\supset\sum_{j=1}^{N_{E}}m_{\phi_{j}}^{2}\phi_{j}^{2}+\sum_{a=2}^{N_{O}}m_{\eta_{a}}^{2}\eta_{a}^{2}\,.$
( 3.4)

Here, we assume
$m_{\phi_{1}}^{2}<m_{\phi_{2}}^{2}<\cdots<m_{\phi_{N_{E}}}^{2}$ and
$m_{\eta_{2}}^{2}<m_{\eta_{3}}^{2}<\cdots<m_{\eta_{N_{O}}}^{2}$.

The neutral scalars $\phi_{j},\eta_{a}$ couple to the SM fermions through the
Yukawa terms in the 2HDM. In addition, CP-even scalars also couple to the $W$
boson. We then have the following Lagrangian terms:

$\displaystyle\mathcal{L}$
$\displaystyle\supset-\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\sum_{f=u,d,e}\frac{m_{f}}{v_{\rm
SM}}\left(\xi_{\phi_{j}}^{f}\bar{f}f\phi_{j}-\xi_{\eta_{a}}^{f}\bar{f}i\gamma^{5}f\eta_{a}\right)+\sum_{j=1}^{N_{E}}\frac{2m_{W}^{2}}{v_{\rm
SM}}\xi_{\phi_{j}}^{W}W_{\mu}^{+}W^{-\,\mu}\phi_{j}\,,$ ( 3.5)

where $\xi_{\phi_{j},\eta_{a}}^{f}$ and $\xi_{\phi_{j}}^{W}$ are related to
$\tan\beta$ and mixing angles for $\phi_{j}$ and $\eta_{a}$ in the visible
sector (for the 2HDM, see Table 2 and Eq. (2.8)). Note that there is no
coupling between a CP-odd scalar and the $W$ boson because our assumption
forbids any tree-level mixing between CP-even and CP-odd scalars in the
visible sector. Furthermore, since we do not extend the fermion and gauge
sectors, couplings among these sectors are unchanged from the SM ones.
Finally, we have the effective mixing between $\phi_{j}$ and $\eta_{a}$, which
is generated in the same manner as in Eq. (2.10). In our analytical
calculation, we define its mixing as the following term in the effective
Lagrangian:

$\displaystyle\mathcal{L}_{\rm eff}\supset F_{\rm
DS}^{\phi_{j}\eta_{a}}(k^{2},M^{2})\phi_{j}\eta_{a}\,,$ ( 3.6)

where $k$ is four-momentum passing through this effective mixing, and $M$ is a
typical mass scale of a dark sector. The gray shaded blobs in Fig. 1
correspond to $F_{\rm DS}^{\phi_{j}\eta_{a}}(k^{2},M^{2})$.

###  3.3 The four-fermion diagram

In contrast to the nucleon operator in Eq. (3.1), we introduce a new four-
fermion operator at the quark level,

$\displaystyle\mathcal{L}_{4\mathchar 45\relax{\rm
fermi}}=C_{ff^{\prime}}\left(\bar{f}f\right)\left(\bar{f^{\prime}}i\gamma^{5}f^{\prime}\right)\,.$
( 3.7)

The coefficient $C_{ff^{\prime}}$ stemming from the diagram $(a)$ can be
readily calculated. Moreover, we can set the momentum of the scalar to be zero
for the current purpose, i.e., $F_{\rm
DS}^{\phi_{j}\eta_{a}}(k^{2},M^{2})=F_{\rm DS}^{\phi_{j}\eta_{a}}(0,M^{2})$.
Consequently, we obtain

$\displaystyle C_{ff^{\prime}}=\frac{m_{f}m_{f^{\prime}}}{v_{\rm
SM}^{2}}\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\frac{\xi_{\phi_{j}}^{f}}{m_{\phi_{j}}^{2}}\frac{\xi_{\eta_{a}}^{f^{\prime}}}{m_{\eta_{a}}^{2}}F_{\rm
DS}^{\phi_{j}\eta_{a}}(0,M^{2})\,,$ ( 3.8)

where $\xi_{\phi_{j}}^{f}$ and $\xi_{\eta_{a}}^{f^{\prime}}$ are defined in
Eq. (3.5).

The $C_{ff^{\prime}}$ term contributes to $C_{S}$ in Eq. (3.1). The
coefficient of the electron-quark interaction $C_{qe}$ leads to

$\displaystyle C_{S}\approx\sum_{q}C_{qe}\langle N|\bar{q}q|N\rangle$
$\displaystyle\approx C_{ue}\frac{16\,{\rm MeV}}{m_{u}}+C_{de}\frac{29\,{\rm
MeV}}{m_{d}}+C_{se}\frac{49\,{\rm MeV}}{m_{s}}$ $\displaystyle\hskip
11.99998pt+C_{ce}\frac{76\,{\rm MeV}}{m_{c}}+C_{be}\frac{74\,{\rm
MeV}}{m_{b}}+C_{te}\frac{77\,{\rm MeV}}{m_{t}}\,,$ ( 3.9)

where numerators arise from the matrix elements given in refs. [42, 45],

$\displaystyle(m_{u}+m_{d})\langle N|\bar{u}u+\bar{d}d|N\rangle\simeq 90\,{\rm
MeV}\,,\quad\langle N|\bar{u}u-\bar{d}d|N\rangle\simeq 0\,,\quad m_{s}\langle
N|\bar{s}s|N\rangle\simeq 49\,{\rm MeV}\,,$ $\displaystyle m_{c}\langle
N|\bar{c}c|N\rangle\simeq 76\,{\rm MeV}\,,\quad m_{b}\langle
N|\bar{b}b|N\rangle\simeq 74\,{\rm MeV}\,,\quad m_{t}\langle
N|\bar{t}t|N\rangle\simeq 77\,{\rm MeV}\,,$ ( 3.10)

with $m_{u}/m_{d}\simeq 0.55$.

###  3.4 The Barr-Zee diagrams

To calculate the BZ diagram in Fig. 1, it is important to handle the momentum
in $F_{\rm DS}^{\phi_{j}\eta_{a}}$ with care, as it is integrated as the outer
loop momentum in the BZ diagram. Additionally, since we have not specified our
dark sector at this stage, we cannot have any explicit formula for $F_{\rm
DS}^{\phi_{j}\eta_{a}}$. Therefore, we present the general result before the
loop momentum integration. Before the calculation, we define the relevant
couplings for SM fermions as

$\displaystyle\mathcal{L}\supset\bar{f}\left(g_{f\Phi}^{S}+ig_{f\Phi}^{P}\gamma^{5}\right)f\Phi+\bar{f}\left(g_{fG}^{V}\gamma^{\mu}+g_{fG}^{A}\gamma^{\mu}\gamma^{5}\right)fG_{\mu}\,,$
( 3.11)

where $\Phi$ and $G_{\mu}$ denote the neutral scalars ($\phi_{j},\eta_{a}$)
and neutral gauge bosons ($Z_{\mu},A_{\mu}$), respectively. Comparing with Eq.
(3.5), $g_{f\Phi}^{S,P}$ are related to $\xi_{\phi_{j},\eta_{a}}^{f}$, and
$g_{fG}^{V,A}$ are the usual gauge boson couplings in the SM, as summarized in
Table 3.

| $\Phi=\phi_{j}$ | $\Phi=\eta_{a}$ |  | $G=Z$ | $G=\gamma$
---|---|---|---|---|---
$g_{f\Phi}^{S}$ | ${\displaystyle-\frac{m_{f}}{v_{\rm SM}}\xi_{\phi_{j}}^{f}}$ | $0$ | $g_{eG}^{V}$ | ${\displaystyle\frac{g}{\cos\theta_{W}}\left(\frac{1}{2}T_{3}^{f}-Q_{f}\sin^{2}\theta_{W}\right)}$ | ${\rm e}Q_{f}$
$g_{f\Phi}^{P}$ | $0$ | ${\displaystyle+\frac{m_{f}}{v_{\rm SM}}\xi_{\eta_{a}}^{f}}$ | $g_{eG}^{A}$ | ${\displaystyle-\frac{g}{2\cos\theta_{W}}T_{3}^{f}}$ | $0$
Table 3: Couplings defined in Eq. (3.11). Here, $g$ denotes the gauge coupling
of $SU(2)_{L}$, $\theta_{W}$ is the weak mixing angle, ${\rm
e}=g\sin\theta_{W}$ is the electromagnetic coupling, $T_{3}^{f}=+1/2\,(-1/2)$
for $f=u\,(d,e)$, and $Q_{f}$ is the electric charge of $f$.

The relevant sub-diagrams for the BZ diagrams are shown in Fig. 2.

$\mu$$e$$e$$q$$k$$f$$p_{1}$$p_{2}-k$$p_{2}$$\mu$$e$$e$$q$$k$$f$$p_{1}$$p_{1}+k$$p_{2}$$\mu$$e$$e$$q$$k$$W$$p_{1}$$p_{2}-k$$p_{2}$$\mu$$e$$e$$q$$k$$W$$p_{1}$$p_{1}+k$$p_{2}$
Figure 2: Focused BZ type diagrams with fermion (top diagrams) and $W$ boson
(bottom diagrams) loops.

The top panel diagrams feature a fermion inner loop, while the bottom panel
diagrams have the $W$ boson in the inner loop. In our calculation, we adopt
the Feynman-’t Hooft gauge for the $W$ boson propagator and utilize the
results from ref. [46]. In general, for $W$ boson contributions to the
electron EDM, careful consideration of gauge dependence is necessary, as the
diagrams shown in Fig. 2 alone do not yield a gauge-independent result. A
comprehensive gauge-independent analysis of the BZ diagrams with $W$ bosons
was conducted in ref. [47]. Although numerous diagrams warrant consideration
for gauge independence, we focus solely on the diagrams presented in Fig. 2.
This is because these diagrams already encompass the leading contributions
among the $W$ boson loop diagrams in the Feynman-’t Hooft gauge.

The sum of four diagrams in Fig. 2 leads to a general result of the electron
EDM as

$\displaystyle d_{e}$ $\displaystyle=\sum_{f}d_{e}^{f}+d_{e}^{W}\,,$ ( 3.12)
$\displaystyle d_{e}^{f,W}$
$\displaystyle=\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\int_{x}\int_{y}\int_{k}\frac{F_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})}{(k^{2}-\Delta_{f,W})(k^{2}-m_{G}^{2})(k^{2}+2k\cdot
q-m_{\phi_{j}}^{2})(k^{2}+2k\cdot q-m_{\eta_{a}}^{2})}$ $\displaystyle\hskip
10.00002pt\times\sum_{\Phi=\phi_{j},\eta_{a}}\left(\frac{A_{\alpha\beta}^{L,\Phi(f,W)}k^{\alpha}k^{\beta}+B_{\alpha}^{L,\Phi(f,W)}k^{\alpha}}{k^{2}-2k\cdot
p_{2}}+\frac{A_{\alpha\beta}^{R,\Phi(f,W)}k^{\alpha}k^{\beta}+B_{\alpha}^{R,\Phi(f,W)}k^{\alpha}}{k^{2}+2k\cdot
p_{1}}\right)\,,$ ( 3.13)

where we have used $p_{1,2}^{2}=m_{e}^{2}$ for the outer electron lines,
$q=p_{1}-p_{2}$ represents the four-momentum of the external photon, and
$\Delta_{f,W}$ is defined as

$\displaystyle\Delta_{f,W}\equiv\frac{m_{f,W}^{2}}{x(1-x)}-\frac{2y}{1-x}k\cdot
q\,.$ ( 3.14)

In addition, we have defined a shorthand notation for each integral as

$\displaystyle\int_{x}$
$\displaystyle\equiv\int_{0}^{1}\\!dx\frac{1}{x(1-x)}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \int_{y}\equiv\int_{0}^{1-x}\\!dy\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \int_{k}\equiv\int\\!\frac{d^{4}k}{(2\pi)^{4}}\,,$ (
3.15)

where $x$ and $y$ are Feynman parameters arising from the inner loop. The
terms $A_{\alpha\beta}^{L,\Phi(f,W)}$, $A_{\alpha\beta}^{R,\Phi(f,W)}$,
$B_{\alpha}^{L,\Phi(f,W)}$, and $B_{\alpha}^{R,\Phi(f,W)}$ represent the
relevant contributions to the electron EDM from fermion and $W$ boson loop
diagrams. The explicit expressions for these terms are summarized in appendix
D.1.

The expressions in Eq. (3.13) are general and applicable to any dark sector
models. Before delving into the numerical calculation, we emphasize the
following points:

  * (1)

The momenta $p_{1,2}$ correspond to the external electron momenta. Thus, terms
such as $k\cdot p_{i}$ and $k\cdot q$ in the expression will contribute sub-
dominantly. Upon integrating the loop momentum $k$, these terms lead to
$p_{i}\cdot p_{j}=m_{e}^{2}$ or $p_{i}\cdot q=0$, as $q\cdot q=0$.

  * (2)

Due to the propagators, Eq. (3.13) exhibits a pole around $k^{2}\sim v_{\rm
SM}^{2}$ if the neutral scalar masses in the visible sector are of
$\mathcal{O}(v_{\rm SM})$. Consequently, the dominant contribution to $d_{e}$
is expected to occur when $k^{2}\sim v_{\rm SM}^{2}$.

####  3.4.1 The heavy dark sector limit

With the property (2), the general formula can be simplified when the unknown
dark sector is significantly heavier than the electroweak scale, $M\gg v_{\rm
SM}$. By power expanding $F_{\rm DS}^{\phi_{j}\eta_{a}}(k^{2},M^{2})$ in terms
of $v_{\rm SM}^{2}/M^{2}(\sim k^{2}/M^{2})$, we can simplify $F_{\rm
DS}^{\phi_{j}\eta_{a}}$ as

$\displaystyle F_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})=a_{0}^{\phi_{j}\eta_{a}}(M^{2})+a_{1}^{\phi_{j}\eta_{a}}(M^{2})\frac{(k+q)^{2}}{M^{2}}+\mathcal{O}\left(\frac{(k+q)^{4}}{M^{4}}\right)\,,$
( 3.16)

where $a_{0,1}^{\phi_{j}\eta_{a}}(M^{2})$ are independent of $k$, and each
order is expected to be suppressed by $\mathcal{O}(v_{\rm SM}^{2}/M^{2})$.
With this expansion, the electron EDM can be explicitly calculated for the
heavy dark sector (HDS). The leading results are (see appendix D.2 for
derivation):

$\displaystyle d_{e}^{f,W}$
$\displaystyle\simeq\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\sum_{\Phi=\phi_{j},\eta_{a}}\int_{x}\frac{g_{eG}^{V}\left(C_{E}^{f,W}g_{e\Phi}^{P}-C_{O}^{f,W}g_{e\Phi}^{S}\right)}{16\pi^{2}m_{\phi_{j}}^{2}}$
$\displaystyle\hskip
30.00005pt\times\left[\frac{a_{0}^{\phi_{j}\eta_{a}}(M^{2})}{m_{\phi_{j}}^{2}}j_{3}^{(1)}(r_{f,W},r_{G},r_{\eta_{a}})+\frac{a_{1}^{\phi_{j}\eta_{a}}(M^{2})}{M^{2}}j_{3}^{(2)}(r_{f,W},r_{G},r_{\eta_{a}})\right]\,,$
( 3.17)

where we omit terms proportional to $m_{e}^{2}/m_{\phi_{j}}^{2}\sim
m_{e}^{2}/v_{\rm SM}^{2}$. The integrals of momentum $k$ and the Feynman
parameter $y$ have been calculated explicitly. Here, $r_{f,W,G,\eta_{a}}$ are
mass squared ratios defined as

$\displaystyle
r_{f,W}\equiv\frac{m_{f,W}^{2}}{m_{\phi_{j}}^{2}x(1-x)}\,,\qquad
r_{G}\equiv\frac{m_{G}^{2}}{m_{\phi_{j}}^{2}}\,,\qquad
r_{\eta_{a}}\equiv\frac{m_{\eta_{a}}^{2}}{m_{\phi_{j}}^{2}}\,.$ ( 3.18)

The coefficients $C_{E}^{f,W}$ and $C_{O}^{f,W}$ encapsulate the information
from the inner fermion and $W$ boson loops, and they are defined as [46, 36]:

$\displaystyle C_{E}^{f}$
$\displaystyle=+\frac{N_{C}^{f}}{4\pi^{2}}g_{f\gamma}^{V}g_{fG}^{V}g_{f\Phi}^{S}m_{f}\Bigl{[}x^{2}+(1-x)^{2}\Bigr{]}\,,$
( 3.19) $\displaystyle C_{O}^{f}$
$\displaystyle=-\frac{N_{C}^{f}}{4\pi^{2}}g_{f\gamma}^{V}g_{fG}^{V}g_{f\Phi}^{P}m_{f}\,,$
( 3.20) $\displaystyle C_{E}^{W}$ $\displaystyle=+\frac{\rm
e}{8\pi^{2}}g_{WW\Phi}g_{WWG}\left\\{\left(4-\frac{m_{G}^{2}}{m_{W}^{2}}\right)-\frac{x(1-x)}{2}\left[6-\frac{m_{G}^{2}}{m_{W}^{2}}+\left(1-\frac{m_{G}^{2}}{2m_{W}^{2}}\right)\frac{m_{\Phi}^{2}}{m_{W}^{2}}\right]\right\\}\,,$
( 3.21) $\displaystyle C_{O}^{W}$ $\displaystyle=0\,,$ ( 3.22)

where $g_{WW\Phi}=\frac{2m_{W}^{2}}{v_{\rm SM}}\xi_{\phi_{j}}^{W}\,(0)$ for
CP-even (CP-odd) scalar, and $g_{WWG}={\rm
e}\cos\theta_{W}/\sin\theta_{W}\,({\rm e})$ for $G=Z(\gamma)$. The function
$j_{3}^{(1,2)}(r_{f,W},r_{G},r_{\eta_{a}})$, arising from the integration of
Feynman parameters, is defined by

$\displaystyle j_{3}^{(n)}(r_{f,W},r_{G},r_{\eta_{a}})$
$\displaystyle\equiv\frac{r_{f,W}^{n}\mathop{\rm
ln}r_{f,W}}{(r_{f,W}-1)(r_{f,W}-r_{G})(r_{f,W}-r_{\eta_{a}})}+\frac{r_{G}^{n}\mathop{\rm
ln}r_{G}}{(r_{G}-1)(r_{G}-r_{f,W})(r_{G}-r_{\eta_{a}})}$ $\displaystyle\hskip
60.00009pt+\frac{r_{\eta_{a}}^{n}\mathop{\rm
ln}r_{\eta_{a}}}{(r_{\eta_{a}}-1)(r_{\eta_{a}}-r_{f,W})(r_{\eta_{a}}-r_{G})}\,.$
( 3.23)

For the case where the gauge boson is the photon, i.e., $G=\gamma$
($m_{G}=0$), the function $j_{3}$ becomes:

$\displaystyle
j_{3}^{(n)}(r_{f,W},0,r_{\eta_{a}})=\frac{r_{f,W}^{n-1}\mathop{\rm
ln}r_{f,W}}{(r_{f,W}-1)(r_{f,W}-r_{\eta_{a}})}+\frac{r_{\eta_{a}}^{n-1}\mathop{\rm
ln}r_{\eta_{a}}}{(r_{\eta_{a}}-1)(r_{\eta_{a}}-r_{f,W})}\,.$ ( 3.24)

##  4 Numerical analysis for benchmark models

In the previous section, we have derived the formula for the electron EDM in a
general CP-violating dark sector. Let us now consider two specific dark sector
models. The first model consists of a complex singlet scalar whose potential
contains CP violation (the scalar dark sector model). The second model
involves a vector-like singlet fermion with a CP-violating Yukawa coupling and
a complex singlet scalar portal to the Higgs sector (the fermion dark sector
model). For both models, we numerically estimate the size of the electron EDM
and discuss the validity of the expression in the HDS limit.

###  4.1 Scalar dark sector model

We introduce a complex singlet scalar denoted as $S$. Its mass term and self-
interactions are described by the potential $V_{\rm DS}$. In addition, we
allow for interactions between the scalar $S$ and the Higgs doublets of the
2HDM through the potential $V_{H\mathcal{O}}$. The expressions for $V_{\rm
DS}$ and $V_{H\mathcal{O}}$ are given by

$\displaystyle V_{\rm DS}$
$\displaystyle=m_{S}^{2}|S|^{2}+\frac{\lambda_{S}}{2}|S|^{4}+\Bigl{[}B_{S}S^{2}+D_{S}S^{4}+{\rm
h.c.}\Bigr{]}\,,$ ( 4.1) $\displaystyle V_{H\mathcal{O}}$
$\displaystyle=\lambda_{1S}H_{1}^{\dagger}H_{1}|S|^{2}+\lambda_{2S}H_{2}^{\dagger}H_{2}|S|^{2}+\Bigl{[}\lambda_{12S}H_{1}^{\dagger}H_{2}S^{2}+\lambda_{21S}H_{2}^{\dagger}H_{1}S^{2}+{\rm
h.c.}\Bigr{]}\,.$ ( 4.2)

In our analysis, we parameterize the complex scalar as
$S=(s_{1}+is_{2})/\sqrt{2}$. To simplify and focus on relevant terms in
$V_{\rm DS}$ and $V_{H\mathcal{O}}$, we impose a discrete $Z_{4}$ symmetry,
and its charges are summarized in Table 4. According to the charge assignment
of $S$, $B_{S}$ is a soft $Z_{4}$ breaking term.333 Due to the remaining
$Z_{2}$ symmetry, $S$ may become a DM candidate.

$H_{1}$ | $H_{2}$ | $S$ | $Q_{L}$ | $u_{R}$ | $d_{R}$ | $L_{L}$ | $e_{R}$
---|---|---|---|---|---|---|---
$i^{2}$ | $1$ | $i$ | $1$ | $1$ | $z_{4}^{d}$ | $1$ | $z_{4}^{e}$
Table 4: We present an example of $Z_{4}$ charges for the scalar dark sector
model. The exact values of $z_{4}^{d}$ and $z_{4}^{e}$ for the right-handed
$d$ quark and electron depend on the specific type of the 2HDM being
considered.

There are 6 complex parameters in the scalar potential of the model:
$m_{12}^{2}$ and $\lambda_{5}$ in Eq. (2.1) and $B_{S}$, $D_{S}$,
$\lambda_{12S}$ and $\lambda_{21S}$ in $V_{\rm DS}+V_{H\mathcal{O}}$. To
explore CP violation, it is advantageous to identify rephasing-invariant
phases, as studied in refs. [48, 49]. The details are discussed in appendix
C.1, where we have five invariant phases. For our analysis, the relevant
invariant phase in the dark sector is

$\displaystyle\theta_{\rm phys}=\arg\left[B_{S}^{2}D_{S}^{*}\right]\,.$ ( 4.3)

While $D_{S}$ is, in principle, a complex parameter, we can select the phase
of $S$ to render it real. Consequently, in this scenario, the physical phase
$\theta_{\rm phys}$ is dictated by the phase of $B_{S}$, resulting in tree-
level CP violation in the mass mixing matrix for the $S$ scalars. In appendix
C.1, another invariant phase
$\arg[B_{S}^{2}\lambda_{12S}^{*}\lambda_{21S}^{*}]$ can also induce the same
tree-level CP violation once $H_{1,2}$ acquire their VEVs. However, given our
assumption of CP violation happening in the dark sector only, we set
$\arg[B_{S}^{2}\lambda_{12S}^{*}\lambda_{21S}^{*}]=0$ to maintain simplicity
in the CP violation pattern. Thus, in our analysis, we turn off all the
physical phases other than Eq. (4.3).

Our assumption (i) in section 2.2 requires that the scalar field $S$ does not
acquire a nonzero VEV. Therefore, the CP-violating phase in $B_{S}$
contributes to the mass matrix of the dark scalars $s_{1,2}$, resulting in CP-
violating mixing. Explicitly, the mass matrix for the complex scalar $S$ is
given in the basis of $(s_{1},s_{2})$ as

$\displaystyle M_{s}^{2}=\overline{m}_{S}^{2}\mathbf{1}_{2\times
2}+\begin{pmatrix}{\rm Re}(\overline{B}_{S})&-{\rm Im}(\overline{B}_{S})\\\
-{\rm Im}(\overline{B}_{S})&-{\rm Re}(\overline{B}_{S})\end{pmatrix}\,,$ (
4.4)

with $\overline{m}_{S}^{2}\equiv
m_{S}^{2}+\frac{1}{2}\lambda_{1S}v_{1}^{2}+\frac{1}{2}\lambda_{2S}v_{2}^{2}$
and $\overline{B}_{S}\equiv 2B_{S}+(\lambda_{12S}+\lambda_{21S})v_{1}v_{2}$.
Note that since the neutral scalars in the 2HDM sector and $s_{1,2}$ do not
mix with each other, the SM Yukawa couplings and the $W$ boson couplings are
not modified from Eq. (2.7) with coefficients $\xi_{h,H,A}^{f}$ in Table 2 and
Eq. (2.8). As a result, we obtain the following new mass eigenstates
$\varphi_{1,2}$:

$\displaystyle\varphi_{1}=s_{1}c_{\theta_{s}}+s_{2}s_{\theta_{s}}\,,\qquad\varphi_{2}=-s_{1}s_{\theta_{s}}+s_{2}c_{\theta_{s}}\,,$
( 4.5)

where $\theta_{s}$ is a mixing angle with
$m_{\varphi_{1}}^{2}<m_{\varphi_{2}}^{2}$, and its value is determined by

$\displaystyle\sin 2\theta_{s}$ $\displaystyle=\frac{{\rm
Im}(\overline{B}_{S})}{|\overline{B}_{S}|}\,.$ ( 4.6)

When we assume the portal couplings $\lambda_{12S}$ and $\lambda_{21S}$ to be
real, $\theta_{s}$ becomes a phase of $B_{S}$ and corresponds to the physical
CP-violating phase. The mass eigenvalues are given by

$\displaystyle
m_{\varphi_{1}}^{2}=\overline{m}_{S}^{2}-|\overline{B}_{S}|\,,\qquad
m_{\varphi_{2}}^{2}=\overline{m}_{S}^{2}+|\overline{B}_{S}|\,.$ ( 4.7)

Incorporating all the mass eigenstates into the potentials, we can derive the
corresponding neutral scalar cubic and quartic couplings as

$\displaystyle V_{\rm scl}=V_{H}+V_{H\mathcal{O}}+V_{\rm DS}\supset\,$
$\displaystyle\frac{v_{\rm
SM}}{2}\Bigl{(}\lambda_{hjk}h^{0}\varphi_{j}\varphi_{k}+\lambda_{Hjk}H^{0}\varphi_{j}\varphi_{k}+\lambda_{Ajk}A^{0}\varphi_{j}\varphi_{k}\Bigr{)}$
$\displaystyle+\frac{1}{2}\Bigl{(}\lambda_{hAjj}h^{0}A^{0}\varphi_{j}\varphi_{j}+\lambda_{HAjj}H^{0}A^{0}\varphi_{j}\varphi_{j}\Bigr{)}\,,$
( 4.8)

where $j,k=1,2$ and each coupling is summarized in Table 5.

$\lambda_{h11}$ | $-\lambda_{1S}s_{\alpha}c_{\beta}+\lambda_{2S}c_{\alpha}s_{\beta}+\left(\lambda_{12S}+\lambda_{21S}\right)c_{\alpha+\beta}c_{2\theta_{s}}$
---|---
$\lambda_{h12}$ | $-2\left(\lambda_{12S}+\lambda_{21S}\right)c_{\alpha+\beta}s_{2\theta_{s}}$
$\lambda_{h22}$ | $-\lambda_{1S}s_{\alpha}c_{\beta}+\lambda_{2S}c_{\alpha}s_{\beta}-\left(\lambda_{12S}+\lambda_{21S}\right)c_{\alpha+\beta}c_{2\theta_{s}}$
$\lambda_{H11}$ | $\lambda_{1S}c_{\alpha}c_{\beta}+\lambda_{2S}s_{\alpha}s_{\beta}+\left(\lambda_{12S}+\lambda_{21S}\right)s_{\alpha+\beta}c_{2\theta_{s}}$
$\lambda_{H12}$ | $-2\left(\lambda_{12S}+\lambda_{21S}\right)s_{\alpha+\beta}s_{2\theta_{s}}$
$\lambda_{H22}$ | $\lambda_{1S}c_{\alpha}c_{\beta}+\lambda_{2S}s_{\alpha}s_{\beta}-\left(\lambda_{12S}+\lambda_{21S}\right)s_{\alpha+\beta}c_{2\theta_{s}}$
$\lambda_{A11}$ | $-\left(\lambda_{12S}-\lambda_{21S}\right)s_{2\theta_{s}}$
$\lambda_{A12}$ | $-2\left(\lambda_{12S}-\lambda_{21S}\right)c_{2\theta_{s}}$
$\lambda_{A22}$ | $\left(\lambda_{12S}-\lambda_{21S}\right)s_{2\theta_{s}}$
$\lambda_{hA11}$ | $\left(\lambda_{12S}-\lambda_{21S}\right)s_{\alpha-\beta}s_{2\theta_{s}}$
$\lambda_{hA22}$ | $-\left(\lambda_{12S}-\lambda_{21S}\right)s_{\alpha-\beta}s_{2\theta_{s}}$
$\lambda_{HA11}$ | $-\left(\lambda_{12S}-\lambda_{21S}\right)c_{\alpha-\beta}s_{2\theta_{s}}$
$\lambda_{HA22}$ | $\left(\lambda_{12S}-\lambda_{21S}\right)c_{\alpha-\beta}s_{2\theta_{s}}$
Table 5: The relevant neutral scalar cubic and quartic couplings. For
instance, $\lambda_{h12}$ represents the trilinear coupling involving the mass
eigenstate of the SM Higgs and the dark sector scalars $\varphi_{1}$ and
$\varphi_{2}$.

With the cubic and quartic couplings determined, we can now explicitly
calculate the one-loop contributions from the dark sector scalars to the two
effective CP-violating mixing terms $F_{\rm DS}^{hA}$ and $F_{\rm DS}^{HA}$,
as shown in the diagrams of Fig. 3,

$\displaystyle F_{\rm DS}^{hA}(k^{2})$ $\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})+{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})\biggr{.}$ $\displaystyle\hskip
135.0002pt\biggl{.}-2{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{2}})+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{k^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\biggr{)}$
$\displaystyle\hskip 60.00009pt+\lambda_{hS}\left({\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})-{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})-\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\right)\Biggl{.}\Biggr{]}$
$\displaystyle\hskip
11.99998pt+\frac{\lambda_{12S}^{-}}{32\pi^{2}}s_{\alpha-\beta}s_{2\theta_{s}}\left[\left(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}\right)\Lambda_{\rm
UV}^{2}+m_{\varphi_{1}}^{2}\mathop{\rm
ln}\frac{\mu^{2}}{m_{\varphi_{1}}^{2}}-m_{\varphi_{2}}^{2}\mathop{\rm
ln}\frac{\mu^{2}}{m_{\varphi_{2}}^{2}}\right]\,,$ ( 4.9) $\displaystyle F_{\rm
DS}^{HA}(k^{2})$ $\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})+{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})\biggr{.}$ $\displaystyle\hskip
135.0002pt\biggl{.}-2{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{2}})+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{k^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\biggr{)}$
$\displaystyle\hskip 60.00009pt+\lambda_{HS}\left({\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})-{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})-\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\right)\Biggl{.}\Biggr{]}$
$\displaystyle\hskip
11.99998pt-\frac{\lambda_{12S}^{-}}{32\pi^{2}}c_{\alpha-\beta}s_{2\theta_{s}}\left[\left(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}\right)\Lambda_{\rm
UV}^{2}+m_{\varphi_{1}}^{2}\mathop{\rm
ln}\frac{\mu^{2}}{m_{\varphi_{1}}^{2}}-m_{\varphi_{2}}^{2}\mathop{\rm
ln}\frac{\mu^{2}}{m_{\varphi_{2}}^{2}}\right]\,.$ ( 4.10)

Here, we have defined the shorthand notations
$\lambda_{12S}^{\pm}\equiv\lambda_{12S}\pm\lambda_{21S}$,
$\lambda_{hS}\equiv-\lambda_{1S}s_{\alpha}c_{\beta}+\lambda_{2S}c_{\alpha}s_{\beta}$
and
$\lambda_{HS}\equiv\lambda_{1S}c_{\alpha}c_{\beta}+\lambda_{2S}s_{\alpha}s_{\beta}$.
The function ${\rm DiscB}(k^{2},m_{1},m_{2})$ is defined as

$\displaystyle{\rm DiscB}(k^{2},m_{1},m_{2})$
$\displaystyle=\frac{\lambda^{1/2}(k^{2},m_{1}^{2},m_{2}^{2})}{k^{2}}\mathop{\rm
ln}\frac{m_{1}^{2}+m_{2}^{2}-k^{2}+\lambda^{1/2}(k^{2},m_{1}^{2},m_{2}^{2})}{2m_{1}m_{2}}\,,$
( 4.11) $\displaystyle\lambda(k^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle=k^{4}+m_{1}^{4}+m_{2}^{4}-2k^{2}m_{1}^{2}-2k^{2}m_{2}^{2}-2m_{1}^{2}m_{2}^{2}\,,$
( 4.12)

and $\Lambda_{\rm UV}$ and $\mu$ are a UV divergent parameter and ’t Hooft
mass parameter, respectively.

${h,H}$$A$$\varphi_{j}$$\varphi_{k}$${h,H}$$A$$\varphi_{1,2}$ Figure 3: One-
loop diagrams relevant to the effective $h$-$A$ and $H$-$A$ mixings induced by
$\varphi_{1,2}$ in the scalar dark sector model.

To obtain a finite physical prediction, we need to deal with the UV divergent
part appropriately. For this purpose, let us consider the wave function
renormalization for $h$-$A$ and $H$-$A$ self-energies [50, 51]. First, we
define the renormalized two-point functions $\hat{\Sigma}_{hA,HA}(k^{2})$ as

$\displaystyle\hat{\Sigma}_{hA}(k^{2})$ $\displaystyle\equiv F_{\rm
DS}^{hA}(k^{2})+\frac{1}{2}(k^{2}-m_{h}^{2})\delta
Z_{hA}+\frac{1}{2}(k^{2}-m_{A}^{2})\delta Z_{Ah}\,,$ ( 4.13)
$\displaystyle\hat{\Sigma}_{HA}(k^{2})$ $\displaystyle\equiv F_{\rm
DS}^{HA}(k^{2})+\frac{1}{2}(k^{2}-m_{H}^{2})\delta
Z_{HA}+\frac{1}{2}(k^{2}-m_{A}^{2})\delta Z_{AH}\,.$ ( 4.14)

In the on-shell subtraction scheme, we can determine $\delta Z_{hA,HA}$ and
$\delta Z_{Ah,AH}$ by

$\displaystyle{\rm
Re}\bigl{\\{}\hat{\Sigma}_{hA}(m_{h}^{2})\bigr{\\}}=0\,,\quad{\rm
Re}\bigl{\\{}\hat{\Sigma}_{HA}(m_{H}^{2})\bigr{\\}}=0\,,\quad{\rm
Re}\bigl{\\{}\hat{\Sigma}_{hA,HA}(m_{A}^{2})\bigr{\\}}=0\,.$ ( 4.15)

By solving these equations, the renormalized two-point functions
$\hat{\Sigma}_{hA,HA}(k^{2})$ are given by

$\displaystyle\hat{\Sigma}_{hA}(k^{2})$ $\displaystyle=F_{\rm
DS}^{hA}(k^{2})+\frac{1}{m_{h}^{2}-m_{A}^{2}}\Bigl{[}\left(m_{A}^{2}-k^{2}\right){\rm
Re}\\!\left\\{F_{\rm
DS}^{hA}(m_{h}^{2})\right\\}+\left(k^{2}-m_{h}^{2}\right){\rm
Re}\\!\left\\{F_{\rm DS}^{hA}(m_{A}^{2})\right\\}\Bigr{]}\,,$ ( 4.16)
$\displaystyle\hat{\Sigma}_{HA}(k^{2})$ $\displaystyle=F_{\rm
DS}^{HA}(k^{2})+\frac{1}{m_{H}^{2}-m_{A}^{2}}\Bigl{[}\left(m_{A}^{2}-k^{2}\right){\rm
Re}\\!\left\\{F_{\rm
DS}^{HA}(m_{H}^{2})\right\\}+\left(k^{2}-m_{H}^{2}\right){\rm
Re}\\!\left\\{F_{\rm DS}^{HA}(m_{A}^{2})\right\\}\Bigr{]}\,.$ ( 4.17)

One can easily check that the divergent parts (and also other constant terms
concerning the momentum $k$) are canceled, and hence, $\hat{\Sigma}_{hA,HA}$
give finite results. For numerical predictions, we need to replace $F_{\rm
DS}^{hA}\to\hat{\Sigma}_{hA}$ and $F_{\rm DS}^{HA}\to\hat{\Sigma}_{HA}$ in the
expressions that we have obtained in the previous section.

In the HDS approximation, the dark sector scalar masses
$m_{\varphi_{1}},m_{\varphi_{2}}$ are much larger than the masses in the 2HDM.
In this case, we can expand $\hat{\Sigma}_{hA(HA)}(k^{2})$ as

$\displaystyle\hat{\Sigma}_{hA(HA)}(k^{2})=a_{0}^{hA(HA)}+a_{1}^{hA(HA)}\frac{k^{2}}{m_{\varphi_{1}}m_{\varphi_{2}}}+\mathcal{O}\left(\frac{k^{4}}{m_{\varphi_{i}}^{4}}\right)\,,$
( 4.18)

where we choose $M^{2}$ in Eq. (3.16) to be $m_{\varphi_{1}}m_{\varphi_{2}}$.
Each coefficient is given explicitly as

$\displaystyle a_{0}^{hA}$ $\displaystyle=\frac{v_{\rm
SM}^{2}\lambda_{12S}^{-}}{64\pi^{2}}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}2+\left\\{\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{m_{h}^{2}}+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{m_{A}^{2}}-\frac{m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2}}{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}\right\\}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}$ $\displaystyle\hskip
90.00014pt+\frac{1}{m_{h}^{2}-m_{A}^{2}}\Bigl{[}m_{h}^{2}\mathcal{F}_{1}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})-m_{A}^{2}\mathcal{F}_{1}(m_{h},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{]}\biggl{.}\biggr{)}$
$\displaystyle\hskip
65.00009pt+\frac{\lambda_{hS}}{m_{h}^{2}-m_{A}^{2}}\Bigl{(}m_{h}^{2}\mathcal{F}_{2}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})-m_{A}^{2}\mathcal{F}_{2}(m_{h},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{)}\Biggl{.}\Biggr{]}\,,$
( 4.19) $\displaystyle a_{1}^{hA}$ $\displaystyle=-\frac{v_{\rm
SM}^{2}\lambda_{12S}^{-}}{64\pi^{2}}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12s}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}\frac{(m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2})(m_{\varphi_{1}}^{4}-8m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}+m_{\varphi_{2}}^{4})}{6m_{\varphi_{1}}m_{\varphi_{2}}(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{2}}$
$\displaystyle\hskip
120.00018pt+\left\\{\frac{(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})m_{\varphi_{1}}m_{\varphi_{2}}}{m_{h}^{2}m_{A}^{2}}+\frac{2m_{\varphi_{1}}^{3}m_{\varphi_{2}}^{3}}{(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{3}}\right\\}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}$ $\displaystyle\hskip
120.00018pt-\frac{m_{\varphi_{1}}m_{\varphi_{2}}}{m_{h}^{2}-m_{A}^{2}}\Bigl{[}\mathcal{F}_{1}(m_{h},m_{\varphi_{1}},m_{\varphi_{2}})-\mathcal{F}_{1}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{]}\biggl{.}\biggr{)}$
$\displaystyle\hskip
35.00005pt-\lambda_{hS}\biggl{(}\biggr{.}\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{6m_{\varphi_{1}}m_{\varphi_{2}}}+\frac{m_{\varphi_{1}}m_{\varphi_{2}}}{m_{h}^{2}-m_{A}^{2}}\Bigl{[}\mathcal{F}_{2}(m_{h},m_{\varphi_{1}},m_{\varphi_{2}})-\mathcal{F}_{2}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{]}\biggl{.}\biggr{)}\Biggl{.}\Biggr{]}\,,$
( 4.20) $\displaystyle a_{0}^{HA}$ $\displaystyle=\frac{v_{\rm
SM}^{2}\lambda_{12S}^{-}}{64\pi^{2}}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}2+\left\\{\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{m_{H}^{2}}+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{m_{A}^{2}}-\frac{m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2}}{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}\right\\}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}$ $\displaystyle\hskip
90.00014pt+\frac{1}{m_{H}^{2}-m_{A}^{2}}\Bigl{[}m_{H}^{2}\mathcal{F}_{1}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})-m_{A}^{2}\mathcal{F}_{1}(m_{H},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{]}\biggl{.}\biggr{)}$
$\displaystyle\hskip
65.00009pt+\frac{\lambda_{HS}}{m_{H}^{2}-m_{A}^{2}}\Bigl{(}m_{H}^{2}\mathcal{F}_{2}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})-m_{A}^{2}\mathcal{F}_{2}(m_{H},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{)}\Biggl{.}\Biggr{]}\,,$
( 4.21) $\displaystyle a_{1}^{HA}$ $\displaystyle=-\frac{v_{\rm
SM}^{2}\lambda_{12S}^{-}}{64\pi^{2}}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12s}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}\frac{(m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2})(m_{\varphi_{1}}^{4}-8m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}+m_{\varphi_{2}}^{4})}{6m_{\varphi_{1}}m_{\varphi_{2}}(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{2}}$
$\displaystyle\hskip
120.00018pt+\left\\{\frac{(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})m_{\varphi_{1}}m_{\varphi_{2}}}{m_{H}^{2}m_{A}^{2}}+\frac{2m_{\varphi_{1}}^{3}m_{\varphi_{2}}^{3}}{(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{3}}\right\\}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}$ $\displaystyle\hskip
120.00018pt-\frac{m_{\varphi_{1}}m_{\varphi_{2}}}{m_{H}^{2}-m_{A}^{2}}\Bigl{[}\mathcal{F}_{1}(m_{H},m_{\varphi_{1}},m_{\varphi_{2}})-\mathcal{F}_{1}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{]}\biggl{.}\biggr{)}$
$\displaystyle\hskip
35.00005pt-\lambda_{HS}\biggl{(}\biggr{.}\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{6m_{\varphi_{1}}m_{\varphi_{2}}}+\frac{m_{\varphi_{1}}m_{\varphi_{2}}}{m_{H}^{2}-m_{A}^{2}}\Bigl{[}\mathcal{F}_{2}(m_{H},m_{\varphi_{1}},m_{\varphi_{2}})-\mathcal{F}_{2}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{]}\biggl{.}\biggr{)}\Biggl{.}\Biggr{]}\,,$
( 4.22)

where we have defined the following functions:

$\displaystyle\mathcal{F}_{1}(m,m_{\varphi_{1}},m_{\varphi_{2}})$
$\displaystyle\equiv{\rm Re}\\!\left\\{{\rm
DiscB}(m^{2},m_{\varphi_{1}},m_{\varphi_{1}})\right\\}+{\rm Re}\\!\left\\{{\rm
DiscB}(m^{2},m_{\varphi_{2}},m_{\varphi_{2}})\right\\}$ $\displaystyle\hskip
25.00003pt-2{\rm Re}\\!\left\\{{\rm
DiscB}(m^{2},m_{\varphi_{1}},m_{\varphi_{2}})\right\\}\,,$ ( 4.23)
$\displaystyle\mathcal{F}_{2}(m,m_{\varphi_{1}},m_{\varphi_{2}})$
$\displaystyle\equiv{\rm Re}\\!\left\\{{\rm
DiscB}(m^{2},m_{\varphi_{1}},m_{\varphi_{1}})\right\\}-{\rm Re}\\!\left\\{{\rm
DiscB}(m^{2},m_{\varphi_{2}},m_{\varphi_{2}})\right\\}\,.$ ( 4.24)

Since $a_{0}^{hA(HA)}=\hat{\Sigma}_{\rm DS}^{hA(HA)}(k^{2}\to 0)$, they are
used to calculate the coefficient of the four-fermion operator
$C_{ff^{\prime}}$ by inserting $a_{0}^{hA,HA}$ into Eq. (3.8). Therefore, we
have obtained the electron EDM in the HDS approximation as well as the
coefficient $C_{S}$. By applying the limit of the JILA experiment in Eq.
(3.3), one can establish bounds on the parameters for the scalar dark sector
model.

####  4.1.1 Numerical results

Let us now perform the numerical analysis for the scalar dark sector model. We
choose the following values for the relevant masses and parameters:

$\displaystyle m_{h}=125.25\leavevmode\nobreak\ {\rm GeV},\leavevmode\nobreak\
\leavevmode\nobreak\ m_{H}=1500\leavevmode\nobreak\ {\rm
GeV},\leavevmode\nobreak\ \leavevmode\nobreak\ m_{A}=1550\leavevmode\nobreak\
{\rm GeV},\leavevmode\nobreak\ \leavevmode\nobreak\
\alpha=-0.02,\leavevmode\nobreak\ \leavevmode\nobreak\ \tan\beta=50,$
$\displaystyle\theta_{s}=\frac{\pi}{4},\leavevmode\nobreak\
\leavevmode\nobreak\ \lambda_{1S}=1,\leavevmode\nobreak\ \leavevmode\nobreak\
\lambda_{2S}=2,\leavevmode\nobreak\ \leavevmode\nobreak\
\lambda_{12S}=2,\leavevmode\nobreak\ \leavevmode\nobreak\ \lambda_{21S}=0.1.$
( 4.25)

Note that the value of $\theta_{s}$ corresponds to the maximal choice for the
physical CP phase. With this parameter choice, we have found that the $C_{S}$
contribution is several orders of magnitude smaller than the electron EDM
$d_{e}$ from the BZ diagrams. Hence, we simply neglect the $C_{S}$ term in the
calculation.

Figure 4: Numerical results of the electron EDM in the scalar dark sector
model. The left panel shows the results of Type-I (dashed lines) and Type-II
(solid lines) 2HDMs, while the right panel shows those of Type-X (solid lines)
and Type-Y (dashed lines) 2HDMs. The color means different choices of $\Delta
m_{\varphi}\equiv m_{\varphi_{2}}-m_{\varphi_{1}}$: $\Delta m_{\varphi}=10$
GeV (red), $\Delta m_{\varphi}=100$ GeV (green) and $\Delta m_{\varphi}=1000$
GeV (blue). The gray shaded region is excluded by the current electron EDM
bound from the JILA experiment [18], and the brown shaded region is
constrained by the SM Higgs invisible decay, $h^{0}\to{\rm inv}$ [52].

Fig. 4 shows the size of the electron EDM $d_{e}$ in terms of
$m_{\varphi_{1}}$, with 3 choices of mass difference for $m_{\varphi_{1,2}}$:
$\Delta m_{\varphi}\equiv m_{\varphi_{2}}-m_{\varphi_{1}}=10$ GeV (red lines),
$100$ GeV (green lines) and $1000$ GeV (blue lines). The left panel
corresponds to the cases of Type-I (dashed lines) and Type-II (solid lines),
while the right panel shows those in Type-X (solid lines) and Type-Y (dashed
lines). The shaded regions are constrained by the JILA experiment for the
electron EDM bound [18] (gray) and the SM Higgs invisible decay
$h\to\varphi_{i}\varphi_{j}$, which leads to the exclusion of
$m_{\varphi_{1}}<m_{h}/2$ (brown).

The value of the electron EDM $d_{e}$ is obtained from the general expression
in Eq. (3.13) by omitting $k\cdot q$ and $k\cdot p_{1,2}$ terms which give
small corrections of $\mathcal{O}\left(m_{e}^{2}/v_{\rm SM}^{2}\right)$. We
take account of the BZ diagrams with the 3rd generation fermions and the $W$
boson in the inner loop.

As expected from Table 2, Type-II and Type-X 2HDMs have enough $\tan\beta$
enhancement from the electron-neutral scalar couplings to reach the current
upper bound on $|d_{e}|$, while Type-I and Type-Y 2HDMs predict several orders
of magnitude smaller $|d_{e}|$, due to the $\tan\beta$ suppression from those
couplings. We can see from Fig. 4 that Type-Y 2HDM predicts slightly larger
$d_{e}$ than Type-I 2HDM, due to the enhancement from the bottom quark
couplings to the neutral scalars. In addition, Type-II and Type-X 2HDM cases
exhibit similar predictions for $|d_{e}|$ because the primary contribution to
$|d_{e}|$ arises from the top quark and $W$ boson loop processes which are the
same for both cases. In the case of a large $\tan\beta$, the value of
$\lambda_{1S}$ becomes irrelevant for $|d_{e}|$, while the values of
$\lambda_{2S}$ and $\lambda_{12S}^{-}$ are crucial to enhance $|d_{e}|$. The
value of $|d_{e}|$ is also enhanced by the hierarchy between $m_{\varphi_{1}}$
and $m_{\varphi_{2}}$, as observed in Fig. 4. This enhancement occurs because
$\hat{\Sigma}_{hA(HA)}(k^{2})$ are amplified when $m_{\varphi_{2}}\gg
m_{\varphi_{1}}$, owing to the logarithmic terms involving the ratio of their
masses, which can be achieved by choosing
$\overline{m}_{S}^{2}\simeq|\overline{B}_{S}|$ (see Eq. (4.7)). As a result,
we obtain a sizable electron EDM for Type-II and Type-X 2HDMs: for instance,
by choosing $(m_{\varphi_{1}},m_{\varphi_{2}})=(300\leavevmode\nobreak\ {\rm
GeV},1300\leavevmode\nobreak\ {\rm GeV})$, $|d_{e}|\sim 4\times
10^{-30}\,e\,{\rm cm}$ is predicted. Even for the case of $\Delta
m_{\varphi}=10$ GeV, a large $|d_{e}|$ very close to the current upper limit
can be obtained for $m_{\varphi_{1}}\simeq 100$ GeV.

####  4.1.2 Discussion on the HDS approximation

If $m_{\varphi_{1,2}}$ are much heavier than all neutral scalars of the 2HDM
sector, the electron EDM is given by the HDS approximation shown in Eq. (3.17)
along with Eqs. (4.19)-(4.22). Fig. 5 shows the numerical estimates of the
electron EDM in Type-II 2HDM for $\Delta m_{\varphi}=100$ GeV, with (red
dotted) and without (solid) the HDS approximation. Similar behaviors can be
seen for other types and/or other $\Delta m_{\varphi}$ cases. We can see from
the figure that for large $m_{\varphi_{1}}$, both curves decrease
monotonically. However, the values of $|d_{e}|$ are not close to each other.
Moreover, even for the mass of $m_{\varphi_{1}}=\mathcal{O}(v_{\rm SM})$, the
HDS approximation looks good. To explain these reasons, we first parameterize
$d_{e}$ as

$\displaystyle d_{e}=d_{e}^{\rm Org.}+d_{e}^{\rm Ren.}\,,$ ( 4.26)

where $d_{e}^{\rm Org.}$ and $d_{e}^{\rm Ren.}$ denote the finite part of
$F_{\rm DS}^{hA,HA}$ and the finite part of terms induced by renormalization
for $F_{\rm DS}^{hA,HA}$ in Eqs. (4.16) and (4.17), whose expressions are
summarized in appendix D.3. For the HDS approximation result, we add the
superscript “HDS”,

$\displaystyle d_{e}^{\rm HDS}=d_{e}^{\rm Org.,HDS}+d_{e}^{\rm Ren.,HDS}\,.$ (
4.27)

Note that $d_{e}^{\rm Ren.}$ and $d_{e}^{\rm Ren.,HDS}$ are the same at the
analytical level. Then, the difference between the results with and without
the HDS approximation comes from the accuracy of calculations in $d_{e}^{\rm
Org.}$ and $d_{e}^{\rm Org.,HDS}$.

Figure 5: Numerical results of the electron EDM with (dotted line) and without
(solid line) the HDS approximation in the scalar dark sector model.

To quantify the accuracy of the HDS approximation, we now define

$\displaystyle\rho_{d_{e}}^{\rm Org.}\equiv\frac{d_{e}^{\rm
Org.,HDS}}{d_{e}^{\rm Org.}}\,,$ ( 4.28)

and plot $|1-\rho_{d_{e}}^{\rm Org.}|$ in the left panel of Fig. 6. We can see
that the HDS approximation agrees within 1% when $m_{\varphi_{1}}\gtrsim
4.5\times 10^{4}$ GeV. However, as parameterized in Eq. (4.26), the physical
prediction is the sum of $d_{e}^{\rm Org.}$ and $d_{e}^{\rm Ren.}$, and
therefore, a cancellation between them may affect the accuracy of the physical
results. The right panel of Fig. 6 shows this feature, by defining

$\displaystyle\epsilon_{d_{e}}\equiv\frac{d_{e}}{d_{e}^{\rm
Org.}}\,,\qquad\epsilon_{d_{e}^{\rm HDS}}\equiv\frac{d_{e}^{\rm
HDS}}{d_{e}^{\rm Org.,HDS}}\,.$ ( 4.29)

From this plot, we can understand that $|\epsilon_{d_{e}}|$ for the HDS
approximation is much smaller than that for the result without the HDS
approximation for the heavy $m_{\varphi_{1}}$ region. This means that the
cancellation in the HDS approximation is strong, and hence, $d_{e}^{\rm HDS}$
is predicted to be smaller than $d_{e}$, as one can see from Fig. 5. To
improve the accuracy of $d_{e}$, we need to take more and more higher accuracy
in the numerical integration in $d_{e}^{\rm Org.}$, but it takes lots of
calculation costs. We would like to emphasize that for such a heavy mass
range, the prediction is far below the current upper limit, and hence, we do
not need to improve the accuracy of the calculation unless the electron EDM
bound is significantly improved.

Figure 6: The left panel indicates the accuracy of the HDS approximation with
$\rho_{d_{e}}^{\rm Org.}\equiv d_{e}^{\rm Org.,HDS}/d_{e}^{\rm Org.}$, while
the right panel shows the degree of cancellation in $d_{e}$ (solid line) and
$d_{e}^{\rm HDS}$ (dotted line) with $\epsilon_{d_{e}}\equiv d_{e}/d_{e}^{\rm
Org.}$ and $\epsilon_{d_{e}^{\rm HDS}}\equiv d_{e}^{\rm HDS}/d_{e}^{\rm
Org.,HDS}$. Here, we focus on the Type-II 2HDM with $\Delta m_{\varphi}=100$
GeV.

Accidentally, the HDS approximation is good even for
$m_{\varphi_{1}}=\mathcal{O}(v_{\rm SM})$. This is because around this mass
range, the dominant contribution to $d_{e}$ and $d_{e}^{\rm HDS}$ is coming
from $d_{e}^{\rm Ren.}$ which is unchanged with and without the HDS
approximation. In the right panel of Fig. 6, we can see this feature directly:
$|\epsilon_{d_{e}^{\rm(HDS)}}|>1$ means that $d_{e}>d_{e}^{\rm Org.}$
($d_{e}^{\rm HDS}>d_{e}^{\rm Org.,HDS}$), and hence, $d_{e}\sim d_{e}^{\rm
Ren.}$ ($d_{e}^{\rm HDS}\sim d_{e}^{\rm Ren.,HDS}=d_{e}^{\rm Ren.}$).

###  4.2 Fermion dark sector model

Let us consider a dark sector model which contains a vector-like singlet
fermion $\psi_{L,R}$ and an SM singlet complex scalar $S$. The dark sector
Lagrangian and the scalar potential include the following terms:

$\displaystyle\mathcal{L}_{\rm DS}$
$\displaystyle=-y_{\psi}\overline{\psi}_{L}\psi_{R}S-m_{\psi}\overline{\psi}_{L}\psi_{R}+{\rm
h.c.}\,,$ ( 4.30) $\displaystyle V_{\rm scl}$ $\displaystyle\supset
V_{H}+m_{S}^{2}|S|^{2}+\Bigl{[}\mu_{1S}H_{1}^{\dagger}H_{1}S+\mu_{2S}H_{2}^{\dagger}H_{2}S+\mu_{12S}H_{1}^{\dagger}H_{2}S+\mu_{21S}H_{2}^{\dagger}H_{1}S+{\rm
h.c.}\Bigr{]}\,,$ ( 4.31)

where $V_{H}$ represents the scalar potential of the 2HDM sector given in Eq.
(2.1).444 For the scalar potential, we have omitted some terms that are
allowed by the charge assignment in Table 6, because these terms are
irrelevant to our analysis for the electron EDM. The full scalar potential and
its analysis are shown in appendix B. We assume a $Z_{2}$ symmetry in which
$H_{1}$ and $\psi_{L,R}$ are odd. The charge assignment is summarized in Table
6. According to this $Z_{2}$ assignment, $\mu_{12S}$ and $\mu_{21S}$ terms are
soft $Z_{2}$ breaking terms. It should be noted that $\mu_{12S}$ and
$\mu_{21S}$ terms are necessary for the mixing among CP-odd components of
$H_{1,2}$ and $S$ (see appendix B for details). Similar to section 4.1, we
parameterize $S$ as $S=(s_{1}^{0}+is_{2}^{0})/\sqrt{2}$ and assume that $S$
does not acquire a nonzero VEV.

$H_{1}$ | $H_{2}$ | $S$ | $Q_{L}$ | $u_{R}$ | $d_{R}$ | $L_{L}$ | $e_{R}$ | $\psi_{L}$ | $\psi_{R}$
---|---|---|---|---|---|---|---|---|---
$-$ | $+$ | $+$ | $+$ | $+$ | $z_{2}^{d}$ | $+$ | $z_{2}^{e}$ | $-$ | $-$
Table 6: The $Z_{2}$ charge assignment for the fermion dark sector model where
$z_{2}^{d}$ and $z_{2}^{e}$ depend on the type of 2HDM.

The model incorporates numerous complex parameters, encompassing 13 invariant
phases as detailed in appendix C.2. Adhering to the presence of CP violation
solely within the dark sector, we select the following invariant phase as non-
zero, contributing to the electron EDM:

$\displaystyle\theta_{\rm phys}=\arg\left[\mu_{1S}^{*}y_{\psi}\right]\,,$ (
4.32)

while assuming the other phases to be zero. Consequently, if we choose the
phase of $S$ to have real $\mu_{1S}$, then the only CP violation happens
exclusively within the Yukawa interaction of the dark sector. It is pertinent
to note that a chiral phase rotation of the dark fermion $\psi_{L,R}$ has
already been applied to render the mass parameter $m_{\psi}$ real. Hence, we
choose $y_{\psi}$ to be complex and parameterize it by using two nonzero real
parameters, $y_{\psi}^{R}$ and $y_{\psi}^{I}$, as
$y_{\psi}=y_{\psi}^{R}+iy_{\psi}^{I}$. Defining
$\theta_{\psi}\equiv\arg(y_{\psi})$, we also have
$y_{\psi}^{R}=y_{\psi}\cos\theta_{\psi}$ and
$y_{\psi}^{I}=y_{\psi}\sin\theta_{\psi}$, and $\theta_{\psi}$ corresponds to
the physical CP phase. Consequently, we have explicit CP violation in the dark
sector Yukawa couplings:

$\displaystyle\mathcal{L}_{\not{\rm CP}}$
$\displaystyle=-\frac{\left(y_{\psi}^{R}s_{1}^{0}-y_{\psi}^{I}s_{2}^{0}\right)}{\sqrt{2}}\overline{\psi}\psi-\frac{\left(y_{\psi}^{I}s_{1}^{0}+y_{\psi}^{R}s_{2}^{0}\right)}{\sqrt{2}}\overline{\psi}i\gamma^{5}\psi-
m_{\psi}\overline{\psi}\psi\,.$ ( 4.33)

The CP-even scalars $(H_{1}^{0},H_{2}^{0},s_{1}^{0})$ and the CP-odd scalars
$(A_{1}^{0},A_{2}^{0},s_{2}^{0})$ are respectively expressed in terms of CP-
even mass eigenstates $\phi_{1,2,3}^{0}$ and CP-odd mass eigenstates
$\eta_{1,2,3}^{0}$ as

$\displaystyle H_{1}^{0}$
$\displaystyle=\sum_{j=1}^{3}(R_{E})_{1j}\phi_{j}^{0}\,,\quad
H_{2}^{0}=\sum_{j=1}^{3}(R_{E})_{2j}\phi_{j}^{0}\,,\quad
s_{1}^{0}=\sum_{j=1}^{3}(R_{E})_{3j}\phi_{j}^{0}\,,$ ( 4.34) $\displaystyle
A_{1}^{0}$ $\displaystyle=\sum_{a=1}^{3}(R_{O})_{1a}\eta_{a}^{0}\,,\quad
A_{2}^{0}=\sum_{a=1}^{3}(R_{O})_{2a}\eta_{a}^{0}\,,\quad\leavevmode\nobreak\
s_{2}^{0}=\sum_{a=1}^{3}(R_{O})_{3a}\eta_{a}^{0}\,.$ ( 4.35)

Here, we choose $\eta_{1}^{0}$ to be the neutral NG mode as mentioned in
section 3.2, and hence, we can parameterize $R_{O}$ by

$\displaystyle
R_{O}=\begin{pmatrix}c_{\beta}&-s_{\beta}c_{\theta_{O}}&s_{\beta}s_{\theta_{O}}\\\
s_{\beta}&c_{\beta}c_{\theta_{O}}&-c_{\beta}s_{\theta_{O}}\\\
0&s_{\theta_{O}}&c_{\theta_{O}}\end{pmatrix}\,,$ ( 4.36)

where $\theta_{O}$ is a mixing angle between $A^{0}$ in the 2HDM sector and
$s_{2}^{0}$. On the other hand, $R_{E}$ is parameterized by three mixing
angles $\alpha,\alpha_{hS},\alpha_{HS}$ as

$\displaystyle
R_{E}=\begin{pmatrix}c_{\alpha}c_{\alpha_{hS}}&-s_{\alpha}c_{\alpha_{HS}}-c_{\alpha}s_{\alpha_{hS}}s_{\alpha_{HS}}&-c_{\alpha}s_{\alpha_{hS}}c_{\alpha_{HS}}+s_{\alpha}s_{\alpha_{HS}}\\\
s_{\alpha}c_{\alpha_{hS}}&c_{\alpha}c_{\alpha_{HS}}-s_{\alpha}s_{\alpha_{hS}}s_{\alpha_{HS}}&-s_{\alpha}s_{\alpha_{hS}}c_{\alpha_{HS}}-c_{\alpha}s_{\alpha_{HS}}\\\
s_{\alpha_{hS}}&c_{\alpha_{hS}}s_{\alpha_{HS}}&c_{\alpha_{hS}}c_{\alpha_{HS}}\end{pmatrix}\,.$
( 4.37)

With these mixing angles, we can rewrite the explicit CP-violating couplings
in Eq. (4.33) in terms of the scalar mass eigenstates,

$\displaystyle\mathcal{L}_{\not{\rm CP}}^{\rm
Yuk}=-\frac{1}{\sqrt{2}}\sum_{j=1}^{3}\overline{\psi}\left(\hat{y}_{\psi,j}^{R}+i\hat{y}_{\psi,j}^{I}\gamma^{5}\right)\psi\phi_{j}^{0}+\frac{1}{\sqrt{2}}\sum_{a=2}^{3}\overline{\psi}\left(\tilde{y}_{\psi,a}^{I}-i\tilde{y}_{\psi,a}^{R}\gamma^{5}\right)\psi\eta_{a}^{0}\,,$
( 4.38)

where we have defined $\hat{y}_{\psi,j}^{R,I}\equiv
y_{\psi}^{R,I}(R_{E})_{3j}$ and $\tilde{y}_{\psi,a}^{R,I}\equiv
y_{\psi}^{R,I}(R_{O})_{3a}$. The SM fermion Yukawa couplings and the $W$ boson
couplings are also rewritten as

$\displaystyle\mathcal{L}_{\rm Yuk}$
$\displaystyle\supset-\sum_{f=u,d,e}\frac{m_{f}}{v_{\rm
SM}}\left[\sum_{j=1}^{3}\xi_{\phi_{j}}^{f}\bar{f}f\phi_{j}-\sum_{a=2}^{3}\xi_{\eta_{a}}^{f}\bar{f}i\gamma^{5}f\eta_{a}\right]\,,$
( 4.39) $\displaystyle\mathcal{L}$
$\displaystyle\supset\sum_{j=1}^{3}\frac{2m_{W}^{2}}{v_{\rm
SM}}\xi_{\phi_{j}}^{W}W_{\mu}^{+}W^{-\,\mu}\phi_{j}\,.$ ( 4.40)

Here, $\xi_{\phi_{j}}^{f}$ and $\xi_{\eta_{a}}^{f}$ are summarized in Table 7,
and $\xi_{\phi_{j}}^{W}\equiv c_{\beta}(R_{E})_{1j}+s_{\beta}(R_{E})_{2j}$.

| $\xi_{\phi_{j}}^{u}$ | $\xi_{\eta_{a}}^{u}$ | $\xi_{\phi_{j}}^{d}$ | $\xi_{\eta_{a}}^{d}$ | $\xi_{\phi_{j}}^{e}$ | $\xi_{\eta_{a}}^{e}$
---|---|---|---|---|---|---
Type-I | $(R_{E})_{2j}/s_{\beta}$ | $(R_{O})_{2a}\cot\beta$ | $(R_{E})_{2j}/s_{\beta}$ | $-(R_{O})_{2a}\cot\beta$ | $(R_{E})_{2j}/s_{\beta}$ | $-(R_{O})_{2a}\cot\beta$
Type-II | $(R_{E})_{2j}/s_{\beta}$ | $(R_{O})_{2a}\cot\beta$ | $(R_{E})_{1j}/c_{\beta}$ | $-(R_{O})_{1a}\tan\beta$ | $(R_{E})_{1j}/c_{\beta}$ | $-(R_{O})_{1a}\tan\beta$
Type-X | $(R_{E})_{2j}/s_{\beta}$ | $(R_{O})_{2a}\cot\beta$ | $(R_{E})_{2j}/s_{\beta}$ | $-(R_{O})_{2a}\cot\beta$ | $(R_{E})_{1j}/c_{\beta}$ | $-(R_{O})_{1a}\tan\beta$
Type-Y | $(R_{E})_{2j}/s_{\beta}$ | $(R_{O})_{2a}\cot\beta$ | $(R_{E})_{1j}/c_{\beta}$ | $-(R_{O})_{1a}\tan\beta$ | $(R_{E})_{2j}/s_{\beta}$ | $-(R_{O})_{2a}\cot\beta$

Table 7: The coefficients $\xi_{\phi_{j}}^{f}$ and $\xi_{\eta_{a}}^{f}$ for SM
fermion Yukawa couplings in the fermion dark sector model.

With new Yukawa couplings for $\psi$ in Eq. (4.38), the CP-violating
contribution to the CP-even scalar $\phi_{j}$ and CP-odd scalar $\eta_{a}$
mixing diagrams in Eq. (2.10) can be calculated at one-loop by running the
vector-like fermion,

$\displaystyle F_{\rm DS}^{\phi_{j}\eta_{a}}(k^{2})$
$\displaystyle=N_{c}^{\psi}\frac{\hat{y}_{\psi,j}^{R}\tilde{y}_{\psi,a}^{I}}{8\pi^{2}}m_{\psi}^{2}\left(\Lambda_{\rm
UV}^{2}-\int_{0}^{1}\\!dx\mathop{\rm
ln}\frac{m_{\psi}^{2}-x(1-x)k^{2}}{\mu^{2}}\right),$ ( 4.41)
$\displaystyle=N_{c}^{\psi}\frac{y_{\psi}^{R}y_{\psi}^{I}}{8\pi^{2}}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}\left(\Lambda_{\rm
UV}^{2}+2+{\rm DiscB}(k^{2},m_{\psi},m_{\psi})+\mathop{\rm
ln}\frac{\mu^{2}}{m_{\psi}^{2}}\right)\,,$

where $\Lambda_{\rm UV}$ denotes the UV divergent parameter, $\mu$ is the ’t
Hooft mass parameter and $N_{c}^{\psi}$ represents some new color factor for
the vector-like fermion $\psi$. In the present study, we set $N_{c}^{\psi}=1$.

Similar to the scalar dark sector model, we need to properly handle the UV
divergent part. In the current case, we consider the wave function
renormalization for the $\phi_{j}$-$\eta_{a}$ self-energy [50, 51]. Then, the
renormalized two-point function $\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})$ is
defined as

$\displaystyle\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})\equiv F_{\rm
DS}^{\phi_{j}\eta_{a}}(k^{2})+\frac{1}{2}(k^{2}-m_{\phi_{j}}^{2})\delta
Z_{\phi_{j}\eta_{a}}+\frac{1}{2}(k^{2}-m_{\eta_{a}}^{2})\delta
Z_{\eta_{a}\phi_{j}}\,,$ ( 4.42)

where no summations over $j$ and $a$ should be understood. In the on-shell
subtraction scheme, we can determine $\delta Z_{\phi_{j}\eta_{a}}$ and $\delta
Z_{\eta_{a}\phi_{j}}$ by

$\displaystyle{\rm
Re}\bigl{\\{}\hat{\Sigma}_{\phi_{j}\eta_{a}}(m_{\phi_{j}}^{2})\bigr{\\}}=0\,,\qquad{\rm
Re}\bigl{\\{}\hat{\Sigma}_{\phi_{j}\eta_{a}}(m_{\eta_{a}}^{2})\bigr{\\}}=0\,.$
( 4.43)

By solving these equations, we can obtain the renormalized two-point function
$\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})$ as

$\displaystyle\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})=\frac{y_{\psi}^{R}y_{\psi}^{I}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}}{8\pi^{2}(m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2})}$
$\displaystyle\Bigl{[}\Bigr{.}(m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2}){\rm
DiscB}(k^{2},m_{\psi},m_{\psi})$ $\displaystyle+(m_{\eta_{a}}^{2}-k^{2}){\rm
Re}\bigl{\\{}{\rm DiscB}(m_{\phi_{j}}^{2},m_{\psi},m_{\psi})\bigr{\\}}$
$\displaystyle+(k^{2}-m_{\phi_{j}}^{2}){\rm Re}\bigl{\\{}{\rm
DiscB}(m_{\eta_{a}}^{2},m_{\psi},m_{\psi})\bigr{\\}}\Bigl{.}\Bigr{]}\,,$ (
4.44)

where the divergent part is properly canceled.

In the HDS approximation, we can further expand
$\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})$ as

$\displaystyle\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})=a_{0}^{\phi_{j}\eta_{a}}+a_{1}^{\phi_{j}\eta_{a}}\frac{k^{2}}{m_{\psi}^{2}}+\mathcal{O}\left(\frac{k^{4}}{m_{\psi}^{4}}\right)\,.$
( 4.45)

Here, $M^{2}$ in Eq. (3.16) corresponds to the vector-like fermion mass
$m_{\psi}^{2}$, and the coefficients are

$\displaystyle a_{0}^{\phi_{j}\eta_{a}}(m_{\psi})$
$\displaystyle=\frac{y_{\psi}^{R}y_{\psi}^{I}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}}{8\pi^{2}(m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2})}\left[-2(m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2})-m_{\phi_{j}}^{2}{\rm
Re}\bigl{\\{}{\rm DiscB}(m_{\eta_{a}}^{2},m_{\psi},m_{\psi})\bigr{\\}}\right.$
$\displaystyle\hskip 130.0002pt\left.+m_{\eta_{a}}^{2}{\rm Re}\bigl{\\{}{\rm
DiscB}(m_{\phi_{j}}^{2},m_{\psi},m_{\psi})\bigr{\\}}\right]\,,$ ( 4.46)
$\displaystyle a_{1}^{\phi_{j}\eta_{a}}(m_{\psi})$
$\displaystyle=\frac{y_{\psi}^{R}y_{\psi}^{I}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}}{8\pi^{2}(m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2})}\Biggl{[}\Biggr{.}\frac{m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2}}{6}-m_{\psi}^{2}{\rm
Re}\bigl{\\{}{\rm DiscB}(m_{\phi_{j}}^{2},m_{\psi},m_{\psi})\bigr{\\}}$
$\displaystyle\hskip 130.0002pt+m_{\psi}^{2}{\rm Re}\bigl{\\{}{\rm
DiscB}(m_{\eta_{a}}^{2},m_{\psi},m_{\psi})\bigr{\\}}\Biggl{.}\Biggr{]}\,.$ (
4.47)

As in the case of the scalar dark sector model,
$a_{0}^{\phi_{j}\eta_{a}}(m_{\psi})$ is used for the calculation of the four-
fermion operator coefficients $C_{ff^{\prime}}$.

####  4.2.1 Numerical results

Similar to the scalar dark sector model, to numerically estimate the size of
the electron EDM generated in the current model, we choose the following
masses and parameters:

$\displaystyle m_{\phi_{1}}=125.25\leavevmode\nobreak\ {\rm
GeV},\leavevmode\nobreak\ m_{\phi_{2}}=1500\leavevmode\nobreak\ {\rm
GeV},\leavevmode\nobreak\ m_{\phi_{3}}=1600\leavevmode\nobreak\ {\rm
GeV},\leavevmode\nobreak\ \alpha=\frac{\pi}{3},\leavevmode\nobreak\
\alpha_{hS}=0.01,\leavevmode\nobreak\ \alpha_{HS}=\frac{\pi}{6},$
$\displaystyle m_{\eta_{2}}=1550\leavevmode\nobreak\ {\rm
GeV},\leavevmode\nobreak\ m_{\eta_{3}}=1650\leavevmode\nobreak\ {\rm
GeV},\leavevmode\nobreak\ \tan\beta=50,\leavevmode\nobreak\
\theta_{O}=\frac{\pi}{4},\leavevmode\nobreak\ \leavevmode\nobreak\
y_{\psi}^{R}=y_{\psi}^{I}=\frac{1}{\sqrt{2}}.$ ( 4.48)

The choices of $y_{\psi}^{R,I}$ can be obtained by $|y_{\psi}|=1$ and
$\theta_{\psi}=\pi/4$, and the value of $\theta_{\psi}$ corresponds to the one
which maximizes the $d_{e}$ with the assumption where all parameters in
$V_{\rm scl}$ are real. We then calculate $|d_{e}|$ as a function of
$m_{\psi}$. The $C_{S}$ contribution is found to be significantly smaller than
the electron EDM contribution $d_{e}$, and hence neglected in the subsequent
calculations.

Figure 7: Numerical results of the electron EDM in the fermion dark sector
model. The left panel shows the cases of Type-I (dashed line) and Type-II
(solid line) 2HDMs, while the right panel shows those of Type-X (solid line)
and Type-Y (dashed line) 2HDMs. The results of the HDS approximation for Type-
II and Type-X 2HDMs are also shown by dotted lines in each panel. The gray-
shaded region is excluded by the current electron EDM bound from the JILA
experiment [18]. In this plot, we choose the parameters that satisfy the
constraint of the SM Higgs invisible decay, $h^{0}\to\psi\bar{\psi}$.

The left panel of Fig. 7 presents curves for $|d_{e}|$ in Type-I (dashed line)
and Type-II (solid line) 2HDMs, and the right panel shows those for Type-X
(solid line) and Type-Y (dashed line) 2HDMs. To get those results, we have
used the general expression in Eq. (3.13) by omitting $k\cdot q$ and $k\cdot
p_{1,2}$ and calculated contributions from the 3rd generation fermions ($t,b$)
and the $W$ boson in the inner loop, as in the case of the scalar dark sector
model. The gray-shaded region is excluded by the constraint from the JILA
experiments [18]. We may have a constraint from the SM Higgs invisible decay,
$h^{0}\to\psi\bar{\psi}$, especially if the vector-like fermion is light. The
branching ratio is directly proportional to
$|y_{\psi}|^{2}\sin^{2}\alpha_{hS}$, making $|y_{\psi}|$ and $\alpha_{hS}$
critical parameters. Larger values of $|y_{\psi}|$ and/or $\alpha_{hS}$ tend
to exclude the lower mass region $m_{\psi}<m_{h}/2$. We have checked that the
parameter choice in Eq. (4.48) can avoid the constraint.

The curves of the electron EDM $|d_{e}|$ in Fig. 7 show significant
complexity, stemming from the $m_{\psi}$ dependence of the renormalized two-
point function as given in Eq. (4.44). Notably, the region that $|d_{e}|$ is
peaked corresponds to $m_{\psi}\simeq m_{\phi_{j}}/2$ or $m_{\eta_{a}}/2$,
where $\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})$ is enhanced. As in the case of
the scalar dark sector model, Type-II and Type-X 2HDMs can reach the current
upper limit in this mass region, while Type-I and Type-Y 2HDMs predict
$|d_{e}|$ that is several orders of magnitude smaller than the limit.
Moreover, $\hat{\Sigma}_{\phi_{j}\eta_{a}}(k^{2})$ approaches to zero as
$m_{\psi}$ goes to zero or infinity. Such a feature is visible in Fig. 7. Note
that for $m_{\psi}\gtrsim 7\times 10^{3}$ GeV, the curves of $|d_{e}|$ change
their behaviors, due to the low accuracy in the numerical integration, which
is explained below.

####  4.2.2 Discussion on the HDS approximation

Remarkably, the HDS approximation results which are represented by the dotted
curves for the Type-II (left panel) and Type-X (right panel) 2HDMs in Fig. 7
align well with the results without using the approximation. Note that the
result for the HDS approximation can be obtained from Eq. (3.17) with Eqs.
(4.46) and (4.47).

Compared with the scalar dark sector model, the HDS approximation is good, due
to the simple expression for $\hat{\Sigma}_{\phi_{j}\eta_{a}}$ in Eq. (4.44)
(see Eqs. (4.16) and (4.17) together with Eqs. (4.9) and (4.10) for the scalar
dark sector model). Hence, we may be able to use the HDS approximation results
to estimate the $d_{e}$ prediction in the heavy $m_{\psi}$ limit.

Figure 8: Plots corresponding to Fig. 6 for the fermion dark sector model.

Fig. 8 shows $|1-\rho_{d_{e}}^{\rm Org.}|$ and $|\epsilon_{d_{e}}|$ which are
defined in Eqs. (4.28) and (4.29), and the expressions for $d_{e}^{\rm Org.}$
and $d_{e}^{\rm Ren.}$ are summarized in appendix D.3. From the left panel of
Fig. 8, we can see that the HDS approximation does not give the correct
results for a light $m_{\psi}$, as observed from Fig. 7. On the other hand,
when $m_{\psi}>1000$ GeV, the accuracy of $d_{e}^{\rm Org.}$ becomes less than
1%, and hence, the HDS approximation can be used to estimate the value of the
electron EDM. However, for $m_{\psi}\gtrsim 7\times 10^{3}$ GeV, the accuracy
of $d_{e}^{\rm Org.}$ is saturated, and the curves with and without the HDS
approximation deviate from each other (see Fig. 7). For this mass range, we
need a more accurate calculation for numerical integration, which is shown in
the right panel of Fig. 8. For $m_{\psi}\sim 7\gtrsim 10^{3}$ GeV, the
cancellation between $d_{e}^{\rm Org.}$ and $d_{e}^{\rm Ren.}$ is also
saturated, at the level of $|\epsilon_{d_{e}}|\simeq 10^{-5}$, while that
between $d_{e}^{\rm Org.,HDS}$ and $d_{e}^{\rm Ren.,HDS}$ decreases
monotonically. Therefore, even if we obtain a good accuracy for
$|1-\rho_{d_{e}}^{\rm Org.}|$, the numerical results with and without the HDS
approximation deviate for the mass range of $m_{\psi}\gtrsim 7\times 10^{3}$
GeV. As a result, for such a heavy $m_{\psi}$, the correct prediction can be
obtained from the HDS approximation, although the size is far below the
current upper limit on $|d_{e}|$.

Around $m_{\psi}\simeq m_{\phi_{j}}/2$ and $m_{\eta_{a}}/2$, it seems that the
HDS approximation is quite good in Fig. 7. The reason is the same as that of
the scalar dark sector model: in this mass region, $d_{e}$ ($d_{e}^{\rm HDS}$)
is dominated by $d_{e}^{\rm Ren.}$ ($d_{e}^{\rm Ren.,HDS}$), as one can see
from the right panel of Fig. 8.

##  5 Conclusion

In the present paper, we have discussed the possibility that a CP-violating
dark sector coupling to the 2HDM generates a detectable electron EDM. The 2HDM
was assumed to preserve the CP symmetry at the beginning, otherwise, the model
is severely constrained, while the dark sector contains CP violation which is
mediated to the 2HDM sector by radiative corrections. We have presented a
general form of the electron EDM induced by BZ diagrams shown in Fig. 2, which
can be applied to a wide range of dark sector models with CP violation.

For benchmark scenarios, we considered the scalar dark sector model and the
fermion dark sector model. In the scalar dark sector model, CP violation
arises in the scalar potential involving neutral Higgs and dark sector
scalars. In the fermion dark sector model, a singlet complex scalar serves as
a portal field, connecting a vector-like singlet fermion to the 2HDM through a
complex CP-violating Yukawa coupling. For both scenarios, we carefully
considered the two-point renormalized scalar self-energy functions to ensure
finite results under the on-shell subtraction scheme. In each model, we found
that Type-II and Type-X 2HDM cases exhibit similar predictions for $|d_{e}|$
and have a chance to reach the current upper limit on the electron EDM.
Therefore, we can constrain the model parameter space: the mass scale of the
dark sector and/or the size of effective mixing for CP-even and CP-odd scalars
in the 2HDMs. The analytical expression of the electron EDM was given in the
Heavy Dark Sector limit, and we have checked that it is valid to use the
approximation for the estimation of the size of the electron EDM.

Since the experimental sensitivity will be improved in the future [53], the
discovery of the electron EDM or its more stringent limit is expected to be
obtained. In either case, our results will give important implications for
physics beyond the SM with CP violation.

## Acknowledgements

We would like to thank Michael Ramsey-Musolf for very helpful discussions. YN
is supported by Natural Science Foundation of China under grant No.
12150610465. JL is supported by Natural Science Foundation of China under
grant No. 12075005 and 12235001.

##  Appendix A The 2HDM

Here, we summarize mass eigenvalues and eigenstates of the 2HDM. The Higgs
potential in Eq. (2.1) leads to $2\times 2$ mass matrices for CP-even and CP-
odd neutral scalars as

$\displaystyle M_{\rm even}^{2}$
$\displaystyle=\begin{pmatrix}m_{12}^{2}\tan\beta+\lambda_{1}v_{1}^{2}&-m_{12}^{2}+\lambda_{345}v_{1}v_{2}\\\
-m_{12}^{2}+\lambda_{345}v_{1}v_{2}&m_{12}^{2}\cot\beta+\lambda_{2}v_{2}^{2}\end{pmatrix}\,,$
( A.1) $\displaystyle M_{\rm odd}^{2}$
$\displaystyle=\begin{pmatrix}m_{12}^{2}\tan\beta-\lambda_{5}v_{2}^{2}&-m_{12}^{2}+\lambda_{5}v_{1}v_{2}\\\
-m_{12}^{2}+\lambda_{5}v_{1}v_{2}&m_{12}^{2}\cot\beta-\lambda_{5}v_{1}^{2}\end{pmatrix}\,,$
( A.2)

in basis of $(H_{1}^{0},H_{2}^{0})$ and $(A_{1}^{0},A_{2}^{0})$, respectively.
For charged scalars, we find

$\displaystyle
M_{\pm}^{2}=\begin{pmatrix}m_{12}^{2}\tan\beta-\frac{\lambda_{4}+\lambda_{5}}{2}v_{2}^{2}&-m_{12}^{2}+\frac{\lambda_{4}+\lambda_{5}}{2}v_{1}v_{2}\\\
-m_{12}^{2}+\frac{\lambda_{4}+\lambda_{5}}{2}v_{1}v_{2}&m_{12}^{2}\cot\beta-\frac{\lambda_{4}+\lambda_{5}}{2}v_{1}^{2}\end{pmatrix}\,,$
( A.3)

in basis of $(H_{1}^{+},H_{2}^{+})$. Here, $m_{12}^{2}$ and $\lambda_{5}$ are
considered as real parameters, and we define
$\lambda_{345}\equiv\lambda_{3}+\lambda_{4}+\lambda_{5}$. The minimization
conditions for the Higgs potential impose relations among model parameters,

$\displaystyle m_{1}^{2}$
$\displaystyle=m_{12}^{2}\tan\beta-\frac{1}{2}\left(\lambda_{1}v_{1}^{2}+\lambda_{345}v_{2}^{2}\right)\,,$
( A.4) $\displaystyle m_{2}^{2}$
$\displaystyle=m_{12}^{2}\cot\beta-\frac{1}{2}\left(\lambda_{2}v_{2}^{2}+\lambda_{345}v_{1}^{2}\right)\,.$
( A.5)

The above mass matrices are diagonalized by orthogonal matrices, and each mass
eigenvalue can be obtained as

$\displaystyle m_{h}^{2}$ $\displaystyle=\frac{v_{\rm
SM}^{2}}{2}\left[\left(\lambda_{1}c_{\beta}^{2}+\lambda_{2}s_{\beta}^{2}+\lambda_{\beta}\right)-\sqrt{(\lambda_{1}c_{\beta}^{2}-\lambda_{2}s_{\beta}^{2}-\lambda_{\beta}c_{2\beta})^{2}+(\lambda_{345}-\lambda_{\beta})^{2}s_{2\beta}^{2}}\right]\,,$
( A.6) $\displaystyle m_{H}^{2}$ $\displaystyle=\frac{v_{\rm
SM}^{2}}{2}\left[\left(\lambda_{1}c_{\beta}^{2}+\lambda_{2}s_{\beta}^{2}+\lambda_{\beta}\right)+\sqrt{(\lambda_{1}c_{\beta}^{2}-\lambda_{2}s_{\beta}^{2}-\lambda_{\beta}c_{2\beta})^{2}+(\lambda_{345}-\lambda_{\beta})^{2}s_{2\beta}^{2}}\right]\,,$
( A.7) $\displaystyle m_{A}^{2}$
$\displaystyle=\left(\lambda_{\beta}-\lambda_{5}\right)v_{\rm SM}^{2}\,,$ (
A.8) $\displaystyle m_{H^{\pm}}^{2}$
$\displaystyle=\left(\lambda_{\beta}-\frac{1}{2}\lambda_{4}-\frac{1}{2}\lambda_{5}\right)v_{\rm
SM}^{2}=m_{A}^{2}-\frac{1}{2}\left(\lambda_{4}-\lambda_{5}\right)v_{\rm
SM}^{2}\,,$ ( A.9)

where $m_{h}$, $m_{H}$ and $m_{A}$ are the masses of the SM-like Higgs, heavy
CP-even and CP-odd scalars, respectively, and $m_{H^{\pm}}$ is the physical
charged scalar mass. Here, we have defined

$\displaystyle\lambda_{\beta}\equiv\frac{m_{12}^{2}}{v_{1}v_{2}}=\frac{1}{s_{\beta}c_{\beta}}\frac{m_{12}^{2}}{v_{\rm
SM}^{2}}\,.$ ( A.10)

Note that each of $M_{\pm}^{2}$ and $M_{\rm odd}^{2}$ has one zero eigenvalue,
which corresponds to NG modes, $G^{\pm}$ and $G^{0}$. The mass matrix for CP-
even neutral scalars $M_{\rm even}^{2}$ is diagonalized by one mixing angle
$\alpha$,

$\displaystyle\sin
2\alpha=\frac{(\lambda_{345}-\lambda_{\beta})s_{2\beta}}{\sqrt{(\lambda_{1}c_{\beta}^{2}-\lambda_{2}s_{\beta}^{2}-\lambda_{\beta}c_{2\beta})^{2}+(\lambda_{345}-\lambda_{\beta})^{2}s_{2\beta}^{2}}}\,,$
( A.11)

while $M_{\pm}^{2}$ and $M_{\rm odd}^{2}$ are diagonalized by using
$\tan\beta$ defined in Eq. (2.3). The mass eigenstates are obtained as in Eqs.
(2.4) and (2.5), and for charged scalars,

$\displaystyle H_{1}^{\pm}$
$\displaystyle=G^{\pm}c_{\beta}-H^{\pm}s_{\beta}\,,\qquad
H_{2}^{\pm}=G^{\pm}s_{\beta}+H^{\pm}c_{\beta}\,,$ ( A.12)

where $H^{\pm}$ is a physical charged scalar.

##  Appendix B Potential analysis for the fermion dark sector model

The scalar potential of the fermion dark sector model allowed by the $Z_{2}$
charge assignments in Table 6 is given by

$\displaystyle V_{\rm scl}$ $\displaystyle=V_{H}+V_{S}+V_{HS}\,,$ ( B.1)
$\displaystyle V_{S}$
$\displaystyle=m_{S}^{2}|S|^{2}+\frac{\lambda_{S}}{2}|S|^{4}+\Bigl{[}\mu_{2}^{2}S^{2}+\mu_{3}S|S|^{2}+\mu^{\prime}_{3}S^{3}+\lambda^{\prime}_{S}S^{2}|S|^{2}+\lambda^{\prime\prime}_{S}S^{4}+{\rm
h.c.}\Bigr{]}\,,$ ( B.2) $\displaystyle V_{HS}$
$\displaystyle=\Bigl{[}\mu_{1S}H_{1}^{\dagger}H_{1}S+\mu_{2S}H_{2}^{\dagger}H_{2}S+\mu_{12S}H_{1}^{\dagger}H_{2}S+\mu_{21S}H_{2}^{\dagger}H_{1}S+{\rm
h.c.}\Bigr{]}$ $\displaystyle\hskip
11.99998pt+\lambda_{1S}H_{1}^{\dagger}H_{1}|S|^{2}+\lambda_{2S}H_{2}^{\dagger}H_{2}|S|^{2}+\Bigl{[}\lambda^{\prime}_{1S}H_{1}^{\dagger}H_{1}S^{2}+\lambda^{\prime}_{2S}H_{2}^{\dagger}H_{2}S^{2}+{\rm
h.c.}\Bigr{]}\,.$ ( B.3)

We assume that all the parameters in $V_{\rm scl}$ are real, and $S$ does not
acquire a nonzero VEV. In the potential $V_{S}$, we can introduce a tadpole
term $T_{S}S+\text{h.c.}$ since $S$ is a singlet under the considered
symmetries. Nevertheless, this term is removable through a shift in $S$ by
appropriately redefining all the relevant couplings. The potential leads to
the following minimization condition:

$\displaystyle\mu_{1S}v_{1}^{2}+\mu_{2S}v_{2}^{2}+\left(\mu_{12S}+\mu_{21S}\right)v_{1}v_{2}=0\,,$
( B.4)

with others shown in Eqs. (A.4) and (A.5). After applying all the minimization
conditions, we find mass matrices for CP-even and CP-odd neutral scalars as

$\displaystyle M_{\rm even}^{2}$
$\displaystyle=\begin{pmatrix}\lambda_{\beta}v_{2}^{2}+\lambda_{1}v_{1}^{2}&\left(-\lambda_{\beta}+\lambda_{345}\right)v_{1}v_{2}&\hat{m}_{1S}^{2}\\\
\left(-\lambda_{\beta}+\lambda_{345}\right)v_{1}v_{2}&\lambda_{\beta}v_{1}^{2}+\lambda_{2}v_{2}^{2}&\hat{m}_{2S}^{2}\\\
\hat{m}_{1S}^{2}&\hat{m}_{2S}^{2}&m_{S,+}^{2}\end{pmatrix}\,,$ ( B.5)
$\displaystyle M_{\rm odd}^{2}$
$\displaystyle=\begin{pmatrix}\left(\lambda_{\beta}-\lambda_{5}\right)v_{2}^{2}&\left(-\lambda_{\beta}+\lambda_{5}\right)v_{1}v_{2}&\tilde{m}_{1S}^{2}\\\
\left(-\lambda_{\beta}+\lambda_{5}\right)v_{1}v_{2}&\left(\lambda_{\beta}-\lambda_{5}\right)v_{1}^{2}&\tilde{m}_{2S}^{2}\\\
\tilde{m}_{1S}^{2}&\tilde{m}_{2S}^{2}&m_{S,-}^{2}\end{pmatrix}\,,$ ( B.6)

where we define

$\displaystyle\hat{m}_{1S}^{2}$
$\displaystyle\equiv-\frac{\sqrt{2}v_{2}^{2}}{v_{1}}\mu_{2S}-\frac{\mu_{12S}+\mu_{21S}}{\sqrt{2}}v_{2}\,,\quad\hat{m}_{2S}^{2}\equiv\sqrt{2}\mu_{2S}v_{2}+\frac{\mu_{12S}+\mu_{21S}}{\sqrt{2}}v_{1}\,,$
( B.7) $\displaystyle\tilde{m}_{1S}^{2}$
$\displaystyle\equiv\frac{\mu_{12S}-\mu_{21S}}{\sqrt{2}}v_{2}\,,\quad\tilde{m}_{2S}^{2}\equiv-\frac{\mu_{12S}-\mu_{21S}}{\sqrt{2}}v_{1}\,,$
( B.8) $\displaystyle m_{S,\pm}^{2}$ $\displaystyle\equiv
m_{S}^{2}+\frac{\lambda_{1S}}{2}v_{1}^{2}+\frac{\lambda_{2S}}{2}v_{2}^{2}\pm\left(\lambda^{\prime}_{1S}v_{1}^{2}+\lambda^{\prime}_{2S}v_{2}^{2}+2\mu_{2}^{2}\right)\,,$
( B.9)

and $\lambda_{\beta}$ is defined in Eq. (A.10). Note that terms in the first
line of Eq. (B.3) are essential for mixings between neutral scalars in the
2HDM and $S$. In particular, if we do not introduce soft $Z_{2}$ breaking
terms $\mu_{12S}$ and $\mu_{21S}$, $\tilde{m}_{1S}^{2}=0$ and
$\tilde{m}_{2S}^{2}=0$, which results in no mixing between CP-odd components
of the 2HDM and $S$. Since $S$ is the neutral scalar, the charged scalar mass
is the same as that of the 2HDM. See appendix A.

Upon careful examination, it’s apparent that the determinant of $M_{\rm
odd}^{2}$ is zero. Consequently, one of the eigenvalues derived from the
matrix corresponds to the neutral NG boson field. By defining one mixing
matrix for $M_{\rm odd}^{2}$ as

$\displaystyle R_{O,1}=\begin{pmatrix}c_{\beta}&-s_{\beta}&0\\\
s_{\beta}&c_{\beta}&0\\\ 0&0&1\end{pmatrix}\,,$ ( B.10)

we obtain

$\displaystyle R_{O,1}^{T}M_{\rm odd}^{2}R_{O,1}=\begin{pmatrix}0&0&0\\\
0&\left(\lambda_{\beta}-\lambda_{5}\right)v_{\rm
SM}^{2}&-\tilde{m}_{1S}^{2}s_{\beta}+\tilde{m}_{2S}^{2}c_{\beta}\\\
0&-\tilde{m}_{1S}^{2}s_{\beta}+\tilde{m}_{2S}^{2}c_{\beta}&m_{S,-}^{2}\end{pmatrix}\,.$
( B.11)

The remaining part is diagonalized by another mixing matrix,

$\displaystyle R_{O,2}=\begin{pmatrix}1&0&0\\\
0&c_{\theta_{O}}&-s_{\theta_{O}}\\\
0&s_{\theta_{O}}&c_{\theta_{O}}\end{pmatrix}\,,$ ( B.12)

with the mixing angle,

$\displaystyle\sin 2\theta_{O}=\frac{2}{\Delta
m_{\eta}^{2}}\left(\tilde{m}_{1S}^{2}s_{\beta}-\tilde{m}_{2S}^{2}c_{\beta}\right)\,.$
( B.13)

Here, $\Delta m_{\eta}^{2}$ is the mass-squared difference of two physical CP-
odd scalars,

$\displaystyle\Delta m_{\eta}^{2}$
$\displaystyle=m_{\eta_{3}}^{2}-m_{\eta_{2}}^{2}$ ( B.14)
$\displaystyle=\sqrt{\left\\{\left(\lambda_{\beta}-\lambda_{5}\right)v_{\rm
SM}^{2}-m_{S,-}^{2}\right\\}^{2}+4\left(-\tilde{m}_{1S}^{2}s_{\beta}+\tilde{m}_{2S}^{2}c_{\beta}\right)^{2}}\,.$
( B.15)

Then, each mass eigenstate can be obtained as

$\displaystyle{\rm
diag}\left(0,m_{\eta_{2}}^{2},m_{\eta_{3}}^{2}\right)=R_{O,2}^{T}R_{O,1}^{T}M_{\rm
odd}^{2}R_{O,1}R_{O,2}=R_{O}^{T}M_{\rm odd}^{2}R_{O}\,,$ ( B.16)
$\displaystyle
m_{\eta_{2}}^{2}=\frac{1}{2}\left[\left(\lambda_{\beta}-\lambda_{5}\right)v_{\rm
SM}^{2}+m_{S,-}^{2}-\Delta m_{\eta}^{2}\right]\,,$ ( B.17) $\displaystyle
m_{\eta_{3}}^{2}=\frac{1}{2}\left[\left(\lambda_{\beta}-\lambda_{5}\right)v_{\rm
SM}^{2}+m_{S,-}^{2}+\Delta m_{\eta}^{2}\right]\,,$ ( B.18)

where $R_{O}$ is in Eq. (4.36).

In contrast to the CP-odd sector, the CP-even sector requires three mixing
angles for diagonalization of $M_{\rm even}^{2}$, as shown in Eq. (4.37). The
mixing matrix $R_{E}$ is obtained as

$\displaystyle R_{E}$ $\displaystyle=R_{E,1}R_{E,2}R_{E,3}\,,$ ( B.19)
$\displaystyle R_{E,1}$
$\displaystyle=\begin{pmatrix}c_{\alpha}&-s_{\alpha}&0\\\
s_{\alpha}&c_{\alpha}&0\\\ 0&0&1\end{pmatrix}\,,\leavevmode\nobreak\
R_{E,2}=\begin{pmatrix}c_{\alpha_{hS}}&0&-s_{\alpha_{hS}}\\\ 0&1&0\\\
s_{\alpha_{hS}}&0&c_{\alpha_{hS}}\end{pmatrix}\,,\leavevmode\nobreak\
R_{E,3}=\begin{pmatrix}1&0&0\\\ 0&c_{\alpha_{HS}}&-s_{\alpha_{HS}}\\\
0&s_{\alpha_{HS}}&c_{\alpha_{HS}}\end{pmatrix}\,.$ ( B.20)

##  Appendix C Rephasing-invariant CP-violating phases

We here summarize rephasing-invariant CP-violating (CPV) phases for the scalar
dark sector model and the fermion dark sector model we focused on in the
present paper. Detailed studies have been done in refs. [48, 49].

###  C.1 CPV phases in the scalar dark sector model

In the scalar potential, we have six complex parameters:
$m_{12}^{2},\lambda_{5},B_{S},D_{S},\lambda_{12S},\lambda_{21S}$, and one of
VEVs of $H_{1,2}$ is generally complex, e.g., $v_{1}=v_{1}^{*}$ and
$v_{2}=|v_{2}|e^{i\xi}$. Then, by considering rephases for all complex scalars
as

$\displaystyle H_{1}\to e^{i\theta_{1}}H_{1}\,,\qquad H_{2}\to
e^{i\theta_{2}}H_{2}\,,\qquad S\to e^{i\theta_{S}}S\,,$ ( C.1)

phases of $\theta_{2}-\theta_{1}$ and $\theta_{S}$ can be absorbed by
redefining all complex parameters and VEVs as

$\displaystyle m_{12}^{2}\to
m_{12}^{2}e^{i(\theta_{1}-\theta_{2})}\,,\quad\lambda_{5}\to\lambda_{5}e^{2i(\theta_{1}-\theta_{2})}\,,\quad
v_{1}v_{2}^{*}\to v_{1}v_{2}^{*}e^{i(\theta_{1}-\theta_{2})}\,,$ ( C.2)
$\displaystyle B_{S}\to B_{S}e^{-2i\theta_{S}}\,,\quad D_{S}\to
D_{S}e^{-4i\theta_{S}}\,,$ ( C.3)
$\displaystyle\lambda_{12S}\to\lambda_{12S}e^{i(\theta_{1}-\theta_{2}-2\theta_{S})}\,,\quad\lambda_{21S}\to\lambda_{21S}e^{i(\theta_{2}-\theta_{1}-2\theta_{S})}\,.$
( C.4)

As a result, the scalar potential in Eqs. (2.1), (4.1) and (4.2) is invariant
under rephasing of Eq. (C.1). Then, according to the discussion in refs. [48,
49], we have the following five rephasing-invariant CPV phases:

$\displaystyle\delta_{\lambda_{5}}$
$\displaystyle=\arg\left[\lambda_{5}^{*}(v_{1}v_{2}^{*})^{2}\right],\quad\delta_{m_{12}^{2}}=\arg\left[(m_{12}^{2})^{*}v_{1}v_{2}^{*}\right],\quad\delta_{S}=\arg\left[B_{S}^{2}D_{S}^{*}\right],$
( C.5) $\displaystyle\delta_{\lambda_{12S}}$
$\displaystyle=\arg\left[\lambda_{12S}^{*}\lambda_{21S}(v_{1}v_{2}^{*})^{2}\right],\quad\delta_{\lambda_{21S}}=\arg\left[B_{S}^{2}\lambda_{12S}^{*}\lambda_{21S}^{*}\right].$
( C.6)

Note that by selecting a CP-conserving scenario within the 2HDM sector, where
$m^{2}_{12}=0$, the minimization of the scalar potential implies that
$\text{Im}\lambda_{5}=0$. In the case where the CPV phase only shows up in the
dark sector, one could opt to assign real values to $\lambda_{12S}$ and
$\lambda_{21S}$, resulting in no invariant phases within the portal sector.
Consequently, the sole source of CP violation would be $\delta_{S}$, with
$D_{S}$ assumed to be real and only $B_{S}$ is complex. This configuration
leads to the remaining physical CP phase being expressed as Eq. (4.3).

###  C.2 CPV phases in the fermion dark sector model

In the relevant Lagrangian of the fermion dark sector model, we have 15
complex parameters:
$y_{\psi},m_{\psi},m_{12}^{2},\lambda_{5},\mu_{2}^{2},\mu_{3},\mu^{\prime}_{3},\mu_{1S},\mu_{2S},\mu_{12S},\mu_{21S},\lambda^{\prime}_{S},\lambda^{\prime\prime}_{S},\lambda^{\prime}_{1S},\lambda^{\prime}_{2S}$.
Among them, $m_{\psi}$ can be taken to be real by using the chiral phase
rotation of $\psi_{L,R}$, as mentioned in the main text. As in the case of the
scalar dark sector model, rephases for all complex scalars in Eq. (C.1) can be
absorbed by redefining complex parameters as

$\displaystyle y_{\psi}\to y_{\psi}e^{-i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ m_{12}^{2}\to
m_{12}^{2}e^{i(\theta_{1}-\theta_{2})}\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\lambda_{5}\to\lambda_{5}e^{2i(\theta_{1}-\theta_{2})}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ v_{1}v_{2}^{*}\to
v_{1}v_{2}^{*}e^{i(\theta_{1}-\theta_{2})}\,,$ ( C.7)
$\displaystyle\mu_{2}^{2}\to\mu_{2}^{2}e^{-2i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \mu_{3}\to\mu_{3}e^{-i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \mu^{\prime}_{3}\to\mu^{\prime}_{3}e^{-3i\theta_{S}}\,,$
( C.8)
$\displaystyle\mu_{1S}\to\mu_{1S}e^{-i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \mu_{2S}\to\mu_{2S}e^{-i\theta_{S}}\,,$ ( C.9)
$\displaystyle\mu_{12S}\to\mu_{12S}e^{i(\theta_{1}-\theta_{2}-\theta_{S})}\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\mu_{21S}\to\mu_{21S}e^{i(\theta_{2}-\theta_{1}-\theta_{S})}\,,$ ( C.10)
$\displaystyle\lambda^{\prime}_{S}\to\lambda^{\prime}_{S}e^{-2i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\lambda^{\prime\prime}_{S}\to\lambda^{\prime\prime}_{S}e^{-4i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\lambda^{\prime}_{1S}\to\lambda^{\prime}_{1S}e^{-2i\theta_{S}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\lambda^{\prime}_{2S}\to\lambda^{\prime}_{2S}e^{-2i\theta_{S}}\,.$ ( C.11)

Then, we find 13 combinations for rephasing-invariant CPV phases as

$\displaystyle\delta_{\lambda_{5}}=\arg\left[\lambda_{5}^{*}(v_{1}v_{2}^{*})^{2}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{m_{12}^{2}}=\arg\left[(m_{12}^{2})^{*}v_{1}v_{2}^{*}\right]$ ( C.12)
$\displaystyle\delta_{\mu_{2}^{2}}=\arg\left[(\mu_{2}^{2})^{*}y_{\psi}^{2}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\mu_{3}}=\arg\left[\mu_{3}^{*}y_{\psi}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\mu^{\prime}_{3}}=\arg\left[(\mu^{\prime}_{3})^{*}y_{\psi}^{3}\right]\,,$
( C.13)
$\displaystyle\delta_{\lambda^{\prime}_{S}}=\arg\left[(\lambda^{\prime}_{S})^{*}y_{\psi}^{2}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\lambda^{\prime\prime}_{S}}=\arg\left[(\lambda^{\prime\prime}_{S})^{*}y_{\psi}^{4}\right]\,,$
( C.14)
$\displaystyle\delta_{\mu_{1S}}=\arg\left[\mu_{1S}^{*}y_{\psi}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\mu_{2S}}=\arg\left[\mu_{2S}^{*}y_{\psi}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\mu_{12S}^{*}\mu_{21S}^{*}}=\arg\left[\mu_{12S}^{*}\mu_{21S}^{*}y_{\psi}^{2}\right]\,,$
( C.15)
$\displaystyle\delta_{\mu_{12S}^{*}\mu_{21S}}=\arg\left[\mu_{12S}^{*}\mu_{21S}(v_{1}v_{2}^{*})^{2}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\lambda^{\prime}_{1S}}=\arg\left[(\lambda^{\prime}_{1S})^{*}y_{\psi}^{2}\right]\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\delta_{\lambda^{\prime}_{2S}}=\arg\left[(\lambda^{\prime}_{2S})^{*}y_{\psi}^{2}\right]\,.$
( C.16)

Again, if $m_{12}^{2}=0$, the minimization condition implies ${\rm
Im}\lambda_{5}=0$ which means $\delta_{\lambda_{5}}=0$. However, we still have
lots of CPV phases in addition to $\theta_{\rm phys}=\delta_{\mu_{1S}}$.

##  Appendix D Details of calculations

In this appendix, we show the details of the calculations of the BZ diagrams.

###  D.1 General dark sector

Using couplings in Eq. (3.11), the left diagrams in Fig. 2 have contributions
as

$\displaystyle\int\\!\frac{d^{4}k}{(2\pi)^{4}}i\Gamma_{f,W}^{\mu\nu}\bar{u}(p_{2})\left\\{i\left(g_{eG}^{V}\gamma^{\alpha}+g_{eG}^{A}\gamma^{\alpha}\gamma^{5}\right)\right\\}\frac{i(p/_{2}-k/+m_{e})}{(p_{2}-k)^{2}-m_{e}^{2}}\left\\{i\left(g_{e\Phi}^{S}+ig_{e\Phi}^{P}\gamma^{5}\right)\right\\}u(p_{1})$
$\displaystyle\hskip
50.00008pt\times\frac{-ig_{\nu\alpha}}{k^{2}-m_{G}^{2}}\frac{i}{(-k-q)^{2}-m_{\eta_{a}}^{2}}\left(iF_{\rm
DS}^{\phi_{j}\eta_{a}}((-k-q)^{2},M^{2})\right)\frac{i}{(-k-q)^{2}-m_{\phi_{j}}^{2}}\,,$
( D.1)

where effective vertices of photon-gauge boson ($\gamma/Z$)-scalar
$\Gamma_{f,W}^{\mu\nu}$ have information about the inner loop with fermion and
$W$ boson, respectively, $G$ represents the photon or $Z$ boson, and $\Phi$ is
$\phi_{j}$ or $\eta_{a}$. Note that we should sum up all contributions for the
final result. For diagrams in Fig. 2, $\Gamma_{f,W}^{\mu\nu}$ can be
parameterized by

$\displaystyle i\Gamma_{f,W}^{\mu\nu}$
$\displaystyle=i\int_{0}^{1}\\!dx\frac{1}{x(1-x)}\int_{0}^{1-x}\\!dy\frac{C_{E}^{f,W}\left(k^{\mu}q^{\nu}-k\cdot
q\,g^{\mu\nu}\right)+C_{O}^{f,W}\epsilon^{\mu\nu\rho\sigma}q_{\rho}k_{\sigma}}{k^{2}-\Delta_{f,W}}\,.$
( D.2)

Here, $C_{E}^{f,W}$ and $C_{O}^{f,W}$ are obtained from the inner loop
calculation as [46, 36]

$\displaystyle C_{E}^{f}$
$\displaystyle=+\frac{N_{C}^{f}}{2\pi^{2}}g_{f\gamma}^{V}g_{fG}^{V}g_{f\Phi}^{S}m_{f}(1-4xy)\,,$
( D.3) $\displaystyle C_{O}^{f}$
$\displaystyle=-\frac{N_{C}^{f}}{2\pi^{2}}g_{f\gamma}^{V}g_{fG}^{V}g_{f\Phi}^{P}m_{f}\,,$
( D.4) $\displaystyle C_{E}^{W}$ $\displaystyle=+\frac{\rm
e}{4\pi^{2}}g_{WW\Phi}g_{WWG}\left\\{\left(4-\frac{m_{G}^{2}}{m_{W}^{2}}\right)-\left[6-\frac{m_{G}^{2}}{m_{W}^{2}}+\left(1-\frac{m_{G}^{2}}{2m_{W}^{2}}\right)\frac{m_{\Phi}^{2}}{m_{W}^{2}}\right]xy\right\\}\,,$
( D.5) $\displaystyle C_{O}^{W}$ $\displaystyle=0\,,$ ( D.6)

and $\Delta^{f,W}$ is defined in Eq. (3.14). Similarly, the right diagrams in
Fig. 2 are

$\displaystyle\int\\!\frac{d^{4}k}{(2\pi)^{4}}i\Gamma_{f,W}^{\mu\nu}\bar{u}(p_{2})\left\\{i\left(g_{e\Phi}^{S}+ig_{e\Phi}^{P}\gamma^{5}\right)\right\\}\frac{i(p/_{1}+k/+m_{e})}{(p_{1}+k)^{2}-m_{e}^{2}}\left\\{i\left(g_{eG}^{V}\gamma^{\alpha}+g_{eG}^{A}\gamma^{\alpha}\gamma^{5}\right)\right\\}u(p_{1})$
$\displaystyle\hskip
50.00008pt\times\frac{-ig_{\nu\alpha}}{k^{2}-m_{G}^{2}}\frac{i}{(-k-q)^{2}-m_{\eta_{a}}^{2}}\left(iF_{\rm
DS}^{\phi_{j}\eta_{a}}((-k-q)^{2},M^{2})\right)\frac{i}{(-k-q)^{2}-m_{\phi_{j}}^{2}}\,.$
( D.7)

Since $\Gamma_{f,W}^{\mu\nu}$ has one loop momentum $k$ in its numerator (see
Eq. (D.2)), the contributions to the electron EDM can be parameterized as

Eq. (D.1):
$\displaystyle\int_{x}\int_{y}\int_{k}\frac{A_{\alpha\beta}^{L,\Phi(f,W)}k^{\alpha}k^{\beta}+B_{\alpha}^{L,\Phi(f,W)}k^{\alpha}}{(k^{2}-\Delta_{f,W})\cdots(k^{2}-2k\cdot
p_{2})}F_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})\sigma^{\mu\nu}\gamma^{5}q_{\nu}\,,$ (
D.8) Eq. (D.7):
$\displaystyle\int_{x}\int_{y}\int_{k}\frac{A_{\alpha\beta}^{R,\Phi(f,W)}k^{\alpha}k^{\beta}+B_{\alpha}^{R,\Phi(f,W)}k^{\alpha}}{(k^{2}-\Delta_{f,W})\cdots(k^{2}+2k\cdot
p_{1})}F_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})\sigma^{\mu\nu}\gamma^{5}q_{\nu}\,,$ (
D.9)

where we omit some propagators for simplicity, and
$\epsilon^{\mu\nu\alpha\beta}\gamma_{\alpha}\gamma_{\beta}=-i[\gamma^{\mu},\gamma^{\nu}]\gamma^{5}=-2\sigma^{\mu\nu}\gamma^{5}$
is used. $A_{\alpha\beta}^{L,\Phi(f,W)}$, $A_{\alpha\beta}^{R,\Phi(f,W)}$,
$B_{\alpha}^{L,\Phi(f,W)}$ and $B_{\alpha}^{R,\Phi(f,W)}$ are independent of
the loop momentum $k$, but depend on $p_{1,2}$ and $q$ as well as $x$ and $y$
through $C_{E,O}^{f,W}$. These expressions are generally obtained as

$\displaystyle A_{\alpha\beta}^{L,\Phi(f,W)}$
$\displaystyle=-\frac{2}{\left(4m_{e}^{2}-q^{2}\right)q^{2}}\left[g_{eG}^{V}\left(C_{E}^{f,W}g_{e\Phi}^{P}+C_{O}^{f,W}g_{e\Phi}^{S}\right)+ig_{eG}^{A}\left(C_{E}^{f,W}g_{e\Phi}^{S}-C_{O}^{f,W}g_{e\Phi}^{P}\right)\right]$
$\displaystyle\hskip
90.00014pt\times\Bigl{\\{}2m_{e}^{2}(p_{1\alpha}p_{1\beta}+p_{2\alpha}p_{2\beta})-(2m_{e}^{2}-q^{2})(p_{1\alpha}p_{2\beta}+p_{1\beta}p_{2\alpha})\Bigr{\\}}$
$\displaystyle\hskip
11.99998pt+C_{O}^{f,W}\left(g_{eG}^{V}g_{e\Phi}^{S}-ig_{eG}^{A}g_{e\Phi}^{P}\right)g_{\alpha\beta}\,,$
( D.10) $\displaystyle A_{\alpha\beta}^{R,\Phi(f,W)}$
$\displaystyle=-\frac{2}{\left(4m_{e}^{2}-q^{2}\right)q^{2}}\left[g_{eG}^{V}\left(C_{E}^{f,W}g_{e\Phi}^{P}+C_{O}^{f,W}g_{e\Phi}^{S}\right)-ig_{eG}^{A}\left(C_{E}^{f,W}g_{e\Phi}^{S}-C_{O}^{f,W}g_{e\Phi}^{P}\right)\right]$
$\displaystyle\hskip
90.00014pt\times\Bigl{\\{}2m_{e}^{2}(p_{1\alpha}p_{1\beta}+p_{2\alpha}p_{2\beta})-(2m_{e}^{2}-q^{2})(p_{1\alpha}p_{2\beta}+p_{1\beta}p_{2\alpha})\Bigr{\\}}$
$\displaystyle\hskip
11.99998pt+C_{O}^{f,W}\left(g_{eG}^{V}g_{e\Phi}^{S}+ig_{eG}^{A}g_{e\Phi}^{P}\right)g_{\alpha\beta}\,,$
( D.11) $\displaystyle B_{\alpha}^{L,\Phi(f,W)}$
$\displaystyle=-\frac{4C_{E}^{f,W}g_{eG}^{V}g_{e\Phi}^{P}m_{e}^{2}}{4m_{e}^{2}-q^{2}}p_{1\alpha}+\left[\frac{2C_{E}^{f,W}(2m_{e}^{2}-q^{2})g_{eG}^{V}g_{e\Phi}^{P}}{4m_{e}^{2}-q^{2}}+2iC_{E}^{f,W}g_{eG}^{A}g_{e\Phi}^{S}\right]p_{2\alpha}\,,$
( D.12) $\displaystyle B_{\alpha}^{R,\Phi(f,W)}$
$\displaystyle=-\left[\frac{2C_{E}^{f,W}(2m_{e}^{2}-q^{2})g_{eG}^{V}g_{e\Phi}^{P}}{4m_{e}^{2}-q^{2}}-2iC_{E}^{f,W}g_{eG}^{A}g_{e\Phi}^{S}\right]p_{1\alpha}+\frac{4C_{E}^{f,W}g_{eG}^{V}g_{e\Phi}^{P}m_{e}^{2}}{4m_{e}^{2}-q^{2}}p_{2\alpha}\,.$
( D.13)

Note that $A_{\alpha\beta}^{L,\Phi(f,W)}$ and $A_{\alpha\beta}^{R,\Phi(f,W)}$
are symmetric under $\alpha\leftrightarrow\beta$.

###  D.2 Heavy dark sector

If we expand $F_{\rm DS}^{\phi_{j}\eta_{a}}$ by $k^{2}/M^{2}$, as in Eq.
(3.16), we can easily perform the loop momentum integration. The HDS
approximation result is then given by

$\displaystyle d_{e}=\,$
$\displaystyle\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\left[\sum_{f}\left(d_{e}^{f}\right)_{\phi_{j}\eta_{a}}+\left(d_{e}^{W}\right)_{\phi_{j}\eta_{a}}\right]\,,$
( D.14)
$\displaystyle\left(d_{e}^{\mathcal{P}}\right)_{\phi_{j}\eta_{a}}\simeq\,$
$\displaystyle\sum_{\Phi=\phi_{j},\eta_{a}}\int_{x}\int_{y}\int_{\bm{z}}\Biggl{\\{}\Biggr{.}a_{0}^{\phi_{j}\eta_{a}}(M^{2})\left[\frac{A_{\alpha\beta}^{L,\Phi(\mathcal{P})}+A_{\alpha\beta}^{R,\Phi(\mathcal{P})}}{32\pi^{2}\Delta_{5}(\mathcal{P})^{2}}g^{\alpha\beta}-\frac{f_{5}^{L,\Phi(\mathcal{P})}+f_{5}^{R,\Phi(\mathcal{P})}}{8\pi^{2}\Delta_{5}(\mathcal{P})^{3}}\right]$
$\displaystyle\hskip
60.00009pt+\frac{a_{1}^{\phi_{j}\eta_{a}}(M^{2})}{M^{2}}\Biggl{[}\Biggr{.}-\frac{3(A_{\alpha\beta}^{L,\Phi(\mathcal{P})}+A_{\alpha\beta}^{R,\Phi(\mathcal{P})})}{32\pi^{2}\Delta_{5}(\mathcal{P})}g^{\alpha\beta}$
$\displaystyle\hskip
75.0001pt+\frac{\bar{B}^{L,\Phi(\mathcal{P})}\cdot\lambda^{L}+\bar{B}^{R,\Phi(\mathcal{P})}\cdot\lambda^{R}}{32\pi^{2}\Delta_{5}(\mathcal{P})^{2}}+\frac{f_{5}^{L,\Phi(\mathcal{P})}+f_{5}^{R,\Phi(\mathcal{P})}}{8\pi^{2}\Delta_{5}(\mathcal{P})^{2}}$
$\displaystyle\hskip
75.0001pt+\left(\frac{A_{\alpha\beta}^{L,\Phi(\mathcal{P})}+A_{\alpha\beta}^{R,\Phi(\mathcal{P})}}{32\pi^{2}\Delta_{5}(\mathcal{P})^{2}}g^{\alpha\beta}-\frac{f_{5}^{L,\Phi(\mathcal{P})}+f_{5}^{R,\Phi(\mathcal{P})}}{8\pi^{2}\Delta_{5}(\mathcal{P})^{3}}\right)C_{p}^{2}m_{e}^{2}\Biggl{.}\Biggr{]}\Biggl{.}\Biggr{\\}}\,,$
( D.15)

where $\mathcal{P}=f$ or $W$, $\mathcal{O}(k^{4}/M^{4})$ terms are omitted,
and we define

$\displaystyle\int_{\bm{z}}$
$\displaystyle\equiv\int_{0}^{1}\\!dz_{1}dz_{2}dz_{3}dz_{4}dz_{5}\delta(1-z_{1}-z_{2}-z_{3}-z_{4}-z_{5})\,,$
( D.16) $\displaystyle\Delta_{5}(f,W)$
$\displaystyle\equiv\frac{z_{1}m_{f,W}^{2}}{x(1-x)}+z_{2}m_{\phi_{j}}^{2}+z_{3}m_{\eta_{a}}^{2}+z_{4}m_{G}^{2}+z_{5}^{2}m_{e}^{2}\,,$
( D.17) $\displaystyle\bar{B}_{\alpha}^{L,\Phi(f,W)}$
$\displaystyle\equiv-2A_{\alpha\beta}^{L,\Phi(f,W)}(C_{q}q^{\beta}-C_{p}p_{2}^{\beta})+B_{\alpha}^{L,\Phi(f,W)}\,,$
( D.18) $\displaystyle\bar{B}_{\alpha}^{R,\Phi(f,W)}$
$\displaystyle\equiv-2A_{\alpha\beta}^{R,\Phi(f,W)}(C_{q}q^{\beta}+C_{p}p_{1}^{\beta})+B_{\alpha}^{R,\Phi(f,W)}\,,$
( D.19) $\displaystyle f_{5}^{L,\Phi(f,W)}$ $\displaystyle\equiv
A_{\alpha\beta}^{L,\Phi(f,W)}(C_{q}q^{\alpha}-C_{p}p_{2}^{\alpha})(C_{q}q^{\beta}-C_{p}p_{2}^{\beta})-B_{\alpha}^{L,\Phi(f,W)}(C_{q}q^{\alpha}-C_{p}p_{2}^{\alpha})\,,$
( D.20) $\displaystyle f_{5}^{R,\Phi(f,W)}$ $\displaystyle\equiv
A_{\alpha\beta}^{R,\Phi(f,W)}(C_{q}q^{\alpha}+C_{p}p_{1}^{\alpha})(C_{q}q^{\beta}+C_{p}p_{1}^{\beta})-B_{\alpha}^{R,\Phi(f,W)}(C_{q}q^{\alpha}+C_{p}p_{1}^{\alpha})\,,$
( D.21) $\displaystyle\lambda^{L}$ $\displaystyle\equiv
2(1-C_{q})q+2C_{p}p_{2}\,,$ ( D.22) $\displaystyle\lambda^{R}$
$\displaystyle\equiv 2(1-C_{q})q-2C_{p}p_{1}\,,$ ( D.23) $\displaystyle C_{q}$
$\displaystyle=\frac{y}{1-x}z_{1}+z_{2}+z_{3}\,,\quad
C_{p}=(1-z_{1}-z_{2}-z_{3}-z_{4})\,.$ ( D.24)

Here, $C_{q}$ and $C_{p}$ correspond to the momentum shift originated from the
$k\cdot p_{i}$ and $k\cdot q$ terms in the denominator,

$\displaystyle k\to\tilde{k}=\begin{cases}k+C_{q}q-C_{p}p_{2}&\text{for
Eq.\leavevmode\nobreak\ \eqref{eq:outerloopcalc1}}\\\
k+C_{q}q+C_{p}p_{1}&\text{for Eq.\leavevmode\nobreak\
\eqref{eq:outerloopcalc2}}\end{cases}\,.$ ( D.25)

Note that due to $q^{2}=0$ which leads to $p_{1,2}\cdot q=0$, we have
$\lambda^{L}\cdot\lambda^{L}=\lambda^{R}\cdot\lambda^{R}=4C_{p}^{2}m_{e}^{2}$.
For the diagram with $\Gamma^{\mu\nu}$ for the inner loop, each term can be
calculated as

$\displaystyle(A_{\alpha\beta}^{L,\Phi(f,W)}+A_{\alpha\beta}^{R,\Phi(f,W)})g^{\alpha\beta}$
$\displaystyle=-4g^{V}_{e}\left(C_{E}^{f,W}g^{P}_{e}-C_{O}^{f,W}g^{S}_{e}\right)\,,$
( D.26)
$\displaystyle\bar{B}^{L,\Phi(f,W)}\cdot\lambda^{L}+\bar{B}^{R,\Phi(f,W)}\cdot\lambda^{R}$
$\displaystyle=-8C_{p}^{2}C_{E}^{f,W}g^{V}_{e}g^{P}_{e}m_{e}^{2}\,,$ ( D.27)
$\displaystyle f_{5}^{L,\Phi(f,W)}+f_{5}^{R,\Phi(f,W)}$
$\displaystyle=-2C_{p}^{2}C_{E}^{f,W}g^{V}_{e}g^{P}_{e}m_{e}^{2}\,.$ ( D.28)

It should be emphasized that the effects of $k\cdot p_{i}$ and $k\cdot q$ in
the denominator of Eqs. (D.1) and (D.7) appear in the final result as sub-
dominant terms: once we set $C_{q,p}=0$ which means no momentum shift for $k$,
Eq. (D.15) becomes the leading results in Eq. (3.17). This can be understood
by the fact that Eqs. (D.27) and (D.28) and the last term in Eq. (D.15) become
zero. This feature remains even for the general result without the HDS
approximation shown in Eq. (3.13).

###  D.3 Decomposed expressions for the electron EDM

In Eq. (4.26), we define two parts of contributions to the electron EDM:
$F_{\rm DS}^{\phi_{j}\eta_{a}}$ and the renormalization part. Here, the
explicit form of each contribution is shown.

Let us start with the case of a general dark sector. In our setup, we have
effective mixing between CP-even ($\phi_{j}$) and CP-odd ($\eta_{a}$) scalars
via a dark sector CP violation, denoted as $F_{\rm DS}^{\phi_{j}\eta_{a}}$. If
$F_{\rm DS}^{\phi\eta_{a}}$ has a UV divergent part, as in the two benchmark
models, the divergence needs to be treated appropriately. In the on-shell
subtraction scheme, we can remove it and obtain a renormalized one:

$\displaystyle\hat{\Sigma}_{\phi_{j}\eta_{a}}=F_{\rm
DS}^{\phi_{j}\eta_{a}}+\delta F_{\rm DS}^{\phi_{j}\eta_{a}}\,,$ ( D.29)

where $\delta F_{\rm DS}^{\phi_{j}\eta_{a}}$ corresponds to the term induced
by renormalization. Then, the UV divergence is canceled, and hence, we find a
finite prediction,

$\displaystyle\hat{\Sigma}_{\phi_{j}\eta_{a}}=\widetilde{F}_{\rm
DS}^{\phi_{j}\eta_{a}}+\widetilde{\delta F}_{\rm DS}^{\phi_{j}\eta_{a}}\,.$ (
D.30)

Here, $\widetilde{F}_{\rm DS}^{\phi_{j}\eta_{a}}$ and $\widetilde{\delta
F}_{\rm DS}^{\phi_{j}\eta_{a}}$ do not have any divergent part. Note that they
also do not have any constant term associated with the momentum passing
through this effective mixing. Using this notation, we can now find the
explicit formulas for $d_{e}^{\rm Org.,Ren.}$ as

$\displaystyle d_{e}^{\rm Org.}$
$\displaystyle=\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\int_{x}\int_{y}\int_{k}\frac{\widetilde{F}_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})}{(k^{2}-\Delta)(k^{2}-m_{G}^{2})(k^{2}+2k\cdot
q-m_{\phi_{j}}^{2})(k^{2}+2k\cdot q-m_{\eta_{a}}^{2})}$ $\displaystyle\hskip
10.00002pt\times\sum_{\Phi=\phi_{j},\eta_{a}}\left(\frac{A_{\alpha\beta}^{L,\Phi}k^{\alpha}k^{\beta}+B_{\alpha}^{L,\Phi}k^{\alpha}}{k^{2}-2k\cdot
p_{2}}+\frac{A_{\alpha\beta}^{R,\Phi}k^{\alpha}k^{\beta}+B_{\alpha}^{R,\Phi}k^{\alpha}}{k^{2}+2k\cdot
p_{1}}\right)\,,$ ( D.31) $\displaystyle d_{e}^{\rm Ren.}$
$\displaystyle=\sum_{j=1}^{N_{E}}\sum_{a=2}^{N_{O}}\int_{x}\int_{y}\int_{k}\frac{\widetilde{\delta
F}_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})}{(k^{2}-\Delta)(k^{2}-m_{G}^{2})(k^{2}+2k\cdot
q-m_{\phi_{j}}^{2})(k^{2}+2k\cdot q-m_{\eta_{a}}^{2})}$ $\displaystyle\hskip
10.00002pt\times\sum_{\Phi=\phi_{j},\eta_{a}}\left(\frac{A_{\alpha\beta}^{L,\Phi}k^{\alpha}k^{\beta}+B_{\alpha}^{L,\Phi}k^{\alpha}}{k^{2}-2k\cdot
p_{2}}+\frac{A_{\alpha\beta}^{R,\Phi}k^{\alpha}k^{\beta}+B_{\alpha}^{R,\Phi}k^{\alpha}}{k^{2}+2k\cdot
p_{1}}\right)\,,$ ( D.32)

where we omit the superscript $f$ and $W$ from the general result in Eq.
(3.13) for simplicity. Note that for $d_{e}^{\rm Ren.}$, a further calculation
can be done as in the case of the HDS limit in section 3.4.1, because the on-
shell subtraction scheme gives additional terms which are proportional to the
squared momentum and constant (see Eqs. (4.16), (4.17) and (4.44)). For the
HDS approximation, we expand $\widetilde{F}_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})$ and $\widetilde{\delta F}_{\rm
DS}^{\phi_{j}\eta_{a}}((k+q)^{2},M^{2})$ with respect to $k^{2}/M^{2}$.

For the scalar dark sector model, we have

$\displaystyle\widetilde{F}_{\rm DS}^{hA}(k^{2})$ $\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})+{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})$ $\displaystyle\hskip
135.0002pt-2{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{2}})+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{k^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\biggl{.}\biggr{)}$
$\displaystyle\hskip 60.00009pt+\lambda_{hS}\left({\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})-{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})\right)\Biggl{.}\Biggr{]}\,,$ (
D.33) $\displaystyle\widetilde{\delta F}_{\rm DS}^{hA}(k^{2})$
$\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\biggl{[}\biggr{.}\frac{k^{2}-m_{A}^{2}}{m_{A}^{2}-m_{h}^{2}}\Bigl{(}\Bigr{.}\lambda_{12S}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\overline{\mathcal{F}}_{1}(m_{h},m_{\varphi_{1}},m_{\varphi_{2}})$
$\displaystyle\hskip
135.0002pt+\lambda_{hS}\mathcal{F}_{2}(m_{h},m_{\varphi_{1}},m_{\varphi_{2}})\Bigl{.}\Bigr{)}$
$\displaystyle\hskip
72.00008pt+\frac{m_{h}^{2}-k^{2}}{m_{A}^{2}-m_{h}^{2}}\Bigl{(}\lambda_{12S}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\overline{\mathcal{F}}_{1}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})$
$\displaystyle\hskip
135.0002pt+\lambda_{hS}\mathcal{F}_{2}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{)}\biggl{.}\biggr{]}\,,$
( D.34) $\displaystyle\widetilde{F}_{\rm DS}^{HA}(k^{2})$
$\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})+{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})$ $\displaystyle\hskip
135.0002pt-2{\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{2}})+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{k^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\biggl{.}\biggr{)}$
$\displaystyle\hskip 60.00009pt+\lambda_{HS}\left({\rm
DiscB}(k^{2},m_{\varphi_{1}},m_{\varphi_{1}})-{\rm
DiscB}(k^{2},m_{\varphi_{2}},m_{\varphi_{2}})\right)\Biggl{.}\Biggr{]}\,,$ (
D.35) $\displaystyle\widetilde{\delta F}_{\rm DS}^{HA}(k^{2})$
$\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\biggl{[}\biggr{.}\frac{k^{2}-m_{A}^{2}}{m_{A}^{2}-m_{H}^{2}}\Bigl{(}\Bigr{.}\lambda_{12S}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\overline{\mathcal{F}}_{1}(m_{H},m_{\varphi_{1}},m_{\varphi_{2}})$
$\displaystyle\hskip
135.0002pt+\lambda_{HS}\mathcal{F}_{2}(m_{H},m_{\varphi_{1}},m_{\varphi_{2}})\Bigl{.}\Bigr{)}$
$\displaystyle\hskip
72.00008pt+\frac{m_{H}^{2}-k^{2}}{m_{A}^{2}-m_{H}^{2}}\Bigl{(}\lambda_{12S}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\overline{\mathcal{F}}_{1}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})$
$\displaystyle\hskip
135.0002pt+\lambda_{hS}\mathcal{F}_{2}(m_{A},m_{\varphi_{1}},m_{\varphi_{2}})\Bigr{)}\biggl{.}\biggr{]}\,,$
( D.36)

where

$\displaystyle\overline{\mathcal{F}}_{1}(m,m_{\varphi_{1}},m_{\varphi_{2}})\equiv\mathcal{F}_{1}(m,m_{\varphi_{1}},m_{\varphi_{2}})+\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{1}}^{2}}{m^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\,,$ ( D.38)

and $\mathcal{F}_{1}(m,m_{\varphi_{1}},m_{\varphi_{2}})$ and
$\mathcal{F}_{2}(m,m_{\varphi_{1}},m_{\varphi_{2}})$ are given in Eqs. (4.23)
and (4.24), respectively. For the HDS approximation, we can obtain the
expressions by expanding each of them in the heavy $m_{\varphi_{1,2}}$ limit
as

$\displaystyle\widetilde{F}_{\rm DS}^{hA,{\rm HDS}}(k^{2})$
$\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}c_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}-2+\frac{m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2}}{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}$ $\displaystyle\hskip
80.00012pt+k^{2}\left\\{\frac{(m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2})(m_{\varphi_{1}}^{4}-8m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}+m_{\varphi_{2}}^{4})}{6m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{2}}\right.$
$\displaystyle\hskip
100.00015pt\left.+\frac{2m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}}{(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{3}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\right\\}\biggl{.}\biggr{)}-\lambda_{hS}\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{6m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}}k^{2}\Biggl{.}\Biggr{]}\,,$
( D.39) $\displaystyle\widetilde{F}_{\rm DS}^{HA,{\rm HDS}}(k^{2})$
$\displaystyle=-\frac{v_{\rm
SM}^{2}}{64\pi^{2}}\lambda_{12S}^{-}s_{2\theta_{s}}\Biggl{[}\Biggr{.}\lambda_{12S}^{+}s_{\alpha+\beta}c_{2\theta_{s}}\biggl{(}\biggr{.}-2+\frac{m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2}}{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}$ $\displaystyle\hskip
80.00012pt+k^{2}\left\\{\frac{(m_{\varphi_{1}}^{2}+m_{\varphi_{2}}^{2})(m_{\varphi_{1}}^{4}-8m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}+m_{\varphi_{2}}^{4})}{6m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{2}}\right.$
$\displaystyle\hskip
100.00015pt\left.+\frac{2m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}}{(m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2})^{3}}\mathop{\rm
ln}\frac{m_{\varphi_{1}}^{2}}{m_{\varphi_{2}}^{2}}\right\\}\biggl{.}\biggr{)}-\lambda_{HS}\frac{m_{\varphi_{1}}^{2}-m_{\varphi_{2}}^{2}}{6m_{\varphi_{1}}^{2}m_{\varphi_{2}}^{2}}k^{2}\Biggl{.}\Biggr{]}\,.$
( D.40)

Here, we omit $\mathcal{O}(k^{4}/m_{\varphi_{i}}^{4})$ terms, and
$\widetilde{\delta F}_{\rm DS}^{hA,{\rm HDS}}(k^{2})$ and $\widetilde{\delta
F}_{\rm DS}^{HA,{\rm HDS}}(k^{2})$ are the same as $\widetilde{\delta F}_{\rm
DS}^{hA}(k^{2})$ and $\widetilde{\delta F}_{\rm DS}^{HA}(k^{2})$.

For the fermion dark sector model, we find

$\displaystyle\widetilde{F}_{\rm DS}^{\phi_{j}\eta_{a}}(k^{2})$
$\displaystyle=\frac{y_{\psi}^{R}y_{\psi}^{I}}{8\pi^{2}}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}{\rm
DiscB}(k^{2},m_{\psi},m_{\psi})\,,$ ( D.41) $\displaystyle\widetilde{\delta
F}_{\rm DS}^{\phi_{j}\eta_{a}}(k^{2})$
$\displaystyle=\frac{y_{\psi}^{R}y_{\psi}^{I}}{8\pi^{2}}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}\biggl{[}\biggr{.}\frac{m_{\eta_{a}}^{2}-k^{2}}{m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2}}{\rm
Re}\bigl{\\{}{\rm DiscB}(m_{\phi_{j}}^{2},m_{\psi},m_{\psi})\bigr{\\}}$
$\displaystyle\hskip
115.00017pt+\frac{k^{2}-m_{\phi_{j}}^{2}}{m_{\phi_{j}}^{2}-m_{\eta_{a}}^{2}}{\rm
Re}\bigl{\\{}{\rm
DiscB}(m_{\eta_{a}}^{2},m_{\psi},m_{\psi})\bigr{\\}}\biggl{.}\biggr{]}\,.$ (
D.42)

For the HDS approximation, they are

$\displaystyle\widetilde{F}_{\rm DS}^{\phi_{j}\eta_{a},{\rm HDS}}(k^{2})$
$\displaystyle=\frac{y_{\psi}^{R}y_{\psi}^{I}}{8\pi^{2}}m_{\psi}^{2}(R_{E})_{3j}(R_{O})_{3a}\left(-2+\frac{k^{2}}{6m_{\psi}^{2}}\right)\,,$
( D.43)

and $\widetilde{\delta F}_{\rm DS}^{\phi_{j}\eta_{a},{\rm HDS}}(k^{2})$
unchanged from $\widetilde{\delta F}_{\rm DS}^{\phi_{j}\eta_{a}}(k^{2})$.

## References

  * [1] ATLAS collaboration, _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1 [1207.7214].
  * [2] CMS collaboration, _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30 [1207.7235].
  * [3] T. D. Lee, _A Theory of Spontaneous T Violation_ , _Phys. Rev. D_ 8 (1973) 1226.
  * [4] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher and J. P. Silva, _Theory and phenomenology of two-Higgs-doublet models_ , _Phys. Rept._ 516 (2012) 1 [1106.0034].
  * [5] K. Inoue, A. Kakuto, H. Komatsu and S. Takeshita, _Low Energy Parameters and Particle Masses in a Supersymmetric Grand Unified Model_ , _Progress of Theoretical Physics_ 67 (1982) 1889 [https://academic.oup.com/ptp/article-pdf/67/6/1889/5333146/67-6-1889.pdf].
  * [6] P. Fayet, _Supergauge invariant extension of the higgs mechanism and a model for the electron and its neutrino_ , _Nuclear Physics B_ 90 (1975) 104.
  * [7] J. F. Gunion and H. E. Haber, _Higgs Bosons in Supersymmetric Models. 1._ , _Nucl. Phys. B_ 272 (1986) 1.
  * [8] R. A. Flores and M. Sher, _Higgs masses in the standard, multi-higgs and supersymmetric models_ , _Annals of Physics_ 148 (1983) 95.
  * [9] N. G. Deshpande and E. Ma, _Pattern of Symmetry Breaking with Two Higgs Doublets_ , _Phys. Rev. D_ 18 (1978) 2574.
  * [10] L. Lopez Honorez, E. Nezri, J. F. Oliver and M. H. G. Tytgat, _The Inert Doublet Model: An Archetype for Dark Matter_ , _JCAP_ 02 (2007) 028 [hep-ph/0612275].
  * [11] E. M. Dolle and S. Su, _The Inert Dark Matter_ , _Phys. Rev. D_ 80 (2009) 055012 [0906.1609].
  * [12] L. Lopez Honorez and C. E. Yaguna, _The inert doublet model of dark matter revisited_ , _JHEP_ 09 (2010) 046 [1003.3125].
  * [13] N. Blinov, J. Kozaczuk, D. E. Morrissey and A. de la Puente, _Compressing the Inert Doublet Model_ , _Phys. Rev. D_ 93 (2016) 035020 [1510.08069].
  * [14] M. Krawczyk, N. Darvishi and D. Sokolowska, _The Inert Doublet Model and its extensions_ , _Acta Phys. Polon. B_ 47 (2016) 183 [1512.06437].
  * [15] M. Carena, M. Quirós and Y. Zhang, _Electroweak Baryogenesis from Dark-Sector CP Violation_ , _Phys. Rev. Lett._ 122 (2019) 201802 [1811.09719].
  * [16] M. Carena, M. Quirós and Y. Zhang, _Dark CP violation and gauged lepton or baryon number for electroweak baryogenesis_ , _Phys. Rev. D_ 101 (2020) 055014 [1908.04818].
  * [17] ACME collaboration, _Improved limit on the electric dipole moment of the electron_ , _Nature_ 562 (2018) 355.
  * [18] T. S. Roussy et al., _An improved bound on the electron’s electric dipole moment_ , _Science_ 381 (2023) adg4084 [2212.11841].
  * [19] L. Caldwell et al., _Systematic and statistical uncertainty evaluation of the HfF+ electron electric dipole moment experiment_ , _Phys. Rev. A_ 108 (2023) 012804 [2212.11837].
  * [20] S. Okawa, M. Pospelov and A. Ritz, _Electric Dipole Moments From Dark Sectors_ , _Phys. Rev. D_ 100 (2019) 075017 [1905.05219].
  * [21] T.-K. Chen, C.-W. Chiang and I. Low, _Simple model of dark matter and CP violation_ , _Phys. Rev. D_ 105 (2022) 075025 [2202.02954].
  * [22] D. Azevedo, P. M. Ferreira, M. M. Muhlleitner, S. Patel, R. Santos and J. Wittbrodt, _CP in the dark_ , _JHEP_ 11 (2018) 091 [1807.10322].
  * [23] D. Egana-Ugrinovic and S. Thomas, _Higgs Boson Contributions to the Electron Electric Dipole Moment_ , 1810.08631.
  * [24] E. J. Chun, J. Kim and T. Mondal, _Electron EDM and Muon anomalous magnetic moment in Two-Higgs-Doublet Models_ , _JHEP_ 12 (2019) 068 [1906.00612].
  * [25] K. Cheung, A. Jueid, Y.-N. Mao and S. Moretti, _Two-Higgs-doublet model with soft $CP$ violation confronting electric dipole moments and colliders_, _Phys. Rev. D_ 102 (2020) 075029 [2003.04178].
  * [26] S. Kanemura, M. Kubota and K. Yagyu, _Aligned CP-violating Higgs sector canceling the electric dipole moment_ , _JHEP_ 08 (2020) 026 [2004.03943].
  * [27] N. Chen, T. Li, Z. Teng and Y. Wu, _Collapsing domain walls in the two-Higgs-doublet model and deep insights from the EDM_ , _JHEP_ 10 (2020) 081 [2006.06913].
  * [28] I. Low, N. R. Shah and X.-P. Wang, _Higgs alignment and novel CP-violating observables in two-Higgs-doublet models_ , _Phys. Rev. D_ 105 (2022) 035009 [2012.00773].
  * [29] S. Kanemura, M. Kubota and K. Yagyu, _Testing aligned CP-violating Higgs sector at future lepton colliders_ , _JHEP_ 04 (2021) 144 [2101.03702].
  * [30] M. Frank, E. G. Fuakye and M. Toharia, _Restricting the parameter space of type-II two-Higgs-doublet models with CP violation_ , _Phys. Rev. D_ 106 (2022) 035010 [2112.14295].
  * [31] S. Kanemura, K. Katayama, T. Mondal and K. Yagyu, _Multi-photon signatures as a probe of CP-violation in extended Higgs sectors_ , 2308.01772.
  * [32] S. L. Glashow and S. Weinberg, _Natural Conservation Laws for Neutral Currents_ , _Phys. Rev. D_ 15 (1977) 1958.
  * [33] A. Djouadi, _The Anatomy of electro-weak symmetry breaking. II. The Higgs bosons in the minimal supersymmetric model_ , _Phys. Rept._ 459 (2008) 1 [hep-ph/0503173].
  * [34] Muon (g-2) collaboration, _An Improved Limit on the Muon Electric Dipole Moment_ , _Phys. Rev. D_ 80 (2009) 052008 [0811.1207].
  * [35] C. Abel et al., _Measurement of the Permanent Electric Dipole Moment of the Neutron_ , _Phys. Rev. Lett._ 124 (2020) 081803 [2001.11966].
  * [36] Y. Nakai and M. Reece, _Electric Dipole Moments in Natural Supersymmetry_ , _JHEP_ 08 (2017) 031 [1612.08090].
  * [37] C. Cesarotti, Q. Lu, Y. Nakai, A. Parikh and M. Reece, _Interpreting the Electron EDM Constraint_ , _JHEP_ 05 (2019) 059 [1810.07736].
  * [38] M. Pospelov and A. Ritz, _Electric dipole moments as probes of new physics_ , _Annals Phys._ 318 (2005) 119 [hep-ph/0504231].
  * [39] J. Engel, M. J. Ramsey-Musolf and U. van Kolck, _Electric Dipole Moments of Nucleons, Nuclei, and Atoms: The Standard Model and Beyond_ , _Prog. Part. Nucl. Phys._ 71 (2013) 21 [1303.2371].
  * [40] T. Chupp, P. Fierlinger, M. Ramsey-Musolf and J. Singh, _Electric dipole moments of atoms, molecules, nuclei, and particles_ , _Rev. Mod. Phys._ 91 (2019) 015001 [1710.02504].
  * [41] B. Graner, Y. Chen, E. G. Lindahl and B. R. Heckel, _Reduced Limit on the Permanent Electric Dipole Moment of Hg199_ , _Phys. Rev. Lett._ 116 (2016) 161601 [1601.04339].
  * [42] J. R. Ellis, J. S. Lee and A. Pilaftsis, _Electric Dipole Moments in the MSSM Reloaded_ , _JHEP_ 10 (2008) 049 [0808.1819].
  * [43] M. Jung and A. Pich, _Electric Dipole Moments in Two-Higgs-Doublet Models_ , _JHEP_ 04 (2014) 076 [1308.6283].
  * [44] S. M. Barr and A. Zee, _Electric dipole moment of the electron and of the neutron_ , _Phys. Rev. Lett._ 65 (1990) 21.