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 the CP-violating mixing angles are related to the Higgs cubic self couplings in this setup. Based on these constraints, we suggest benchmark models for the future high-energy collider searches for the Higgs pair productions. The e+e− colliders operating at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = (500 GeV, 1 TeV) are capable of measuring the Higgs cubic self couplings of the benchmark models directly. Afterwards, we estimate the cross sections of the resonance contributions to the Higgs pair productions for the benchmark models at the future LHC and SppC/Fcc-hh runs. Other possible decay modes for the heavy Higgs bosons are also discussed.


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 of bbγγ. Some of the detailed studies at the LHC can be found in refs. [3][4][5][6][7][8][9][10][11][12][13]. From the experimental side, it is well-known that several future high-energy collider programs, such as the International Linear Collider (ILC) [14] in Japan, the Future eplus-eminus/hadronhadron Cicular Collider (Fcc-ee/Fcc-hh) [15] at CERN, and the Circular electron-positron -1 -

JHEP09(2016)069
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 section 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 section 5.

The CPV 2HDM potential
In the general 2HDM, two Higgs doublets of (Φ 1 , Φ 2 ) ∈ 2 +1 are introduced in the scalar sector. For simplicity, we consider the soft breaking of a discrete Z 2 symmetry, under which two Higgs doublets transform as (Φ 1 , Φ 2 ) → (−Φ 1 , Φ 2 ). The corresponding Lagrangian is expressed as 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

JHEP09(2016)069
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 JHEP09(2016)069 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 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.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 of c 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 eigenstates 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.

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 figure 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 JHEP09(2016)069 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 figure 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  singly-charged a + 0 : 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 The stability constraints of (3.8) bound the soft mass term of m soft from above. 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  where

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 σ[pp → where we consider the leading production channel of ggF obtained from eq. (3.9). The decay branching ratios are obtained from the partial decay widths of eqs. (3.12), (3.13), (3.14a), and (3.14b) evaluated at the LO. We find the most stringent constraint to the heavy Higgs boson searches are from the recent CMS searches for the resonances with two SM-like Higgs bosons in ref. [91]. By converting all heavy Higgs boson constraints to the (M , m soft ) plane, we find the mass regions of M 2 ,3 600 GeV are excluded for |α b | = 0.1, or M 2 ,3 500 GeV are excluded for |α b | = 0.05, respectively. The current B-physics data also excludes the charged Higgs boson mass greater than M ± ∼ 340 GeV for 2HDM-II [92,93]. Combining with the previous unitarity and stability constraints, we display the allowed parameter regions of (M , m soft ) in figure 4. Accordingly, we consider two scenarios of  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 degenerate masses of M 2 = M 3 = M ± , and α c = 0, they read 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 figure 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 -14 -

JHEP09(2016)069
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 figure 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.

4.2
The precise measurement of λ 111 at the future e + e − colliders 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 figure 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 lightgreen 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  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 , G ) 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 C hh = (c q ,1 ) 2 , (4.10a) C AA = (c q ,1 ) 2 , (4.10b) C hA = c q ,1cq ,1 , (4.10c) 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), which are resonances :   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 sections for the bb + W W of our benchmark models are ∼ O(10) fb at the LHC, and they increase to O(1) pb at the SppC.

Other channels
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 figure 9, we display the -20 -cross sections 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.