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

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.


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 O(10 −29 − 10 −30 ) ecm.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.Moreover, various researchers have explored the implications of the new dark gauge symmetry coupled with the scalar, generating CP-violation at two-loop levels as well [22].This implies that EDM originating from the dark sector can be either highly suppressed at higher loops or, if present at the tree level, plays a vital role in understanding the CP violation.
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. 1  A new CP-violating phase, which can be 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.

The framework
Let us start with an overview of the features of our CP-conserving 2HDMs and CPviolating dark sectors under consideration in the present paper.The 2HDM is a wellstudied 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.

Two Higgs Doublet Model
The potential of two Higgs doublets H 1,2 can be expressed as
Type-Y (d R : −, e R : +) −Y (2) e L L H 2 e R + h.c.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 .
remains invariant.This phase is severely constrained by the electron EDM measurement, as previously studied in refs.[24][25][26][27][28][29][30][31][32].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[λ * 5 (m2 12 ) 2 ] = 0, and both λ 5 and m 2 12 are real 2 .The two SU (2) L doublets H 1,2 are parametrized as 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 2 1 + v 2 2 = v 2 SM = (246.22GeV) 2 .We define their ratio as Upon diagonalizing the mass matrices derived from the scalar potential in Eq. (2.1), we obtain the mass eigenstates for the neutral scalars: with mixing angle α 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.The Yukawa couplings for SM fermions in the 2HDM are represented by the following Lagrangian term: where 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 [33].For example, H 1 behaves as odd while H 2 as even under a discrete Z 2 symmetry. 3Without loss of generality, one can select a scenario in which 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 where m f is the SM fermion mass, and the coefficients ξ f h,H,A are summarized in Table 2.It is notable that the lepton couplings to the heavy scalars are significantly enhanced by a large tan β for Type-II and Type-X 2HDMs, leading to a larger effect on the electron EDM.
Due to the additional Higgs doublet, the W boson couplings to CP-even scalars are modified as 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 [34].

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 where 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: where the gray shaded blob represents loops of dark sector particles.Nonzero mass mixing between CP-even (h 0 , H 0 ) and CP-odd (A 0 ) scalars in the 2HDM is induced by dark sector CP violation.
We have comment on a renormalization group (RG) running of couplings in 2HDM sector.When we consider a UV complete model that incorporates our dark sector, the RG running induces a phase for λ 5 as well as m 2  12 , due to the coupling between 2HDM and dark sectors in Eq. (2.9).In addition, there is the possibility of generating terms such as . through the RG running, which have complex parameters of λ 6,7 .These terms and λ 5 are additional sources of CP violation, inducing a CP-violating nature in our 2HDM sector where CP-even and CP-odd scalars mix.As mentioned in section 2.1, a CP-violating 2HDM is severely constrained by the bound on d e , and the current focus is electron EDM prediction from dark sector CP violation.Therefore, we specifically consider a scenario in which the phases of λ 5,6,7 induced by the RG running are sufficiently small.This ensures that contributions from these CP phases to the electron EDM are sub-dominant, allowing our EDM prediction to be predominantly influenced by CP violation in the dark sector.The details of effect of this RG running will be performed in the future work.

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.

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,35,36], which put stringent constraints on models incorporating CP violation, such as the complex 2HDM and supersymmetry (see e.g.refs.[37,38]).For an in-depth understanding of EDMs and their measurement techniques, along with theoretical aspects, one can refer to comprehensive reviews such as refs.[39][40][41].
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 ThO | < 1.1 × 10 −29 e 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: where σ µν = i 2 [γ µ , γ ν ], 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]   In addition to EDMs of paramagnetic molecules, there are experiments to measure EDMs of diamagnetic atoms, such as the mercury EDM [42] and the neutron EDM [36].These experiments, in principle, can impose constraints on the 2HDM or supersymmetry, and the details can be found in refs.[43,44].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.

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 ēiγ 5 e N N , while the right corresponds to the Barr-Zee (BZ) diagram that results in the electron EDM operator d e ēσ µν γ 5 eF µν [45].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 3 e : 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.
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 ϕ j and η a , respectively.For example, the 2HDM can be identified with (N E , N O ) = (2, 2) and (ϕ 1 , ϕ 2 , η 1 , η 2 ) = (h 0 , H 0 , G 0 , A 0 ).Note that we assume η 1 always represents the NG field G 0 , and thus, the sum of the index a should start from 2 for η a scalars.ϕ j and η a are considered as mass eigenstates, and hence, masses for these scalars are described as Here, we assume m 2 . The neutral scalars ϕ j , η 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: where ξ f ϕ j ,ηa and ξ W ϕ j are related to tan β and mixing angles for ϕ j and η 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 ϕ j and η 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: 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 ϕ j ηa DS (k 2 , M 2 ).

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, The coefficient C f f ′ 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 . Consequently, we obtain Couplings defined in Eq. (3.11).Here, g denotes the gauge coupling of SU (2) L , θ W is the weak mixing angle, e = g sin θ W is the electromagnetic coupling, where ξ f ϕ j and ξ f ′ ηa are defined in Eq. (3.5).The C f f ′ term contributes to C S in Eq. (3.1).The coefficient of the electron-quark interaction C qe leads to where numerators arise from the matrix elements given in refs.[43,46], with m u /m d ≃ 0.55.

