Higgs pair productions in the CP-violating two-Higgs-doublet model

In this work, we study the SM-like Higgs pair productions in the framework of the general CP-violating two-Higgs-doublet model. Several constraints are imposed to the model sequentially, including the SM-like Higgs boson signal fits, the precise measurements of the electric dipole moments, the perturbative unitarity and stability bounds to the Higgs potential, and the most recent LHC searches for the heavy Higgs bosons. We show how are the CP-violating mixing angles related to the Higgs cubic self couplings in this setup. Afterwards, we estimate the cross sections of the future LHC/SppC searches for the Higgs pair productions, as well as other possible decay modes for the heavy Higgs bosons.


Introduction
The discovery of the 125 GeV Higgs boson [1,2] at the LHC runs at 7 ⊕ 8 TeV validate Higgs mechanism for the spontaneous breaking of the electroweak gauge symmetry (EWSB). The current LHC measurements of the Higgs boson couplings to the SM fermions, gauge bosons, and loop-induced couplings to photons and gluons reach the precision of ∼ 10 − 20 % level. Besides, it is important to probe the Higgs self couplings to confirm the mechanism of the EWSB. This can be done by looking for the Higgs pair productions at both high-energy e + e − and pp colliders. The current LHC searches for the Higgs pair productions focus on the leading production channel of gluon-gluon fusion (ggF), as well as the promising final states review the setup of the CPV 2HDM. With the assumptions of the degenerate heavy Higgs boson mass spectrum, we take the simplified parameter sets of α = −π/4. We also obtain the gauge couplings, Yukawa couplings, and the self couplings for Higgs bosons in the physical basis. In Sec. 3, we impose series of constraints to the CPV 2HDM parameter space. The combined constraints of 125 GeV Higgs signals and the eEDM bounds point to the t β ∼ 1 parameter choice. The size of the CPV mixing angle |α b | is also bounded from above. For the CPV 2HDM-I, the CPV mixing is stringently constrained to be |α b | 5×10 −3 , which is quite approaching to the CPC limit. For the CPV 2HDM-II, the constraints to the CPV mixing are much relaxed, and we focus on this case for the Higgs pair productions. The constraints from the unitarity, the stability, and the current LHC 8 TeV searches for the heavy Higgs bosons further restrict the allowed mass ranges of the heavy Higgs bosons and the soft Z 2 -breaking mass term of m soft . The main results of the Higgs pair productions in the CPV 2HDM are presented in Sec. 4. By combining the current constraints, we show that the variations of the Higgs cubic self couplings are controlled by the size of the CPV mixing angle |α b | and the soft mass term m soft in the 2HDM potential. A set of benchmark models are given with the fixed CPV mixing angles and the maximally allowed soft mass terms. Under the small CPV limit, the Higgs cubic self coupling of λ 111 for the SM-like Higgs boson tends to the SM predicted value of λ SM hhh 32 GeV, and the resonance contributions become negligible as well. The corresponding Higgs pair production cross sections will tend to the predictions for the SM case. We estimate the physical opportunities of the precise measurement of the SM-like Higgs cubic self coupling λ 111 at the future high-energy e + e − colliders, with focus on the e + e − → hhZ process at the √ s = 500 GeV run. On the other hand, the heavy resonance contributions to the Higgs pair productions can become dominant at the pp colliders. The cross sections for the possible experimental search modes of h 1 h 1 → (bbγγ , bbW W ) are estimated for both LHC 14 TeV and SppC/Fcc-hh 100 TeV runs. In addition, several other possible search modes of (W + W − , ZZ , hZ) are also mentioned. The conclusions and discussions are given in Sec. 5.

