Collapsing domain walls in the two-Higgs-doublet model and deep insights from the EDM

We study the domain wall solutions in the general two-Higgs-doublet model (2HDM) with a CP-violating phase. The 2HDM with the spontaneouse CP violation is found to have domain wall solutions whose tensions are $\mathcal{O}(10^6)\,{\rm GeV}^3$, which are excluded by the Zel'dovich-Kobzarev-Okun bound. With the explicit CP-violating (CPV) terms as the so-called biased term in the scalar potential, domain walls can collapse in the early Universe. The sizes of the explicit CP violation can be constrained from the Big Bang nucleosynthesis. This constraint is converted to the CPV mixing of $\alpha_c$, and is mostly sensitive to the mass splittings between two heavy neutral Higgs bosons. We estimate the possible gravitational wave signals and the electric dipole moment (EDM) predictions due to the domain wall collapsing. It turns out that the peak spectrum of the GW from the domain wall collapsing cannot be probed in any future program. In contrast, the untenable regions with very tiny explicit CPV parameter in the Higgs potential has been partially excluded by the latest electron EDM measurements at the ACME-II and will be further confirmed or excluded by the future ACME-III projection.


Introduction
The Standard Model (SM) of particle physics has been experimentally verified to be successful, with the discovery of the 125 GeV Higgs boson [1,2] in 2012. Yet, the SM itself cannot address the long-lasting puzzle of the baryon asymmetry in the Universe (BAU). Thus, many efforts have been made to realize the three Sakharov conditions [3] by extending the SM. It is quite general that the extensions to the SM lead to larger symmetries. During the phase transitions of the early Universe, different symmetry breaking patterns are likely to produce topological defects. This is known as the Kibble-Zurek mechanism [4][5][6]. The domain walls, which are two-dimensional topological defects, can exist when the 0-th homotopy group associated with the symmetry breaking of G → H is nontrivial, i.e., Π 0 (G/H) = 1. The existence of the domain wall solutions in any new physics can be problematic, since their energy density scale as t −1 with respect to time t. This means once they are formed in the early Universe, they can dominate the energy densities over radiation and matter and thus spoil the standard cosmology. Therefore, such domain wall solutions are usually avoided in the new physics model building by invoking an inflation phase during their formation. Alternatively, one may constrain some parameters in the new physics models such that the possible domain wall tensions are subject to the so-called Zel'dovich-Kobzarev-Okun bound [7].
In this work, we study the possible domain wall solutions arising from the CP symmetries of the two-Higgs-doublet model (2HDM). Originally, the study of the electroweak theory with two Higgs doublets in the scalar sector was motivated to achieve the spontaneous CP violations (SCPV) [8]. However, the existence of two degenerate vacua with opposite-sign CP phases can lead to the domain wall solutions. The domain wall solutions in the 2HDM were previously studied in Refs. [9][10][11], and recently revisited in Refs. [12,13]. With the typical energy scale of few hundred GeV for the heavy scalars, the tensions of the stable domain walls in the 2HDM are found of the size σ ∼ O(10 6 ) GeV 3 . Therefore, this is excluded by the Zel'dovich-Kobzarev-Okun bound. One possible way to evade the domain wall problem is to include a small symmetry breaking term, that is the so-called biased term in the Higgs potential. In this perspective, the symmetry of the CP transformation is approximate. With such a mechanism [14][15][16], the possible domain walls were unstable and collapsed before they overclose the Universe. Indeed, such biased terms are nothing but the explicit CPV (ECPV) parameters in the 2HDM potential.
With the ECPV parameters as the biased terms in the 2HDM potential, one may expect gravitational wave (GW) signals associated with the domain wall collapses. This was previously discussed in Refs. [17][18][19][20][21][22]. Therefore, it is natural to ask if such GW signals can be probed at the future programs, such as the satellite-based interferometers of LISA [23,24], Taiji [25], and Tianqin [26], or the radio telescope of square kilometer arrays (SKA) [27] and the Japanese space GW antenna (DECIGO) [28]. Note that such GW signals are different from those produced due to the bubble collisions during the strongly first-order phase transitions. The constraint on the ECPV parameters is imposed so that the corresponding biased terms are sufficiently large to collapse the domain walls before the epoch of the the big bang nucleosynthesis (BBN). In Ref. [21], we found such signals can be sufficiently strong when the corresponding new physics scale is above about 10 TeV. However, as we shall show, the related GW peak spectrum here is typically O(10 −23 ) and below the future search sensitivities of SKA [27] and DECIGO [28]. Alternatively, the precise measurements of the electric dipole moments (EDM) play an intriguing role in constraining the small ECPV parameters. The ongoing and upcoming EDM measurements include the electron EDM (eEDM) from the ACME collaboration [29,30], the mercury EDM [31], and the radium EDM [32]. The EDM measurements are usually interpreted to constrain CPV mixings in the context of the 2HDM [33][34][35][36][37][38][39][40][41], or the SUSY-breaking scale [42]. Usually, the constraint on CP-violating 2HDM from EDM was performed in the absence of spontaneous CP violating phase to avoid the problematic SCPV domain walls [34]. With the future improvements of the eEDM precision measurements from the ACME-III, given a few hundred GeV heavy scalar masses, we find that they can be used to exclude sizable regions of the ECPV parameters in the 2HDM potential jointly with the BBN constraint. In this sense, the EDM experiments provide us deep insights to the very tiny ECPV term that triggered the domain wall collapse in the early Universe. This is completely different from what was usually discussed in the context of SUSY, such that the precision measurements of the EDM will constrain the corresponding new physics scale above O(10) TeV [42].
The rest of the paper is organized as follows. In Sec. 2, we review the general 2HDM. We focus on the minimization conditions and the mass spectrum in both the SCPV and the ECPV scenarios. The parameters between the generic basis and the physical basis are built for solving the domain wall. The theoretical constraints of the perturbative unitarity and the vacuum stabilities are imposed to the scalar self couplings. In Sec. 3, we present the domain wall solutions in the CPV 2HDM. The domain walls from the SCPV typically have tensions of O(10 6 ) GeV 3 , which should be collapsed with sufficiently sizable biased terms. The ECPV parameter in the 2HDM, which we choose to be Imλ 5 , will be bounded from below by considering the BBN constraint. By estimating the corresponding GW signals from the domain wall collapses, they are beyond the search sensitivities of any future satellite-based observation programs. In Sec. 4, we turn to the estimations of the eEDM in the domain wall collapsing scenario. The current and the future projections of the eEDM measurements from the ACME are used to set upper limits to the CPV mixing angle, as well as exclude the ECPV parameter directly. We summarize our results in Sec. 5.