The Barr-Zee diagrams
To calculate the BZ diagram in Fig. 1, it is important to handle the momentum in F ϕ j ηa DS 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 ϕ j ηa DS .Therefore, we present the general result before the loop momentum integration.Before the calculation, we define the relevant couplings for SM fermions as where Φ and G µ denote the neutral scalars (ϕ j , η a ) and neutral gauge bosons (Z µ , A µ ), respectively.Comparing with Eq. (3.5), g S,P f Φ are related to ξ f ϕ j ,ηa , and g V,A f G are the usual gauge boson couplings in the SM, as summarized in Table 3.
The relevant sub-diagrams for the BZ diagrams are shown in Fig. 2. 4 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. [47].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. [48].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 ) where we have used p 2 1,2 = m 2 e for the outer electron lines, q = p 1 − p 2 represents the four-momentum of the external photon, and ∆ f,W is defined as In addition, we have defined a shorthand notation for each integral as where x and y are Feynman parameters arising from the inner loop.The terms 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 •p i and k •q in the expression will contribute sub-dominantly.Upon integrating the loop momentum k, these terms lead to p i • p j = m 2 e or p i • q = 0, as q • q = 0.
(2) Due to the propagators, Eq. (3.13) exhibits a pole around k 2 ∼ v 2 SM if the neutral scalar masses in the visible sector are of O(v SM ).Consequently, the dominant contribution to d e is expected to occur when k 2 ∼ v 2 SM .

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 ≫ v SM .By power expanding F where a ϕ j ηa 0,1 (M 2 ) are independent of k, and each order is expected to be suppressed by O(v 2 SM /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): where we omit terms proportional to m 2 e /m 2 ϕ j ∼ m 2 e /v 2 SM .The integrals of momentum k and the Feynman parameter y have been calculated explicitly.Here, r f,W,G,ηa are mass squared ratios defined as The coefficients C f,W E and C f,W O encapsulate the information from the inner fermion and W boson loops, and they are defined as [37,47]: where for CP-even (CP-odd) scalar, and g W W G = e cos θ W / sin θ W (e) for G = Z(γ).The function j (1,2) 3 (r f,W , r G , r ηa ), arising from the integration of Feynman parameters, is defined by For the case where the gauge boson is the photon, i.e., G = γ (m G = 0), the function j 3 becomes:

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.

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 DS .In addition, we allow for interactions between the scalar S and the Higgs doublets of the 2HDM through the potential V HO .The expressions for V DS and V HO are given by In our analysis, we parameterize the complex scalar as S = (s 1 + is 2 )/ √ 2. To simplify and focus on relevant terms in V DS and V HO , 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. 5here are 6 complex parameters in the scalar potential of the model: m 2 12 and λ 5 in Eq. (2.1) and B S , D S , λ 12S and λ 21S in V DS + V HO .To explore CP violation, it is advantageous to identify rephasing-invariant phases, as studied in refs.[49,50].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 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 θ 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 2 S λ * 12S λ * 21S ] can also induce the same treelevel 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 2 S λ * 12S λ * 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 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 ξ f h,H,A in Table 2 and Eq.(2.8).As a result, we obtain the following new mass eigenstates φ 1,2 : where θ s is a mixing angle with m 2 φ 1 < m 2 φ 2 , and its value is determined by When we assume the portal couplings λ 12S and λ 21S to be real, θ s becomes a phase of B S and corresponds to the physical CP-violating phase.The mass eigenvalues are given by Incorporating all the mass eigenstates into the potentials, we can derive the corresponding neutral scalar cubic and quartic couplings as Table 5: The relevant neutral scalar cubic and quartic couplings.For instance, λ h12 represents the trilinear coupling involving the mass eigenstate of the SM Higgs and the dark sector scalars φ 1 and φ 2 .
where j, k = 1, 2 and each coupling is summarized in Table 5.
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 hA DS and F HA DS , as shown in the diagrams of Fig. 3, Here, we have defined the shorthand notations and Λ UV and µ are a UV divergent parameter and 't Hooft mass parameter, respectively.
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 [51,52] 6 .First, we define the renormalized two-point functions ΣhA,HA (k 2 ) as 6 To deal with the UV divergent in these self-energies, there is another technique for renormalization, so-called "CP-odd tadpole renoramlization" proposed in ref. [53] (see also refs.[54,55] for its application to the supersymmetric models).This methods is basically equivalent to the wave function renormalization, since there is an explicit correspondence in contributions from CP-odd counterterm T A and renormalization constants δZ hA , δZ Ah , δZ HA , δZ AH in Eqs.(4.15) and (4.16) as where each left-hand side is obtained by the CP-odd tadpole renormalization.We found that the UV divergent parts of each contribution are totally same.Note that in this model, we have loop-induced CPodd tadpole contributions by using λ A11 and λ A22 couplings in Table 5, which are UV divergent ones.
In the on-shell subtraction scheme, we can determine δZ hA,HA and δZ Ah,AH by Re ΣhA (m 2 h ) = 0 , Re ΣHA (m 2 H ) = 0 , Re ΣhA,HA (m 2 A ) = 0 .(4.17) By solving these equations, the renormalized two-point functions ΣhA,HA (k 2 ) are given by One can easily check that the divergent parts (and also other constant terms concerning the momentum k) are canceled, and hence, ΣhA,HA give finite results.For numerical predictions, we need to replace F hA DS → ΣhA and F HA DS → ΣHA in the expressions that we have obtained in the previous section.
In the HDS approximation, the dark sector scalar masses m φ 1 , m φ 2 are much larger than the masses in the 2HDM.In this case, we can expand ΣhA(HA) (k 2 ) as where we choose M 2 in Eq. (3.16) to be m φ 1 m φ 2 .Each coefficient is given explicitly as where we have defined the following functions: Since a hA(HA) 0 = ΣhA(HA) DS (k 2 → 0), they are used to calculate the coefficient of the fourfermion operator C f f ′ by inserting a hA,HA 0 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.

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: 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 → inv [62].
Note that the value of θ 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.Fig. 4 shows the size of the electron EDM d e in terms of m φ 1 , with 3 choices of mass difference for m φ 1,2 : ∆m φ ≡ m φ 2 − m φ 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 → φ i φ j , which leads to the exclusion of m φ 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 • q and k • p 1,2 terms which give small corrections of O (m 2 e /v 2 SM ).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 β 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 β 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 β, the value of λ 1S becomes irrelevant for |d e |, while the values of λ 2S and λ − 12S are crucial to enhance |d e |.The value of |d e | is also enhanced by the hierarchy between m φ 1 and m φ 2 , as observed in Fig. 4.This enhancement occurs because ΣhA(HA) (k 2 ) are amplified when m φ 2 ≫ m φ 1 , owing to the logarithmic terms involving the ratio of their masses, which can be achieved by choosing m 2 S ≃ |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 φ 1 , m φ 2 ) = (300 GeV, 1300 GeV), |d e | ∼ 4 × 10 −30 e cm is predicted.Even for the case of ∆m φ = 10 GeV, a large |d e | very close to the current upper limit can be obtained for m φ 1 ≃ 100 GeV.Note that the behavior of d e around m φ i ≃ m h,H,A /2 is originated from the function Re DiscB(m 2 h,H,A , m φ i , m φ i ) , which has following properties: • For m φ i = m/2, Re DiscB(m 2 , m φ i , m φ i ) = 0 and non-differentiable with respect to m or m φ i .