1b)
with (m 2 12 , λ 5 ) being complex and all other parameters being real for the CPV 2HDM. After the EWSB, two Higgs doublets Φ 1 and Φ 2 in the unitarity gauge can be expressed as The ratio between two Higgs VEVs is parametrized as and ξ represents the relative phase between two Higgs doublets. The imaginary components of m 2 12 and λ 5 are the source of CP violation, which lead to the mixings among three neutral states as ( The angle α parametrizes the mixing between two CP-even states of (H 0 1 , H 0 2 ). The CPV mixing angles of α b and α c parametrize the CP mixings between (H 0 1 , A 0 ) and (H 0 2 , A 0 ), respectively. Their ranges are taken as (2.5) In the CPC limit, one has α b = α c = 0. Correspondingly, R becomes block diagonal, and (h 1 , h 2 ) are purely CP-even states. By minimizing the CPV 2HDM potential, one obtains the following relations for the mass parameters The physical masses of (M 1 , M 2 , M 3 , M ± ) in the scalar spectrum are obtained from the 2HDM potential together with the minimization conditions given in Eqs. (2.6). The charged Higgs boson mass squared reads The mass squared matrix for the neutral sector can be expressed as with the short-handed notations of By diagonalizing the mass squared matrix with the mixing matrix in Eq. (2.4), one has from which one further obtains the relations to trade the quartic Higgs self couplings into the physical inputs as follows For simplicity, we can always work in the basis where ξ = 0 by using the rephasing invariance. We also assume that Re(m 2 12 ) ≥ 0, and use the notation for the soft mass term as The elements of (M 2 0 ) 13 and (M 2 0 ) 23 in Eq. (2.8) provide the CPV mixings, which are related via t β as This leads to one additional constraint between mixing angles and mass eigenvalues as follows [50] ( (2.14) In the analysis below, we always identify h 1 as the SM-like Higgs boson with mass of 125 GeV. We further simplify the parameter inputs by requiring all heavy Higgs boson masses are degenerate, i.e., M 2 = M 3 = M ± ≡ M . This was usually taken to relax the constraints from the electroweak precision measurements. The constraint of Eq. (2.14) among the mixing angles becomes Below, we will always take α = −π/4. 3 The input parameters of (β , α b ) will be determined through other constraints. Since α c determines the size of the CPV mixing between two mass-degenerate Higgs bosons of h 2 and h 3 in our setup, one can anticipate that α c becomes unphysical in physical processes to be studied below. Without loss of generality, we always take α c = 0 for simplicity. Thus, the set of input parameters can be summarized as follows Analogous to the CPC version of the general 2HDM, the parameter choice of β − α = π/2 corresponds to the so-called "alignment limit". This can be achieved when taking into account the signal fit to the 125 GeV SM-like Higgs boson h 1 , as shown later. By further combining with the eEDM constraints, we will fix the parameters of t β and α b and constrain two other mass parameters of M and m soft for our later discussions.

The couplings in the CPV 2HDM
For simplicity, we focus on the 2HDMs 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 2HDM-I and 2HDM-II read (and suppressing the CKM mixing), where 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 CP-odd 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.4), which were also previously obtained in Refs. [52][53][54]. Here, we summarize their explicit expressions under the alignment limit in Table. 1. In this alignment limit of β − α = π/2, the Higgs 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. [43,[55][56][57][58][59][60][61][62]. By extracting the cubic terms in the scalar potential Eq. (2.1b), we can obtain the Higgs cubic self-interacting terms. The neutral part of the cubic terms are expressed as follows in the basis of (H 0 From these terms, one can readily obtain the cubic interactions in terms of the mass eigen-states of (h 1 , h 2 , h 3 ) by using the orthogonal mixing matrix R from Eq. (2.4). Throughout our discussions, we define the Higgs cubic self couplings of λ ijk (i , j , k = 1 , 2 , 3) to be the coefficients of the h i h j h k term from Eq. (2.19) where the symmetry factors are such that S! = 3! = 6 for i = j = k, S! = 2 for i = j = k, and S = 1 for i = j = k. A general derivation of the Higgs cubic self couplings in the CPV 2HDM was previously studied in Refs. [63][64][65]. The explicit expressions of λ ijk are tedious, while they can be greatly simplified with the fixed parameters through the following discussions.
3 The Constraints in The CPV 2HDM

The 125 GeV SM-like Higgs boson constraint
In the CPV 2HDM, the productions and decay rates of the 125 GeV SM-like Higgs boson h 1 are controlled by both CP-even couplings of c f ,1 and CP-odd couplings ofc f ,1 . The production cross sections and decay rates are rescaled from the SM one as follows, For the production cross sections and decay rates of the SM Higgs boson, we use the results from the LHC Higgs Working Group given in Refs. [66,67]. The LHC signal strengths of the SM-like Higgs boson in the presence of the CPV were discussed in Refs. [52][53][54][68][69][70][71][72][73][74][75]. From Table. 1, one notes that the relevant Yukawa couplings of (c f ,1 ,c f ,1 ) and the Higgs gauge couplings of a 1 are only controlled by the Higgs VEV ratio of t β as well as the CPV mixing angle of α b . The heavy Higgs bosons in the spectrum are either irrelevant or negligible for the signal fit of h 1 . Based on the most recent LHC measurements of the 125 GeV signal strengths [76][77][78][79][80][81], we fit the signal strength of h 1 on the (t β , |α b |) plane and present the results with the eEDM constraints later.

The eEDM constraints
The ACME experiment [48],which searches for an energy shift of ThO molecules due to an external electric field, set stringent experimental bound to the eEDM. 4 The bound reads The eEDM constraints to the CPV 2HDM-II were previously studied in the Refs. [40,53]. The effective Lagrangian term is given as follows  In the CPV 2HDM, the Wilson coefficient δ e are contributed by the two-loop Barr-Zee type h i γγ(h i Zγ) diagrams [82], and the H ± W ∓ γ diagrams, as depicted in Fig. 1. The h i γγ(h i Zγ) diagrams include the contributions from: (i) the top-quark loops, (ii) the Wboson and the NGB loops, and (iii) the charged Higgs boson loops. The total contributions can be summarized as follows Here, the superscripts of h i γγ, h i Zγ, and H ± W ∓ γ represent the operators for the specific Barr-Zee type diagrams. The subscripts of (t , W , H ± , h i ) represent the particles in the loops. Explicit expression for each term can be found in Refs. [83][84][85], and summarized in the appendix of Ref. [53]. Numerically, the leading contributions to the Wilson coefficient δ e are mainly due to the (δ e ) h 1 γγ and (δ e ) h 1 Zγ terms, while the contributions from the other heavy Higgs bosons of (h 2 ,3 , H ± ) can be safely neglected. These terms are proportional to the CP-odd couplings ofc f ,1 , and further proportional to the CPV mixing angle α b according to the Yukawa couplings listed in Table. 1.
The eEDM upper bound from the ACME is converted to the constraints to the CPV 2HDM parameters on the (t β , |α b |) plane. The combined 125 GeV Higgs boson signal constraints and the eEDM constraints are shown in Fig. 2. It is clear that the eEDM bound is the leading one to set upper bounds to the CPV mixing angle of |α b |, as compared to the fits of the SM-like Higgs boson signal strengths. For the CPV 2HDM-I (left panel), the size of CPV mixing angle is significantly bound as |α b | 5 × 10 −3 , and the 1 σ allowed range of t β is within (0.9 , 1.7). For the CPV 2HDM-II (right panel), the allowed region of the CPV mixing angle can be extended to |α b | 0.1, while the 1 σ allowed range of t β is basically around 1.0. It has been noted in Ref. [40] that the maximal cancellations between the h i F µν V µν operator and the h i F µνṼ µν operator can be achieved with the input of t β ∼ 1 in the CPV 2HDM-II. In order to highlight the CPV effects in the Higgs self couplings in the following discussions, we will focus on the CPV 2HDM-II with the fixed inputs of α = −π/4 and t β = 1.0. Furthermore, we also find that the Higgs cubic self couplings almost approach to the SM limit when the CPV mixing angle can be constrained as small as |α b | 0.01. As stated in the previous paragraph, the Wilson coefficient of δ e depends on the CPV mixings almost linearly. Therefore, if the future measurements of the eEDM can improve the precisions to an order of magnitude or more, they can be very useful to constrain the benchmark models for the Higgs pair productions in this setup.

The unitarity and stability constraints
To have a self-consistent description of the 2HDM potential, two other theoretical constraints should be taken into account, namely, the perturbative unitarity and the stability.
Very roughly speaking, the perturbative unitarity constraint means that the theory cannot be strongly coupled. According to the relations listed in Eqs. (2.11), the constraints to the self couplings of λ i can be converted to upper bounds to the Higgs boson masses and the soft mass term of m soft in the 2HDM. 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 [86,87]. Due to the Nambu-Goldstone theorem, the S-matrices can be expressed in terms of 2HDM quartic couplings λ i . Explicitly, the unitarity conditions to be satisfied are that the eigenvalues of each S-wave amplitude matrix should be ∈ (−1/2 , 1/2). The S-wave amplitude matrices are due to fourteen neutral, eight singly-charged, and three doubly-charged scalar channels. They read The S-wave amplitude matrices for three different channels are expressed as where the expressions for the submatrices of (X 4×4 , Y 4×4 , Z 3×3 ) are given in the Ref. [87]. The stability constraints require a positive 2HDM potential for large values of Higgs fields along all field space directions. Collectively, they lead to the following conditions As seen from Eqs. (2.11), very large values of m soft will pull λ 1 ,2 into the negative regions, which violate the conditions described by Eqs. (3.8). Later, we will find that the Higgs cubic self couplings, such as λ 113 in our case, become enhanced with the large soft mass inputs of m soft when they are close to the stability boundary.

The LHC searches for heavy Higgs bosons
The constraints to the signal strengths of the 125 GeV SM-like Higgs boson h 1 and the eEDM put bounds to the parameters of (|α b | , t β ). The unitarity and stability constraints put upper bounds to the mass input parameters of (M , m soft ). Below, we take into account the constraints from the 7 ⊕ 8 TeV LHC searches for the heavy Higgs bosons in the 2HDM spectrum. Such constraints were previously given in Ref. [54], where authors included the constraints from h 2 ,3 → W W/ZZ and h 2 ,3 → Zh 1 → + − bb final states. Additionally, there have been recent experimental searches to the hh → bbγγ final states from both ATLAS and CMS collaborations, which are included in our studies.

The heavy Higgs productions
The cross sections of the heavy Higgs bosons via the ggF channel can be rescaled from the SM-like Higgs production with the same mass as with the variable of The cross sections of the heavy Higgs bosons via the VBF channel can be rescaled from the SM-like Higgs production with the same mass as

The heavy Higgs decays
Here, we list the partial decay widths of the heavy neutral Higgs bosons at the leading order (LO). The partial decay widths into the gauge bosons are with V = (W ± , Z). The partial decay widths into the SM fermions are We also consider the non-standard decay modes of the heavy Higgs bosons, which include

The experimental search bounds
The current LHC experimental searches for the heavy Higgs bosons are performed via the (W W , ZZ) final states [88,89], the H → hh → bb + γγ [90,91], and A → hZ → (bb + + − /τ + τ − + + − ) [94,95]. Since we always assume that M 2 = M 3 , the constraints to the heavy Higgs boson searches at the LHC are imposed to the cross sections of σ  Table. 2, where the soft mass terms of m soft are chosen to be close to the stability boundary for each heavy Higgs boson mass. In the next section, we will study the Higgs pair productions at the future e + e − and pp collider experiments based on these benchmark models.

The EW precision constraints
The Peskin-Takeuchi parameters of (S , T ) for the EW precision tests were obtained in Refs. [30,[96][97][98][99][100] for the 2HDM. In our simplified case with the alignment limit, the de- for a reference value of the SM Higgs boson mass M H ,ref = 125 GeV. Here, we denote ∆ 1 ≡ M ± − M 1 , and the functions of G(x, y, z),Ĝ(x, y) are given in [96]. By employing the current Gfitter fit to the EW data [101], the parameters are founded to be constrained by T parameter mostly for the CPV parameter α b allowed by Fig. 2, and the degenerate masses of heavy Higgs bosons relax the constraints again.

Higgs Pair Productions at The Colliders
In this section, we study the SM-like Higgs pair productions in the framework of the CPV 2HDM. The SM-like Higgs cubic self coupling of λ 111 are modified due to the varying inputs of the soft mass term and the CPV mixing angle. Therefore, we will discuss the precision measurement of λ 111 at the future e + e − colliders for the benchmark models in Table. 2. We will also focus on the most dominant channel for the resonance contributions, namely the ggF process at the hadron colliders, which include the LHC 14 TeV and the future SppC/Fcc-hh 100 TeV runs.

The Higgs cubic self couplings
Before evaluating the cross sections of the Higgs pair productions, it is necessary to look at the behaviors of the relevant Higgs cubic self couplings of λ 11i (i = 1, 2, 3). Following the previous constraints, we fix the parameters of (α , t β ) = (−π/4 , 1.0), and keep the input of α c = 0. With these assumptions, we find that only λ 111 and λ 113 survive, and λ 112 is always vanishing. Their explicit expressions in the mass eigenbasis read Since the CPV mixing angle of α b is typically small by imposing the eEDM constraints, it is also useful to expand the cubic couplings in terms of the α b angle as follows The Higgs cubic self coupling of λ 111 starts with the SM predicted values of λ SM hhh 32 GeV, plus the higher order corrections of O(α 2 b ). The overall magnitude of λ 113 is controlled by the size of the CPV mixing angle α b . Hence, one can expect that the improvement in the precisions of the eEDM measurements will reduce the size of the heavy resonance contributions to the Higgs pair productions via the ggF process.
In Fig. 5, we plot the Higgs cubic self couplings of λ 111 and λ 113 for the M 2 = M 3 = 600 GeV case with different CPV mixing angles of α b in the CPV 2HDM-II. The lower and upper bounds of the soft mass inputs m soft in these plots are from the perturbative unitarity and the stability constraints, respectively. For a fixed input of α b , the Higgs cubic self coupling of λ 111 becomes smaller than the SM predicted value with the increasing inputs of m soft . On the other hand, when the CPV mixing angle becomes as small as α b = 0.01, the Higgs cubic self coupling of λ 111 is basically the same as λ SM hhh 32 GeV. The other Higgs cubic self coupling of λ 113 also decreases from positive regions to negative regions with the increasing inputs of m soft . Its variation is also controlled by the size of the CPV mixing angle of α b , as seen from its behaviors with the different inputs of the CPV mixing angle of α b = (0.1 , 0.05 , 0.01). For the M 2 = M 3 = 600 GeV case, λ 113 tends to zero when the soft mass term is m soft 220 GeV, as can be evaluated from Eq. (4.2b). Thus, one would expect the corresponding resonance contributions to vanish. When the soft mass deviates from this value of m soft 220 GeV, either increases to the stability boundary or decreases to zero, |λ 113 | increases. Correspondingly, one can expect large resonance contributions for such parameter inputs. The future high-energy e + e − colliders provide opportunities of measuring the SM-like Higgs cubic self couplings. The direct measurements can be achieved via the e + e − → hhZ process with the center-of-mass energy of √ s = 500 GeV, or via the vector boson fusion process of e + e − → hhν eνe with the center-of-mass energy of √ s = 1 TeV [14,15,17]. The first advantage of the e + e − colliders is that the relevant Higgs-gauge couplings for these processes can be precisely measured to the percentage level at the √ s = 240 − 250 GeV runs [14][15][16].
For the CPV 2HDM, one has the Higgs-gauge couplings of with δg 1ZZ = g 1ZZ − g SM hZZ < O(1 %) after imposing the eEDM constraints. The second advantage of the e + e − colliders is that the contributions to the total cross section from the heavy resonance of h 3 are typically less than O(10 −4 ), hence they are negligible. Therefore, it is a good approximation to assume the SM predicted values for the Higgs-gauge couplings, and only vary the Higgs cubic self coupling of λ 111 . The ratio of the total cross section of σ[e + e − → hhZ] to its SM counterpart can be parametrized as follows at the TLEP and ILC 500 GeV runs, with ξ 111 ≡ λ 111 /λ SM hhh . The total cross sections at the TLEP and ILC 500 GeV runs versus the ratios of different Higgs cubic self couplings λ 111 /λ SM hhh are displayed on the left panel of Fig. 6. The ranges of λ 111 in two set of benchmark models with |α b | = 0.1 and |α b | = 0.05 are also shown in the light-blue and light-green shaded regions, respectively. From the results given in Table. 2 for the benchmark models, the Higgs cubic self couplings of λ 111 are always smaller than the SM predicted values. Thus, the corresponding cross sections of σ[e + e − → h 1 h 1 Z] are smaller than the SM predictions at the TLEP and the ILC. On the right panel of Fig. 6, we display the expected accuracies on the Higgs cubic self The parton-level differential cross sections of the Higgs pair production for both SM Higgs and BSM Higgs bosons via the ggF process were previously derived in Refs. [102][103][104][105]. For the productions of the SM-like Higgs boson pairs, its differential cross section reads where the dominant contributions are due to the top-quark loops. The form factors of (F , F 2 , G 2 ) are from the loop integrals of the triangle diagrams, the J = 0 partial wave of the box diagrams, and the J = 2 partial wave of the box diagrams. Their explicit expressions are summarized in the appendix of Ref. [103]. The relevant coefficients are given by For the most general case in the CPV 2HDM, all neutral Higgs bosons of h i have both CP-even and CP-odd Yukawa couplings. Furthermore, the heavy resonances enter into the Higgs pair productions. The corresponding differential cross sections at the parton level can be generalized from the results in the appendix of Ref. [103] for the different CP combinations of the final-state h 1 h 1 , which are expressed as follows (4.9) The relevant couplings are with g 111 = 6λ 111 and g 113 = 4λ 113 . To evaluate the cross sections, we implement all couplings given in Eqs. (4.10) into the FeynRules [107], and pass the UFO model files into the Madgraph 5. Now we present the results of the Higgs pair productions in the CPV 2HDM, by combining all previous constraints. As one can learn from Eq. (4.9), the cross sections of σ[pp → h 1 h 1 ] get modified from their SM counterparts due to: (i) the modification of the Higgs cubic self coupling λ 111 , (ii) the modifications of the top quark Yukawa couplings, and (iii) the additional resonance contributions. Through the signal fit to the 125 GeV SM-like Higgs boson h 1 and the eEDM constraints, the dimensionless Higgs Yukawa couplings are bounded such that δc f ,1 < 1 % andc f ,1 ∼ −0.1. Therefore, the box diagram contributions are envisioned to approach to the SM predicted values. From the previous estimation of the Higgs cubic self couplings for the M 2 = M 3 = 600 GeV case, we may either have the large resonance contributions or go to the regions with the vanishing resonance contributions of (λ 111 , λ 113 ) → (λ SM hhh , 0). For these two limiting scenarios, further simplifications can be made for Eq. (4.9), (4.11b) In Fig. 7, we display the LO cross sections of σ[pp → h 1 h 1 ] at the LHC and the SppC/Fcchh for the M 2 = M 3 = 600 GeV case. The solid curves represent the full results by combining every term in Eq. (4.9). We also show the hypothetical cross sections of Eq. (4.11b), where we turn off the Higgs cubic self coupling of λ 113 while modify λ 111 according to Eq. (4.2a). Thus, it is evident that the total cross sections approach to the SM-like Higgs pair productions with the modified cubic self couplings. On the other hand, the LO cross sections at the LHC (SppC) can be as large as ∼ O(100) fb (∼ O(6) pb) when the soft mass approaches to the stability boundary for this case.
Furthermore, we evaluate the LO cross sections for the benchmark models listed in Table. 2. The typical cross sections subject all constraints in the previous context are O(10)−O(100) fb at the LHC or O(1) pb at the SppC for the allowed mass ranges. The corresponding results are displayed in Fig. 8, for benchmark models with |α b | = 0.1 and |α b | = 0.05, respectively. We display the cross sections with the h 1 h 1 → bb + γγ and h 1 h 1 → bb + W + W − final states. From the experimental side, the bb + γγ final states are the leading one to look for the Higgs pair productions at the hadron colliders, in that the relevant SM background is under control. The LO cross sections for the bb+γγ of our benchmark models are ∼ O(0.1) fb at the LHC, and they increase to O(10) fb at the SppC. In addition, one may also consider the bb + W h W final states with the aid of the jet substructure technique [8]. The LO cross  Besides the Higgs pair productions, we also have the other search modes for the heavy Higgs bosons of h 2 /h 3 , such as di-bosons and Higgs plus Z. In Fig. 9, we display the cross sec-tions of the other search modes of the heavy Higgs bosons, including pp → h i → (W + W − → 2 2ν , ZZ → 4 , hZ → bb + + − ). The current LHC searches for the heavy Higgs bosons via these channels can be found in Refs. [88,89,94,95,[108][109][110]. The cross sections for these benchmark models are typically ∼ O(0.01) − O(0.1) fb at the LHC, and enhanced to O(1) − O(10) fb at the SppC. Analogous to the Higgs pair production process at the resonance region, the decay branching ratios of Br[h i → W W/ZZ/hZ] ∝ α 2 b . Therefore, the improvements of the precise measurements of the future eEDM experiments can also suppress the expected cross sections for these final states.

Conclusion
The extended Higgs sector is a general setup with rich physical ingredients to address the issues that are beyond the SM. Particularly, the spontaneous CPV can be achieved with the general 2HDM setup. In this work, we study the Higgs pair productions in the framework of the CPV 2HDM, with the focus on the leading production channel of the ggF. The set of constraints to the CPV Higgs sector are taken into account, including the SM-like Higgs signal fit, the eEDM constraint, the perturbative unitarity and stability constraints, and the current LHC searches for the heavy Higgs bosons. Together with the simplification to the model, we focus on the CPV 2HDM-II, where a relatively large size of CPV mixing is possible at t β ∼ 1.
The Higgs cubic self couplings play the most crucial role for the Higgs pair production. For our case, two relevant cubic self couplings are λ 111 and λ 113 , which are controlled by the soft mass term m soft and the CPV mixing angle of α b . The precise measurement of the SM-like Higgs cubic coupling of λ 111 can be achieved via the e + e − → h 1 h 1 Z and e + e − → h 1 h 1 ν eνe processes at the future high-energy e + e − colliders. The benchmark models in our discussions typically predict totally cross sections of σ[e + e − → h 1 h 1 Z] smaller than the SM predictions. The largest deviations of the SM-like Higgs cubic couplings λ 111 are likely to be probed at the future TLEP 500 GeV and ILC 1 TeV runs. At the future high-energy pp collider runs, the Higgs pair productions are very likely to be controlled by the heavy resonance contributions. In the allowed mass range of the heavy Higgs bosons, we find the total production cross sections to be σ[pp → h 1 h 1 ] ∼ O(10) − O(100) fb at the LHC 14 TeV runs. They can be as large as ∼ O(10 3 ) fb at the future SppC 100 TeV runs. Other search modes of di-bosons and Higgs plus Z that are currently probed at the LHC 7 ⊕ 8 TeV experiments are also estimated at the future LHC 14 TeV and the SppC 100 TeV experiments. The discovery of all these channels will manifest the structure of the Higgs sector. Therefore, it will be very helpful to further study the higher-order QCD corrections as well as the collider search capabilities for such heavy resonance contributions to the Higgs pairs.