The general 2HDM with the CPV
In this section, we review the general CP-violating 2HDM and focus on the global symmetry from the CP transformation of the potential. Previous studies of the symmetries and topological structures of the 2HDM include Refs. [12,34,[43][44][45][46][47][48][49].

(2.2)
In order to obtain the domain wall solutions below, we treat (ϕ 1 , ϕ 2 , Θ) as background fields. When the electroweak symmetry breaking (EWSB) occurred, they obtain the vacuum expectation values (VEVs) as ϕ 1 ,2 = v 1 ,2 , Θ = θ, with θ to be solved below for both the SCPV and the ECPV scenarios. As usual, two VEVs are parametrized as v 1 = vc β and v 2 = vs β 1 . The minimizations of the most general 2HDM poential (2.1) are the following In the special case of vanishing spontaneous CPV phase θ = 0, the minimization conditions are reduced to In this work, we focus on the SCPV1 scenario defined in Ref. [49], and use the notion of SCPV throughout the context. Thus, it is equivalent to take 2δ 1 −δ 2 = nπ, with n being an integer. Without loss of generality, we shall take n = 0 for our later discussions. In order to simplify the calculation, we make a rephase of Φ 2 → e −iδ 1 Φ 2 to eliminate the phases in the coefficients, so that Imm 2 12 = Imλ 5 = 0. The minimization conditions for the SCPV scenario are then expressed as Two degenerate solutions of with two opposite signs lead to the domain wall solutions. Note that this expression of θ is not our solution in the physical basis. With the minimization conditions of the SCPV scenario in Eqs. (2.7), the 2HDM potential can be expressed as This is consistent with what was obtained in Ref. [12].

The mass spectrum and self couplings with the SCPV
For the SCPV scenario, we have the charged Higgs boson mass squared matrix of and the corresponding eigenvalues are The gauge eigenstates of (H ± 1 , H ± 2 ) are transformed into mass eigenstates of charged Nambu-Goldstone bosons G ± and Higgs bosons H ± by Similarly, two pseudoscalars of (A 1 , A 2 ) are transformed into neutral Nambu-Goldstone bosons G 0 and a pseudoscalar A 0 by In the basis of (H 1 , H 2 , A 0 ), the 3 × 3 mass squared matrix for the neutral states are diagonalized by [34,51,52] (2.14c) Each term in the mass squared matrix of M 2 0 is listed below In turn, the quartic scalar self couplings are expressed in terms of the Higgs boson masses and mixings as There is a well-known constraint from the (M 2 0 ) 13 and (M 2 0 ) 13 terms in Eq. (2.15): This leads to one additional constraint between mixing angles and mass eigenvalues as follows [53] ( .

(2.18)
In practice, we use three mixing angles of (α , α c , β) as inputs. With the special limit of and there is a singularity at t β = 1 under this limit. In the small α c limit, we can further approximate this relation as In Fig. 1, we display the values of |α b | versus the input CPV mixing angle of α c ∈ (10 −10 , 1) with various mass splittings of ∆M = m 3 − m 2 between two heavy neutral Higgs bosons.
The decreasing values of α c also lead to smaller values of |α b |. Also, the sizes of the CPV mixing angles are suppressed with very degenerate mass splitting of ∆M , as we vary ∆M from 100 GeV to 1 MeV in the plots. This pattern can be found straightforwardly from Eq. (2.20) as well. Therefore, the joint effects of smaller input parameter α c and the smaller mass splitting of ∆M can lead to suppressions to the eEDM. Note that the relations of α b versus α c as given in Eq. (2.18) and Fig. 1 are independent of the quartic scalar self coupling terms in the 2HDM potential, thus they hold for both the SCPV and the ECPV scenarios. Since the α b dependences on α c and the mass splitting ∆M are close between the t β = 0.5 and t β = 2 cases, we shall always use the t β = 0.5 as our benchmark throughout our discussions below. In the physical basis, the SCPV phase of θ is obtained from the relation of One may expressμ 2 1A andμ 2 AA in terms of masses and mixing angles explicitly as followŝ By combining with the dependences of |α b | on α c as depicted in Fig. 1 under the β − α = π/2 and small α c limit. Here, we have also made used of Eq. (2.20). Clearly, one can envision thatμ 2 1A μ 2 AA . Accordingly, the SCPV phase becomes Obviously, the solutions of the SCPV are approaching to ± π 2 when the CPV mixing angle α c is suppressed and/or two heavy neutral Higgs bosons are very mass degenerate.
We summarize the parameter inputs for the SCPV scenario in both the physical basis and the generical basis in Table. 1. The CPV mixing angle α b and the SCPV phase θ can be obtained by Eq. (2.18) and Eq. (2.21), respectively. Therefore, one has eight independent parameters in both basis.

Generical basis
Physical basis Table 1. The parameter inputs in both the generical basis and the physical basis for the SCPV scenario.

The mass spectrum and self couplings with the ECPV
The most general 2HDM potential with complex parameters of m 2 12 and λ 5 violates the CP symmetry explicitly. The corresponding minimization conditions were previously given in Eqs. (2.3).
For the mass spectrum, the same conventions of the mass squared matrix and mixing angles are adopted as in the SCPV scenario. The charged Higgs boson masses are expressed as Schematically, the neutral mass squared matrix is the same as in Eq. (2.15a), with each element listed as beloŵ From the expressions ofμ 2 1A andμ 2 AA , one finds the relation between the imaginary components of Imm 2 12 and Imλ 5 aŝ which will be used for solving the relative CPV phase of θ in general when replacingμ 2 AA andμ 2 1A by physical inputs according to Eq. (2.14). In the ECPV scenario, the phase transformation of Φ 2 can remove one of the three phases. The third minimization condition in Eq. (2.3c) will help to relate the two remaining phases. Thus, we have only one free CPV phase again. There are several situations in solving the CPV phase of θ from Eq. (2.27) for the ECPV scenario: 1. Imm 2 12 = 0, the above equation is a quadratic equation of t θ . This is equivalent to rephase Φ 2 as Φ 2 → e −iδ 1 Φ 2 .
2. Imλ 5 = 0, the above equation is a quartic equation of s θ . This is equivalent to rephase 3. If one keeps both Imm 2 12 and Imλ 5 non-zero, the above equation is a quartic equation of t θ/2 .
To simplify our discussions, we take the special case of Imm 2 12 = 0, where t θ satisfies The exact solution for t θ is expressed as Given that we haveμ 2 1A μ 2 AA in the small α c limit, the ∓ signs in the exact solutions in Eq. (2.29) correspond to the limits of θ → ±π/2 and θ → 0, respectively. To restore the exact solution of t θ in the SCPV scenario, one should take the − sign in Eq. (2.29). It turns out that the relative CPV phase θ will enter into the biased term to make the domain wall collapse in the form of s 2θ . In both limit of θ → ±π/2, we have s 2θ → 0.
Besides of the CPV mixing angle of α c , the ECPV parameter of Imλ 5 also plays a role. In the limit of |Imλ 5 with the leading contributions the same as the expression in the SCPV scenario Eq. (2.21). Correspondingly, the relative CPV phase of θ is expected to depend on CPV mixing angle α c significantly. On the other hand, in the limit ofμ 2 AA |Imλ 5 | μ 2 1A /s β , we have the approximation to t θ as As we have shown above thatμ 2 AA → m 2 3 /v 2 in the small α c limit, the value of θ will be independent of the CPV mixing angle α c in such a limit. In Fig. 2, we depicted |s 2θ | versus the CPV mixing angle of α c ∈ (10 −10 , 1) with various inputs of Imλ 5 . Indeed, for relatively sizeable values of Imλ 5 = 10 −1 and Imλ 5 = 10 −10 , the solutions of |s 2θ | are plateaued when α c drops to certain threshold. Particularly for the Imλ 5 = 10 −1 and ∆M = 1 MeV, the values of s 2θ is basically invariant with the varying α c . When the ECPV parameters are suppressed to Imλ 5 = 10 −15 and even Imλ 5 = 10 −20 , one finds that values of |s 2θ | are (almost) always decreasing with respect to the α c inputs in the ranges of our consideration. Together, they justify our approximations in Eqs. (2.30) and (2.31). With the solution of θ obtained in Eq. (2.29), we can solve for the quartic scalar self couplings in terms of the Higgs boson masses and mixings as below with the choice of Imm 2 12 = 0. All quartic scalar self couplings can easily return to those in the SCPV scenario as listed in Eqs. (2.16), by taking Imλ 5 = 0. The Higgs potential with the ECPV scenario can be extended from Eq. (2.9) as below Generical basis Physical basis λ 1 ,2 ,3 ,4 , Reλ 5 , Imλ 5 m 1 ,2 ,3 , m ± , v m 2 11 , m 2 22 , Rem 2 12 , Imm 2 12 α , α b , α c , β , θ , Imλ 5 , Imm 2 12 Table 2. The parameter inputs in both the generical basis and the physical basis for the ECPV scenario.
We summarize the parameter inputs for the ECPV scenario in both the physical basis and the generical basis in Table. 2. The CPV mixing angle α b is obtained from Eq. (2.18) from other masses and mixing angles. The relative CP phase can be obtained from Eq. (2.29). As we have argued previously, one can always remove one of two ECPV parameters by rephasing the second doublet Φ 2 . Without loss of generality, we choose to take Imm 2 12 = 0. Hence, one has nine independent input parameters in both basis. We shall describe the procedures of converting the input parameters in the physical basis to the generic basis.
1. In the physical basis, all necessary input parameters are masses of three neutral Higgs bosons m 1 ,2 ,3 , the mass of charged Higgs bosons m ± , the Higgs VEV of v = ( √ 2G F ) −1/2 246 GeV, and three mixing angles of (α , α c , β), as well as the ECPV parameter Imλ 5 . The state of h 1 will be regarded as the SM-like Higgs boson, with m 1 = 125 GeV. To simplify the discussion, we shall take the special limit of α = β− π 2 .
3. The relative CPV phase of θ is derived from the exact solution in Eq. (2.29) with all above parameter plus the ECPV parameter of Imλ 5 .
4. With all Higgs masses of m 1 ,2 ,3 ,± , mixing angles of (α , α b , α c , β), the relative CPV phase of θ, and the ECPV parameter of Imλ 5 , we obtain the quartic scalar self couplings from Eqs. (2.32), and mass squared parameters of (m 2 11 , m 2 22 , Rem 2 12 ) from the minimization conditions in Eqs. (2.3). These parameters in the generical basis will be used for solving the domain wall profiles and determine the sizes of the biased terms for the domain wall collapse.

The unitarity and stability constraints
It is well-known that the perturbative unitarity and stability constraints to the Higgs potential should be imposed. The unitarity bounds of the 2HDM were previously studied in Refs. [54,55]. The perturbative unitarity constraint means that the theory cannot be strongly coupled. In practice, the necessary and sufficient condition of the tree-level unitarity bounds can be obtained by evaluating the eigenvalues of the S-matrices for the scattering processes of the scalar fields in the 2HDM [54,55]. Due to the Nambu-Goldstone theorem, the S-matrices can be expressed in terms of 2HDM quartic scalar self couplings λ i . The S-wave amplitude matrices are due to fourteen neutral, eight singly-charged, and three doubly-charged scalar channels in the 2HDM. They read neutral a 0 0 : The S-wave amplitude matrices for three different channels are expressed as The expressions for the submatrices of (X 4×4 , Y 4×4 , Z 3×3 ) are given as follows The unitarity requires |a i 0 | ≤ 1. The stability constraints require a positive 2HDM potential for large values of Higgs fields along all field space directions. Collectively, they read In Fig. 3, we present the theoretical constraints in the (m 2 , t β ) plane, with other parameters fixed to be ∆M = 1 GeV and α c = 0.1. The heavy neutral Higgs boson masses are found to be bounded from above, and large t β 5.0 or small t β 0.2 are disfavored.

The Yukawa couplings in the CPV 2HDM
We focus on the CPV 2HDM where the Yukawa sector has a Z 2 symmetry and Φ 1 and Φ 2 each only gives mass to up-type quarks or down-type quarks and charged leptons. This is sufficient to suppress tree-level flavor changing processes mediated by the neutral Higgs bosons. The Yukawa couplings for the Type-I and Type-II 2HDM read (and suppressing the CKM mixing), where Q T L = (u L , d L ) andΦ 2 ≡ iσ 2 Φ * 2 . For both cases, the charged lepton Yukawa coupling has the same form as that of the down-type quarks. Therefore, we can express the couplings between neutral Higgs bosons and the fermions and gauge bosons in the mass eigenbasis When c f,icf,i = 0 or a icf,i = 0, the mass eigenstate h i couples to both CP-even and CPodd operators, so the CP symmetry is violated. The coefficients of c f,i ,c f,i and a i can be derived from the elements of the rotation matrix R defined in Eq. (2.14), which were also previously obtained in Refs. [33,34,36]. Here, we summarize their explicit expressions in Table. 3. In the special limit of β − α = π/2, the Yukawa couplings and Higgs gauge couplings are determined by the CPV mixing angles of (α b , α c ) and t β . By taking the CPC limit of α b = α c = 0, it is evident that (h 1 , h 2 ) have the purely CP-even Yukawa couplings of c f ,i , while h 3 has the purely CP-odd Yukawa couplings ofc f ,i . The previous studies of the collider measurements of the CPV in the Higgs Yukawa couplings can be found in Refs. [56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73].

Domain walls in the CPV 2HDM
In this section, we study the domain wall solutions in the 2HDM with the SCPV vacuum solution 2 . Such solutions arise from the CP transformations of Φ j → Φ * j in the 2HDM.

The domain wall solutions
Under the discrete CP transformations of two Higgs doublets one has Θ → −Θ for the background fields. This means the CP transformation to two Higgs doublets is equivalent to a Z 2 transformation to their relative phase. Correspondingly, the CPC part of Eq. (2.1) is invariant, while the CPV terms in Eq. (2.1) is manifestly odd. The vacuum manifold and the corresponding nontrivial homotopy group [48,74] is where one uses the fact that the CP symmetry is homeomorphic to the Z 2 symmetry, and the vacuum manifold of SU(2) L × U(1) Y /U(1) em is homeomorphic to S 3 . Therefore, the SCPV part of the 2HDM potential leads to a domain wall solution.
The domain wall solution is obtained in the 'Euclidean basis' of Two domains correspond to φ = (v 1 , c θ v 2 , ±s θ v 2 ). The energy density is given as follows: where V 0 is a pure constant to make the potential of electroweak vacuum zero. In our real calculation of the domain wall profile, V ECPV should be taken into account, as it shifts the potential as well as the local minima positions. However, when considering the case |Imλ 5 | 1, we can safely ignore the ECPV part when solving the domain wall profile. The tension of the domain wall is then the integral of the total energy density: The corresponding domain wall profile is obtained by solving the equations of motion (EOM) of φ: with the boundary conditions being The EOM is solved using the path deformation algorithm 3 [21,75]. In Fig. 4, we display a sample of the domain wall profiles and the energy density of E(z). Since we expect that the relative CPV phase to be θ → ±π/2, thus, the imaginary component of ϕ 2 takes much larger value comparing to the real component. In our numerical estimation, we find that the domain wall tensions from Eq. (3.5) are typically ∼ O(10 6 ) GeV 3 . One can quickly find this result from the energy density plot of our sample in Fig. 4, where the local energy density is roughly E ∼ O(10 8 ) GeV 4 and the domain wall width is around ∼ O(0.01) GeV −1 .

The biased term from the ECPV
The ECPV component of the potential leads to additional contribution to the energy density as below which becomes the biased term for the SCPV domain walls. After the EWSB, the corresponding biased term reads with the relative CPV phase of θ being solved from Eq. (2.29).

The cosmological constraints
The observation of the CMBR leads to the following condition to the domain wall tension which was known as the Zel'dovich-Kobzarev-Okun bound. Our numerical solutions found that the domain wall solutions in the SCPV case generally lead to tensions of σ ∼ O(10 6 ) GeV 3 , which is ∼ 10 15 times above the Zel'dovich-Kobzarev-Okun bound. Therefore, the biased term from the ECPV is necessary for the domain wall collapse.
To have the large scale domain wall to form, it is known that the energy difference between two vacua should be sufficiently small: with |V 0 | representing the height of the potential barrier between two minima in Eq. (3.4f), and the critical value of p c = 0.311 predicted from the percolation theory [76]. Lower bounds can be also obtained to the energy difference. The domain wall cannot exist too long to spoil the known constraints from the BBN [19,77,78]. This leads to an lower bound to the energy difference ∆V 1/4 5.07 × 10 −4 GeV C 1/4 with A ∼ 0.8 ± 0.1 [79], and C ann = 2 for the Z 2 symmetry.σ ≡ σ/(1 TeV) 3 denotes the dimensionless domain wall tension. Assuming that the domain wall collapse occurred during the radiation dominated era, the corresponding temperature is given by with g * (T ann ) counting the relativistic degrees of freedom at the annihilation temperature, and ∆V ≡ ∆V /(1 MeV) 4 . In Fig. 5, we present the size of the biased term ∆V versus the CPV mixing angle α c ∈ (10 −10 , 1), for different mass splittings of ∆M = (100 GeV , 10 GeV , 1 GeV , 1 MeV) between two neutral heavy Higgs bosons. The upper bound from Eq. (3.11) and lower bound from Eq. (3.12) are presented in terms of dotted (in purple) and dashed (in black) lines, respectively. It turns out that the upper bound can be always satisfied with the parameter choices of α c in our considerations. Meanwhile, the lower bounds to the domain wall collapse can become sensitive to the CPV mixing angle of α c only when the explicit CPV parameter of Imλ 5 is sufficiently small. This can be expected from our previous discussions about the t θ dependences on the physical inputs of α c and Imλ 5 . For the Imλ 5 = 10 −10 case, one also finds that ∆V becomes plateaued similar to the corresponding |s 2θ | plot in Fig. 2. On the other hand, too small ∆M and/or Imλ 5 values are excluded by the BBN lower limit as shown in the lower panels of Fig. 5.

The GW signals
The collapsing domain walls can lead to GW signals [17-19, 79, 80], while we find that such signals arising from the CPV 2HDM are impossible to be probed in any of the future satellite observations. We shall briefly discuss the signal estimations below.
The peak frequency of the GWs at the annihilation time of domain walls is proportional to the annihilation temperature T ann in Eq. (3.13), and is given by Here, g * (T ann ) and g * s (T ann ) count the relativistic degrees of freedom contributing to the energy density and the entropy density. They are both 10.75 for 1 MeV T ann 100 MeV. For GWs with peak frequencies in the range of O(10 −4 ) − O(10 −1 ) Hz, they may be probed by the future satellite-based interferometers, such as the LISA [23,24], Taiji [25], and Tianqin [26] programs. The GWs with very small peak frequencies of few nano Hz may be probed at the future radio telescope of SKA [27] and the DECIGO [28] with the latter having wider range of typical frequencies of ∼ O(0.1) − O(10) Hz. With the lower limit of the ∆V in Eq. (3.12), there is a lower limit to the peak frequency of f peak 0.94 × 10 −9 Hz g * (T ann ) 10 1/4 g * s (T ann ) 10 Thus, the peak frequencies of the GW signals are expected to be higher than order of several nano Hz. The peak energy density spectrum of the GW is with˜ GW 0.7 ± 0.4 in the scaling regime [79]. By using the annihilation temperature in Eq. (3.13), the peak energy density spectrum becomes By taking the lower limit of the ∆V in Eq. (3.12) into account, we find an upper limit to the peak energy density spectrum as Ω peak GW h 2 (t 0 ) 1.36 × 10 −17˜ GW A 2 g * s (T ann ) 10 −4/3 g * (T ann ) 10 σ 2 .  The SM estimations of the eEDM were of size |d e /e| ∼ 10 −44 cm from the four-loop contributions of the CKM phase, and |d e /e| ∼ 10 −38 − 10 −39 cm by considering the CP-odd electron-nucleon interaction [81][82][83][84].  The dominant contributions to the eEDMs come from the two-loop Barr-Zee diagrams [85] in the CPV 2HDM. There are three types of dimension-five operators involved: (i) the h i γγ operator, (ii) the h i γZ operator, and (iii) the H ± W ∓ γ operator. Expressed in terms of the Wilson coefficient of these operators, we summarize the total contributions as follows and the corresponding diagrams are depicted in Fig. 6. Explicit expressions of these Wilson coefficients can be found in the appendices of Refs. [34,40]. By combining the total contributions in Eq. The Wilson coefficients are generally related to the normalized scalar or pseudoscalar Yukawa couplings as The pseudoscalar couplings are all proportional to the CPV mixing angles asc f ,i ∝ α c , which can be found in Tab. 3 and the relation of Eq. (2.20). One can expect that future improvements of the eEDM precisions by an order of magnitude will further constrain the CPV mixing angle by an order of magnitude. We display the eEDM predictions of |d e /e| versus the CPV mixing angle of α c in Fig. 7. The latest upper bound to the eEDM from the ACME-II from Eq. (4.1) and the future projected upper bound from the ACME-III in Eq. (4.2) are displayed in horizontal dashed lines. The evaluations of the eEDM depend on the mixing angles of (α , β , α b , α c ). By using the constraint for masses and mixing angles in Eq. (2.18), the mass splitting of The eEDM predictions versus the CPV mixing angle α c in the range of α c ∈ (10 −10 , 1). The lower limits to α c are shown in arrows for the Imλ 5 = 10 −15 (upper panels) and Imλ 5 = 10 −20 (lower panels) cases from the BBN bound to the biased term, as obtained from Eq. (3.12). For each case, we display the eEDM predictions for different mass splittings of ∆M = (100 GeV , 10 GeV , 1 GeV , 1 MeV).
∆M can also play a role in the size of the eEDM predicitons. With a fixed input of α c , one generally has a suppressed value of α b with smaller inputs of ∆M , as was displayed in Fig. 1. Concequently, the eEDM predictions will be suppressed as well. This was previously discussed in Ref. [37]. The evaluations of the eEDMs are independent of the size of ECPV parameter Imλ 5 , as one can visualize between two upper panels and two lower panels in Fig. 7. Meanwhile, different inputs of Imλ 5 lead to different cosmological constraints to the α c through Eq. (3.12). Explicitly, these lower bounds to α c are denoted by dashed vertical lines with arrows in each plot. For the Imλ 5 = 10 −15 case, there are lower limits to the CPV mixing angle of α c for all four ∆M inputs. However, when such lower limits of α c are saturated, the corresponding eEDM predictions are |d e /e| ∼ O(10 −35 ) cm, which is another five orders of magnitude below the future precision aimed at the ACME-III. When one further reduces the ECPV parameter to Imλ 5 = 10 −20 , the constraint of Eq. (3.12) has already ruled out the situation with small mass splitting of ∆M = 1 MeV. This can be also observed in the lower-left panel of Fig. 5. Thus, we denote the ∆M = 1 MeV with Imλ 5 = 10 −20 cases by dashed lines, indicating that their exclusion from the BBN constraint. We also checked that for the Imλ 5 = 10 −20 case, the lower limits of α c from Eq. (3.12) lead to the eEDM predictions of |d e /e| ∼ O(10 −30 ) cm. Therefore, the situation with very tiny ECPV parameter of Imλ 5 = 10 −20 is expected to be confirmed or excluded with the joint BBN constraints and the improved experimental precision of the eEDM from the future ACME-III.
In Fig. 8, we further present the joint BBN constraint from Eq. (3.12) and the eEDM measurements from the current ACME-II limits (4.1) (orange shaded regions) and the future ACME-III projections (4.2) (vertical dashed lines). Two different mass splittings of ∆M = 100 GeV (left panels) and ∆M = 1 GeV (right panels) are displayed. With the future improvements of the eEDM precision by an order of magnitude, the corresponding constraints to the CPV mixing of α c will be improved by an order of magnitude accordingly. For a relatively large mass splittings of ∆M = 100 GeV, the upper limits to α c from the future eEDM measurements can be as small as ∼ O(10 −4 ). While for a suppressed mass splittings of ∆M = 1 GeV, the upper limits to α c become ∼ O(10 −2 ). This is due to the relation between |α b | versus α c as given in Eq. (2.20) and Fig. 1.
Furthermore, with the BBN constraints to the biased domain wall term, the sizes of the ECPV parameter Imλ 5 are constrained with various inputs of the physical CPV mixing angle α c . Such constraints are becoming more stringent with smaller inputs of α c . For a fixed value of Imλ 5 , a lower limit to α c is given in the light green region. We find that Imλ

Conclusion and Discussion
In this work, we focus on the vacuum with the SCPV in the 2HDM, which can lead to a domain wall structure. Through our numerical studies, such domain walls typically lead to incredibly large tensions of O(10 6 ) GeV 3 , which is well above the Zel'dovich-Kobzarev-Okun bound. Therefore, the complex parameters in the 2HDM potential are necessary in playing the role as the biased terms to collapse these domain walls. With reasonable sizes of the biased terms, such domain walls could have been formed in the early Universe. In order not to spoil the BBN constraints, we find the direct constraints to the sizes of the ECPV terms, for which we choose to be Imλ 5 , to be O(10 −14 − 10 −12 ) ( O(10 −24 − 10 −22 ) with the CPV mixing angles of α c ∼ O(10 −10 ) (α c ∼ O(1)). Although the related domain wall collapse does not lead to sufficiently strong signals for the future GW probes at the SKA, we find that the eEDM measurements will play a role to probe the deep echos from this process. The future projection of the eEDM measurements from the ACME-III can set upper limits to the CPV mixing angle of α c O(10 −2 ) − O(10 −4 ), which depends on the types of the Yuakwa couplings and mass spliting of two neutral heavy Higgs bosons. For the first time, we find the current and future eEDM measurements have excluded or can be used to probe the very tiny regions of the ECPV parameter of Imλ 5 ∼ O(10 −24 )−O(10 −20 ) (Type-I) or Imλ 5 ∼ O(10 −24 )−O(10 −21 ) (Type-II). In other words, we find that the eEDM measurements can look deep into the possible domain wall collapses in the early Universe. This is different from the discussions in the context of SUSY where the improved eEDM measurements were thought to set lower limit to the SUSY-breaking scale. In the context of the domain wall collapsing triggered by the complex parameters in the 2HDM potential, the future improvements of the eEDM measurements are found to further set up upper limits to the CPV mixing angles and the sizes of the ECPV term.
Some future efforts can be envisioned from this work.
1. The 2HDM, along with other new physics models, are known to produce other topological defects besides of the domain wall solutions from the CP symmetry, such as vortices and monopoles [48,51,74,[86][87][88][89]. The solutions to these structure can lead to GW signals as well. Therefore, it will be useful to probe the complete spectrum of the GWs for given new physics model.

2.
We did not consider the electroweak phase transition and the possibility of achieving the BAU within this frame. There have been some recent progresses [90][91][92] along this direction. It is therefore to ask if a successful BAU can be achieved when the new physics models have non-trivial topological solutions.
3. From the perspective of the EDM measurements, it will be greatly useful to perform the estimations to the atomic EDMs and find the future experimental search projections as well.