Discussion on the HDS approximation
If m φ 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.21)-(4.24).Fig. 5 shows the numerical estimates of the electron EDM in Type-II 2HDM for ∆m φ = 100 GeV, with (red dotted) and without (solid) the HDS approximation.Similar behaviors can be seen for other types and/or other ∆m φ cases.We can see from the figure that for large m φ 1 , both curves decrease monotonically.However, the values of |d e | are not close to each other.Moreover, even for the mass of m  de | in the left panel of Fig. 6.We can see that the HDS approximation agrees within 1% when m φ 1 ≳ 4.5 × 10 4 GeV.However, as parameterized in Eq. (4.28), the physical prediction is the sum of d Org.From this plot, we can understand that |ϵ de | for the HDS approximation is much smaller than that for the result without the HDS approximation for the heavy m φ 1 region.This means that the cancellation in the HDS approximation is strong, and hence, d HDS e 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 Org.
e , 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.
Accidentally, the HDS approximation is good even for m φ 1 = O(v SM ).This is because around this mass range, the dominant contribution to d e and d HDS e is coming from d Ren.

Fermion dark sector model
Let us consider a dark sector model which contains a vector-like singlet fermion ψ L,R and an SM singlet complex scalar S. The dark sector Lagrangian and the scalar potential include the following terms: where V H represents the scalar potential of the 2HDM sector given in Eq. (2.1). 7We assume a Z 2 symmetry in which H 1 and ψ L,R are odd.The charge assignment is summarized in Table 6.According to this Z 2 assignment, µ 12S and µ 21S terms are soft Z 2 breaking terms.It should be noted that µ 12S and µ 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 0 1 + is 0 2 )/ √ 2 and assume that S does not acquire a nonzero VEV.
The model incorporates numerous complex parameters, encompassing 13 invariant phases as detailed in appendix C.2. Adhering to the presence of CP violation solely where θ O is a mixing angle between A 0 in the 2HDM sector and s 0 2 .On the other hand, R E is parameterized by three mixing angles α, α hS , α HS as With these mixing angles, we can rewrite the explicit CP-violating couplings in Eq. (4.35) in terms of the scalar mass eigenstates, where we have defined ŷR,I ψ,j ≡ y R,I ψ (R E ) 3j and ỹR,I ψ,a ≡ y R,I ψ (R O ) 3a .The SM fermion Yukawa couplings and the W boson couplings are also rewritten as Here, ξ f ϕ j and ξ f ηa are summarized in Table 7, and . With new Yukawa couplings for ψ in Eq. (4.40), the CP-violating contribution to the CP-even scalar ϕ j and CP-odd scalar η a mixing diagrams in Eq. (2.10) can be calculated at one-loop by running the vector-like fermion, where Λ UV denotes the UV divergent parameter, µ is the 't Hooft mass parameter and N ψ c represents some new color factor for the vector-like fermion ψ.In the present study, we set N ψ c = 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 ϕ j -η a self-energy [51,52].Then, the renormalized two-point function Σϕ j ηa (k 2 ) is defined as where no summations over j and a should be understood.In the on-shell subtraction scheme, we can determine δZ ϕ j ηa and δZ ηaϕ j by Re Σϕ j ηa (m By solving these equations, we can obtain the renormalized two-point function Σϕ j ηa (k 2 ) as where the divergent part is properly canceled.
In the HDS approximation, we can further expand Σϕ j ηa (k 2 ) as Here, M 2 in Eq. (3.16) corresponds to the vector-like fermion mass m 2 ψ , and the coefficients are As in the case of the scalar dark sector model, a 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 → ψ ψ.

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: The choices of y R,I ψ can be obtained by |y ψ | = 1 and θ ψ = π/4, and the value of θ ψ corresponds to the one which maximizes the d e with the assumption where all parameters in V scl are real.We then calculate |d e | as a function of m ψ .The C S contribution is found to be significantly smaller than the electron EDM contribution d e , and hence neglected in the subsequent calculations.
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 • q and k • 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 → ψ ψ, especially if the vector-like fermion is light.The branching ratio is directly proportional to |y ψ | 2 sin 2 α hS , making |y ψ | and α hS critical parameters.Larger values of |y ψ | and/or α hS tend to exclude the lower mass region m ψ < m h /2.We have checked that the parameter choice in Eq. (4.50) can avoid the constraint.The curves of the electron EDM |d e | in Fig. 7 show significant complexity, stemming from the m ψ dependence of the renormalized two-point function as given in Eq. (4.46).For the region m ψ < 200 GeV except for m ψ ≃ m ϕ 1 /2 = m h /2, the sharp valleys is due to the sign flipping on d e , while the behavior around m ψ ≃ m ϕ j /2, m ηa /2 is originated from properties of functions Re DiscB(m 2 ϕ j , m ψ , m ψ ) and Re DiscB(m 2 ηa , m ψ , m ψ ) .Notably, the region that |d e | is peaked corresponds to m ψ ≃ m ϕ j /2 or m ηa /2, where Σϕ j η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, Σϕ j ηa (k 2 ) approaches to zero as m ψ goes to zero or infinity.Such a feature is visible in Fig. 7.Note that for m ψ ≳ 7 × 10 3 GeV, the curves of |d e | change their behaviors, due to the low accuracy in the numerical integration, which is explained below.

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 Compared with the scalar dark sector model, the HDS approximation is good, due to the simple expression for Σϕ j ηa in Eq. (4.46) (see Eqs. (4.18) and (4.19) 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 ψ limit.e becomes less than 1%, and hence, the HDS approximation can be used to estimate the value of the electron EDM.However, for m ψ ≳ 7 × 10 3 GeV, the accuracy of d Org. e 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 ψ ∼ 7 ≳ 10 3 GeV, the cancellation between d Org.  de |, the numerical results with and without the HDS approximation deviate for the mass range of m ψ ≳ 7 × 10 3 GeV.As a result, for such a heavy m ψ , 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 ψ ≃ m ϕ j /2 and m η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 HDS

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 [63], 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.in S by appropriately redefining all the relevant couplings.The potential leads to the following minimization condition: 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 where we define and λ β 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 µ 12S and µ 21S , m2 1S = 0 and m2 2S = 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 2 odd 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 2 odd as The remaining part is diagonalized by another mixing matrix, with the mixing angle, Here, ∆m 2 η is the mass-squared difference of two physical CP-odd scalars, Then, each mass eigenstate can be obtained as where R O is in Eq. (4.38).
In contrast to the CP-odd sector, the CP-even sector requires three mixing angles for diagonalization of M 2 even , as shown in Eq. (4.39).The mixing matrix R E is obtained as

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.[49,50].

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 where effective vertices of photon-gauge boson (γ/Z)-scalar Γ µν f,W have information about the inner loop with fermion and W boson, respectively, G represents the photon or Z boson, and Φ is ϕ j or η a .Note that we should sum up all contributions for the final result.For diagrams in Fig. 2, Γ µν f,W can be parameterized by Here, C f,W E and C f,W O are obtained from the inner loop calculation as [37,47] C and ∆ f,W is defined in Eq. (3.14).Similarly, the right diagrams in Fig. 2 are f,W 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): Eq. (D.7): where we omit some propagators for simplicity, and are independent of the loop momentum k, but depend on p 1,2 and q as well as x and y through C f,W E,O .These expressions are generally obtained as

D.3 Decomposed expressions for the electron EDM
In Eq. (4.28), we define two parts of contributions to the electron EDM: F ϕ j ηa DS 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 (ϕ j ) and CP-odd (η a ) scalars via a dark sector CP violation, denoted as F ϕ j ηa DS .If F ϕηa DS 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: Σϕ j ηa = F where we omit the superscript f and W from the general result in Eq. (3.13) for simplicity.Note that for d Ren.
e , 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.For the HDS approximation, we expand F ϕ j ηa DS ((k + q) 2 , M 2 ) and δF ϕ j ηa DS ((k + q) 2 , M 2 ) with respect to k 2 /M 2 .
For the scalar dark sector model, we have

) with k ThO ≃ 1 . 8 ×
10 −15 GeV 2 e cm, where we have used E ThO eff ≈ 78 GVcm −1 for the effective electric field and W ThO S = −282 kHz × 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 HfF + | < 4.1 × 10 −30 e cm at 90% CL, which can be interpreted as |d e + k HfF + C S | < 4.1 × 10 −30 e cm , (3.3) with k HfF + ≃ 1.1 × 10 −15 GeV 2 e cm, where we have used E HfF + eff ≈ 23 GVcm −1 for the effective electric field and W HfF + S = −51 kHz × h for the molecule specific constant.

Figure 1 :
Figure 1: The two diagrams relevant to the CP-violating Lagrangian in Eq. (3.1).Left: The fourfermion 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 ϕ j or CP-odd scalars η a , and the gray shaded blob symbolizes the effective mixing between ϕ j and η a arising from the CP-violating dark sector contribution, which corresponds to F ϕj ηa DS (k 2 , M 2 ) in Eq. (3.6).

Figure 3 :
Figure 3: One-loop diagrams relevant to the effective h-A and H-A mixings induced by φ 1,2 in the scalar dark sector model.

Figure 4 :
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 ∆m φ ≡ m φ2 − m φ1 : ∆m φ = 10 GeV (red), ∆m φ = 100 GeV (green) and ∆m φ = 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 → inv[62].

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

e
and d Org.,HDS e .To quantify the accuracy of the HDS approximation, we now define ρ Org.plot |1 − ρ Org.

e
and d Ren.

e,
and therefore, a cancellation between them may affect the accuracy of the physical results.The right panel of Fig.6shows this feature, by definingϵ de ≡ d e d Org.

Figure 6 :
Figure 6: The left panel indicates the accuracy of the HDS approximation with ρ Org. de ≡ d Org.,HDS e

ϕ j ηa 0 (Figure 7 :
Figure 7: Numerical results of the electron EDM in the fermion dark sector model.The left panelshows 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 → ψ ψ.

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

e
and d Ren.

e
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 ψ , as observed from Fig.7.On the other hand, when m ψ > 1000 GeV, the accuracy of d Org.

e
and d Ren.

e
is also saturated, at the level of |ϵ de | ≃ 10 −5 , while that between d Org.,HDS e and d Ren.,HDS e decreases monotonically.Therefore, even if we obtain a good accuracy for |1 − ρ Org.

e
) is dominated by d Ren. e (d Ren.,HDS e ), as one can see from the right panel of Fig. 8.

28 )
(f,W) • λ L + BR,Φ(f,W) • λ R = −8C 2 p CIt should be emphasized that the effects of k •p i and k •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).
DScorresponds to the term induced by renormalization.Then, the UV divergence is canceled, and hence, we find a finite prediction, Σϕ j ηa = F DS 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 Org.,Ren.

Table 1 :
The four types of Yukawa couplings in the 2HDM, showing the Z 2 charges of d R and e R

Table 2 :
The coefficients between SM fermions and physical neutral scalars.

Table 4 :
We present an example of Z 4 charges for the scalar dark sector model.The exact values of z d 4 and z e 4 for the right-handed d quark and electron depend on the specific type of the 2HDM being considered.

Table 7 :
The coefficients ξ f ϕj and ξ f ηa for SM fermion Yukawa couplings in the fermion dark sector model.