Future prospects of mass-degenerate Higgs bosons in the $CP$-conserving two-Higgs-doublet model

The scenario of two mass-degenerate Higgs bosons within the general two-Higgs-doublet model (2HDM) is revisited. We focus on the global picture when two $CP$-even Higgs bosons of $h$ and $H$ are nearly mass-degenerate. A global fit to the signal strength of the 125 GeV Higgs measured at the LHC is performed. Based on the best-fit result of the 2HDM mixing angles $(\alpha,\beta)$, theoretical constraints, charged and $CP$-odd Higgs boson direct search constraints and the electroweak precision constraints are imposed to the 2HDM parameter space. We present the signal predictions of the $(4b\,, 2b\,2\gamma)$ channels for the benchmark models at the LHC 14 TeV runs. We also study the direct Higgs boson pair productions at the LHC, and the Z-associated Higgs boson pair production search at the ILC 500 GeV runs, as well as the indirect probes at the CEPC 250 GeV run. We find that the mass-degenerate Higgs boson scenario in the Type-II 2HDM can be fully probed by these future experimental searches.


I. INTRODUCTION
When the 125 GeV Higgs boson was discovered at the LHC 7 ⊕ 8 TeV runs [1,2], the experimental measurements of the γγ signal rates from both ATLAS and CMS collaborations were both enhanced relative to the standard model (SM) predictions. It was suggested in Refs. [3][4][5][6] that the observed signals at ∼ 125 GeV may arise from two mass-degenerate Higgs bosons. 1 Even though the obvious enhancement of the γγ rate is not shown at the LHC run-II, the scenario of two mass-degenerate Higgs bosons around 125 GeV still deserves investigation. After all, it is really challenging to distinguish this possibility from the single Higgs boson case by direct measurements of the Higgs boson mass, given that the energy resolutions of photons and leptons are typically of ∼ O(1) GeV at the LHC [10][11][12][13][14]. The current mass uncertainties from the CMS measurements are ∼ 0.24 GeV [15,16]. The direct measurements of the Higgs boson(s) at 125 GeV involve their gauge couplings and Yukawa couplings at the leading order (LO). Alternatively, one may constrain such a scenario from a global point of view, by imposing various theoretical and experimental constraints.
In this work, we study the future experimental prospects of probing the mass-degenerate 125 GeV Higgs bosons at both high-luminosity (HL) LHC runs and the future high-energy colliders. Our discussions are made in the context of the CP -conserving (CPC) general 2HDM. This scenario can be constrained from the current LHC searches for the CP -odd Higgs boson A via the h/H +Z decay channel. Furthermore, we suggest to distinguish the h/H mass-degenerate case from the single resonance case through the probes of the Higgs boson self couplings. This can be done by searching for the Higgs boson pair production processes at both LHC and the future high-energy colliders. For the new physics (NP) models involving a single 125 GeV Higgs boson, it is quite often that the modified Higgs cubic self couplings (including the additional resonances) are the only sources to modify the Higgs pair production cross sections. The signal rates for various final states can be estimated by using the SM-like Higgs boson decay branching fractions (see Ref. [17] for the summary). There have been extensive discussions of the Higgs boson pair productions in various beyond standard model (BSM) NP models [18][19][20][21][22][23]. Currently, two most sensitive search modes for the Higgs boson pair productions at the LHC 13 TeV run are (4b , 2b 2γ) [24][25][26][27][28][29]. For the Higgs pair productions with two mass-degenerate Higgs bosons, one may expect: (i) the deviation of Higgs cubic self-couplings from the SM predictions, and (ii) the existence of multiple Higgs cubic self-couplings, hence, multiple processes contributing to each final state.
The layout of this paper is described as follows. In Sec. II, we review the scenario of degenerate Higgs bosons in the framework of CPC 2HDM. The LHC measurements of the Higgs signal strengths are used for the global fit, where we simplify the discussion with negligible quantum interference and mixing effects. This can be achieved by assuming sufficiently large mass splitting. Other constraints, such as the perturbative unitarity and stability of the 2HDM potential, the EW precision tests, as well as the LHC direct searches for the CP -odd Higgs boson A, are also considered for the 2HDM with mass-degenerate h/H. The benchmark points are suggested for both Type-I and Type-II 2HDM. In Sec. III, we study the gluon-gluon fusion (ggF) productions of Higgs pairs at the LHC for the degenerate Higgs scenario. We compare the signal predictions from various final states of (4b , 2b 2γ) with the corresponding SM predictions at the LO. Their cross sections are generally varying with different soft mass terms of m 12 in the 2HDM potential. In particular, we find that the signal rates of (4b , 2b 2γ) final states are always moderately enhanced with respect to the SM predictions. The corresponding significances are estimated for the h/H mass-degenerate case as well. In Sec. IV, we discuss the capability of distinguishing the h/H mass-degenerate scenario at the future high-energy e + e − colliders. We show the indication from the precise measurement of cross sections of mass-degenerate Higgs bosons with Z-boson for this scenario at the circular electronpositron collider (CEPC). Furthermore, the direct production of Higgs boson pairs associated with Z-boson at the ILC can probe the h/H mass-degenerate scenario in the Type-II 2HDM. The summaries are given in Sec. V.

II. THE MASS-DEGENERATE HIGGS BOSONS IN THE 2HDM
A. The global fit to the mass-degenerate Higgs boson signals at the LHC In the CPC 2HDM, there are five Higgs bosons of (h , H , A , H ± ) in the scalar mass spectrum. The review of the 2HDM setup and the related LHC phenomenology can be found in Refs. [7,30]. The Lagrangian for the general 2HDM is written as follows All parameters are assumed to be real for the CPC case. Very often, a softly broken Z 2 symmetry, under which two Higgs doublets transform as ( , is also assumed to eliminate the possible λ 6 ,7 couplings in the 2HDM potential. One has m 12 = 0 when the Z 2 symmetry is exact. Two Higgs doublets of Φ 1, 2 can be expressed in terms of components as with two Higgs vacuum expectation values (VEVs) and their ratios being Here, π 0 1 ,2 are pseudoreal components, whose linear combinations of A = −s β π 0 1 + c β π 0 2 and G = c β π 0 1 + s β π 0 2 are CP -odd Higgs boson and neutral Nambu-Goldstone boson, respectively. π + 1 ,2 (and their complex conjugates) are complex scalar fields, whose linear combinations of H ± = −s β π ± 1 + c β π ± 2 and G ± = c β π ± 1 + s β π ± 2 are charged Higgs bosons and charged Nambu-Goldstone bosons, respectively.
For further discussion, we list the dimensionless Higgs gauge couplings and Yukawa couplings as follows with Type-II : Here, α represents the mixing angle between two CP -even Higgs bosons of (h , H). The overall signal rates are controlled by the input parameters of (α , β) in the CPC 2HDM, which is manifest from the couplings in Eqs. (5). A global fit to the h/H degenerate scenario respect to (α , β) is thus performed, and this is done by the χ 2 fit to the LHC data defined as where the current LHC measurements of signal strengths and errors of (µ PD exp , σ PD exp ) are summarized in Tables I and II for the run-I and run-II data, respectively. This scenario was previously explored in Refs. [3][4][5][6][7], where the total signal rates for 125 GeV Higgs boson were estimated by simple summation of σ × Br from the individual contribution of h and H as  [41] · · · · · · τ + τ − VBF · · · · · · 0.94 ± 0.41 [42] τ + τ − VH · · · · · · -0.33 ± 1.02 [42]  simplify our discussions, we still use the simple summation method in Eq. (7) to estimate the total signal rates for various channels. The global fit results on the (α , β) plane are shown in Fig. 1. The corresponding benchmark points for the mass-degenerate M h ≈ M H = 125 GeV cases are listed in Table. III for both Type-I and Type-II 2HDM. We also observe that a shift of α → α − π/2 leads to equally minimal χ 2 values in both Type-I and Type-II cases. This corresponds to an interchange between h and H in the mass-degenerate scenario. Besides the 2HDM parameters of (α , β), the total decay widths of Γ h + Γ H and the main decay branching fractions are also listed for the massdegenerate case. Since the total decay widths of Γ h + Γ H are smaller than ∼ O(0.1) GeV, our simplification in Eq. (7) is valid. The alignment parameters are c β−α = 0.21 for the Type-I 2HDM and c β−α = 0.01 for the Type-II 2HDM, respectively. A sizable deviation from the alignment limit is observed in the Type-I benchmark point. For the Type-II case, meanwhile, the H is gaugephobic and decays mostly into fermionic final states of (bb , τ + τ − ). Throughout the context below, we shall always use the best-fit points of (α , β) in Table. III for the phenomenology discussions in ZZ VH 0 ± 1.9 [44,46] 0 ± 2.83 , or 0 ± 2.66 [47] ZZ ttH 0 ± 3.9 [44,46] 0 ± 1.19 [47] W + W − ggF · · · · · · 1.02 ± 0.27 [48] · · · · · · W + W − VBF+VH · · · · · · 0.89 ± 0.67 [48] bb VH 1.20 +0.42 −0.36 [50] · · · · · · τ + τ − ggF · · · · · · 0.84 ± 0.89 [51] τ + τ − VBF · · · · · · 1.11 +0.34 −0.35 [51] τ + τ − ttH · · · · · · 0.72 +0.62 −0.53 [52]  the mass-degenerate scenario. Besides the best-fit points from the current LHC Higgs data, we shall further impose theoretical and experimental constraints to the 2HDM mass spectrum of (M A , M ± , m 12 ) in the following context, for the mass-degenerate Higgs boson scenario. We shall show that mass-degenerate Higgs boson scenario has allowed 2HDM parameter space with all these constraints imposed.

B. The charged Higgs boson and EW precision constraints to the 2HDM
It is known that the charged Higgs bosons of H ± contribute to the flavor-changing neutral current (FCNC) rare decay processes, such as b → sγ transition. The latest measurement is from the Belle Collaboration [53], and the implication to the CPC 2HDM was carried out in Refs. [54][55][56][57]. By imposing the FCNC constraints to the benchmark models in Table. III, we get M ± 590 GeV in the Type-II 2HDM, while the lower mass bound in the Type-I 2HDM is negligible, as compared to the direct collider constraints. Besides, the direct searches for the charged Higgs bosons at the LHC were performed in Refs. [58][59][60]. Here, we shall only consider the FCNC constraints to the   shall take M ± = M A for the Type-I 2HDM 2 , and fix M ± = 600 GeV for the Type-II 2HDM below 3 . We consider the constraints from the EW precision tests [61][62][63] to the 2HDM with massdegenerate h/H. The most general expressions for (∆S , ∆T ) in the CPC 2HDM [61] read where we explicitly split these expressions into terms independent of or dependent on the alignment parameter of c β−α . The relevant auxiliary functions read The current Gfitter fit [64] to the EW data gives S = 0.05 ± 0.11 , T = 0.09 ± 0.13 .
2 As we shall see below, the specific mass ranges of (M± , MA) do not play a role in the Higgs boson pair productions at the LHC or ILC. Without loss of generality, we make such simplification of M± = MA. 3 When taking the unitarity bound into account, it turns out that the charged Higgs boson mass cannot exceed ∼ 625 GeV. Thus, a fixed M± = 600 GeV is taken to compromise the joint constraints from the FCNC rare decay and unitarity for the Type-II case.
For benchmark models in both Type-I and Type-II cases, the alignment parameters were found to be small as from Table. III. Thus, the 2HDM contributions to the (∆S , ∆T ) are mainly controlled by leading terms in the first lines of Eqs. (8). By using the definitions of auxiliary functions of (9a) and (9e), the (∆S , ∆T ) can be suppressed with degenerate mass inputs of M A = M ± . Indeed, by using the best-fit (α , β) inputs for the h/H degenerate cases in Table. III,  The joint constraints of the perturbative unitarity and tree-level stability conditions to the 2HDM potential turns out to be powerful to bound the heavy scalar masses. The conditions to be satisfied for the unitarity constraints to the 2HDM potential are that the absolute values of the following linear combinations of the quartic scalar couplings [65,66]: should be smaller than or equal to 8π. The tree-level vacuum stability conditions for the general 2HDM potential come from the requirement that the scalar potential being bounds from below, which read [67] 5 The quartic self couplings of λ i are related to the Higgs boson masses, the mixing angles, and the soft mass term as follows 4 Here, we take the mass range of MA ∈ (200 , 300) GeV by considering the perturbative unitarity constraint and the direct search limit of A → hZ at the LHC. 5 Recently, it was suggested in Ref. [68] to apply the global minimum condition to constrain the 2HDM potential.
In Ref. [69], the loop effects to the vacuum stability conditions were found to alleviate the tree-level conditions.
By combining the constraints in Eq. (11) and Eq. (12) and using the quartic self couplings given in Eqs. (13), we show the joint unitarity and stability constraints in Fig. 2 for Type-I and Type-II 2HDM. The best-fit points of (α , β) for Type-I and Type-II cases are used as in Table. III. As mentioned in the previous subsection, a fixed charged Higgs boson mass of M ± = 600 GeV is always taken in the Type-II case to evade the B-physics constraints. For the best-fit points of (α, β) in Table. III, a larger m 12 leads to a larger negative scalar quartic couplings λ 1 as indicated by Eq. (13a), therefore results in the perturbative unitarity (mostly from the |a − | ≤8π) and stability bounds on the m 12 , as depicted by the two panels of Fig. 2. For the fixed m 12 , one can expect a larger positive λ 3 for larger heavy Higgs boson masses. This breaks the perturbative unitarity through quartic coupling combination of e 2 , therefore sets the upper bounds on the heavy Higgs boson masses for both Type-I and Type-II 2HDM. For the fixed heavy Higgs boson masses of M A = M ± in the Type-I case, a smaller m 12 leads to a relatively larger positive λ 1 . This results in a larger |a + |, which in turn gives the lower bounds on m 12 on the left panel of Fig. 2. The upper bounds to the mass mixing of m 12 in the Higgs potential set by the unitarity constraints and the stability constraints, which mainly come from the fact that the mass of the second CP -even Higgs boson M H is fixed. Correspondingly, the quartic scalar couplings of λ 1 ,2 are determined by m 12 for the best-fit points. Since the m 12 enters into the Higgs cubic self couplings, we take their ranges to be Here, we also limit the heavy Higgs boson masses of (M A , M ± ) in the ranges that are consistent with the current LHC run-II searches for the CP -odd Higgs boson A in the mass-degenerate scenario, as indicated in Fig. 4 below.
The coefficients of the Higgs cubic self couplings in the physical basis will be used in the calculation of the direct Higgs pair productions at the LHC and the ILC, which are listed as follows In the alignment limit of β − α = π/2, they become In Fig. 3, we display the Higgs cubic self couplings versus the soft mass term m 12 , with the unitarity/stability constraints in Eq. (14) taken into account for two best-fit points listed in Table. III. It turns out that the Higgs cubic self couplings of λ hhh are very close to the SM value of λ SM    Table. III. The ranges of m 12 are taken as in Eq. (14) for Type-I and Type-II cases.

D. The constraints from the CP -odd Higgs boson A searches
Before we discuss the degenerate Higgs boson pair searches at the LHC, we impose the constraints via the LHC searches for the CP -odd Higgs boson A in the mass-degenerate h/H scenario. The decay modes and the corresponding partial decay widths of CP -odd Higgs boson A are with N c,f = 3 (1) for quarks (leptons). The three-body phase space factor reads For the mass-degenerate h/H scenario, one cannot discriminate two decay channels of A → hZ and A → HZ for specific final states, such as bb + + − . The signal rates for this scenario should be evaluated as  Table. III.
The LHC searches for a CP -odd Higgs boson via the A → hZ → bb + + − were previous carried out by both ATLAS [71] and CMS [72] collaborations at the 8 TeV run. The most recent results from the ATLAS searches at the LHC 13 TeV run with integrated luminosity of 36.1 fb −1 is given in Ref. [70]. We estimate the production cross sections of σ(gg → A) at the NLO by using the package of SusHi [73], by using the best-fit points as in Table. III. From Fig. 4, we find that the h/H mass-degenerate scenario for either Type-I or Type-II has been tightly constrained by the current exclusion limits from the LHC 13 TeV run with an integrated luminosity of 36.1 fb −1 . Combining with the joint perturbative unitarity and stability bounds shown in Fig. 2, the h/H mass-degenerate scenario can still exist in the mass ranges of M A 280 GeV for Type-I 2HDM and M A 260 GeV for Type-II 2HDM, respectively. Previously, we further restrict the mass ranges of M A in Eq. (14), in order to maximize our parameter choices of m 12 for the Higgs pair productions. It does not mean the mass ranges of M A in Eq. (14) are constrained by the current LHC search result. The specific mass of the CP -odd Higgs boson A will not play any role in our discussion of the future experimental tests below. On the other hand, the future searches for the CP -odd Higgs boson A via this channel at the LHC 14 TeV run will play a decisive role in justifying or falsifying this scenario.

III. THE LHC SEARCHES FOR DEGENERATE HIGGS BOSONS: CONSTRAINTS AND PAIR PRODUCTIONS
A. The total cross section of gluon-gluon fusion to Higgs pairs By using the best-fit points in Table. III and the range of m 12 in Eq. (14) after the set of constraints, we are ready to present the main results of the Higgs pair productions at the LHC. The cross sections we need to evaluate are where the individual cross section of σ[gg → h i h j ] was first obtained in Ref. [74] for both SM and MSSM cases. The differential cross section at the LO reads with c ij = 1/2 (1) for i = j (i = j). C and C 2 represent the coefficients of the triangle and box diagrams, respectively. The Higgs cubic self couplings contribute to the C 's, and they read The coefficients of C 2 are determined by the dimensionless Yukawa couplings of the Higgs bosons for q = (t , b). The form factors of the triangle and box diagrams (F , F 2 , G 2 ) in Eq. (21) can be found in the Appendix of Ref. [74]. The asymptotic behaviors of these form factors in the large quark mass and small quark mass limits read In practice, we evaluate the corresponding Passarino-Veltman (PV) integrals are evaluated by using the LoopTools package [75]. At the LHC, the differential cross section in the lab frame is obtained by convoluting the partonlevel cross section in Eq. (21) with the gluon PDFs and σ[gg → HH] are sub-leading ones, yet they play a role in determining the signal rates. A dip is shown in the cross section of σ[gg → hH] versus m 12 , which roughly matches the position where the Higgs cubic self coupling λ hHH flips sign. The next-to-leading order contributions to the Higgs pair productions at the LHC are known to be significant [77][78][79][80][81]. Our estimation below focus on the future experimental significances via various channels, where a same K-factor can be assumed for both h/H mass-degenerate case and the SM Higgs boson case, as in Ref. [82]. Therefore, the LO results are sufficient for our estimation below.

B. The cross sections of Higgs pairs to various final states
Next, we proceed to present the cross sections of the Higgs pairs to specific final states with the mass-degenerate Higgs benchmark points in Table III. Two leading final states of 4b and 2b 2γ will be taken into account [83]. The signal rates are estimated as follows with κ XY = 1 (2) for X = Y (X = Y ). For a single SM-like Higgs boson with mass ∼ 125 GeV, the ratio of signal rates between the 4b and 2b 2γ final states is fixed to be σ[4b] : σ[2b 2γ] ≈ 127 : 1. This always holds, no matter how one modifies the SM-like Higgs cubic self couplings and includes the additional resonance contributions. We find the ratios of signal rates between these two channels are generally different from the single SM-like Higgs boson case, which reads σ   Table. III. In comparison, the LO cross sections of σ[gg → h SM h SM ] to the corresponding final states (dotted lines) are displayed as well. Note that the allowed regions obtained from theoretical constraints in m 12 are also indicated by the gray and pink arrows for Type-I and Type-II respectively.
The LHC 14 TeV cross sections for σ[gg → (hh , hH , HH) → (4b , 2b2γ)] final states versus the soft mass term m 12 are shown in Fig. 6. We find enhancements of both 4b and 2b 2γ signal rates, in comparison to the SM Higgs boson pair productions. This kind of enhancements were previously investigated in Ref. [84]. However, the parameter region of m 12 is severely restricted by the theoretical constraints as shown in Fig. 2. The LO cross sections for the 4b final states are moderately enhanced to ∼ 7.7 fb with m 12 = 50 GeV in the Type-I case, or ∼ 7.6 fb with m 12 = 60 GeV in the Type-II case. The corresponding LO cross section for the 4b channel in the SM is 6.8 fb. By extrapolating the current LHC run-II results from 13 TeV, the ATLAS and CMS give conservative estimations to the significances for the SM Higgs boson pairs at the HL LHC runs as [83] ATLAS : 1.05 σ for 2b 2γ CMS : (0.39 σ , 1.6 σ) for (4b , 2b 2γ) We summarize the significance of the Type-I and Type-II h/H mass-degenerate Higgs boson pairs via the 4b and 2b 2γ channels in Table. IV. We note that the current LHC run-II results are not optimal for the HL-LHC runs, future improvements to the significance via the 2b 2γ channel should be expected. The future plans of the high-energy e + e − colliders include the CEPC [85], ILC [86], and TLEP [87]. They will play a role as Higgs factory to produce millions of SM-like Higgs bosons for the precise measurements, with the running at center-of-mass energy of √ s = 250 GeV, which will provide excellent opportunity for us to examine the Higgs properties in many NP models [88]. It was pointed out in Ref. [89] that the precise measurement of the Higgs production cross section at the CEPC or ILC is sensitive to the Higgs self couplings at the NLO. At the ILC, it is likely to upgrade the center-of-mass energy up to √ s = 500 GeV, so that it can directly produce 125 GeV Higgs boson pairs associated with a Z-boson.

A. The CEPC measurements of the degenerate Higgs scenario
The circular electron-positron collider (CEPC) will operate at the center-of-mass energy of √ s = 250 GeV. A key physical goal of CEPC is to measure the Higgs boson mass precisely, which can be as small as ∼ 5.9 MeV with an integrated luminosity of 5 ab −1 . The precision on σ(ZH) is about 0.51% combining all decay modes for the Z-boson with the same luminosity [85]. The resolution for the recoil mass peak measurement is about 400 MeV for ILC [90][91][92].
With the integrated luminosity of 5 ab −1 , CEPC can measure the cross section of σ[e + e − → hZ] with an accuracy of ∼ 0.51 % [85], by combining both the leptonic and hadronic decay modes of Z-bosons. The LO cross section of the SM Higgs boson production associated with Z-boson at the CEPC reads σ[e + e − → h SM Z] 221.54 fb 6 . The cross sections for the best-fit points in Table. [85], there will be roughly 7 σ deviation. Therefore, a decrease of the cross section for the σ[e + e − → hZ] will be a first indication from the best-fit points in Table.   The ILC can directly produce Higgs boson pairs associated with a Z-boson, when it runs at the 6 See Refs. [93][94][95][96][97] for more precisely prediction including higher order corrections.
center-of-mass energy of √ s = 500 GeV. The Feynman diagrams for the corresponding processes are depicted in Fig. 7. The cross sections at the ILC can be expressed as [98] σ = 1 4 (1 + P e − )(1 + P e + )σ RR + (1 − P e − )(1 − P e + )σ LL where σ LR denotes the cross section at beam polarization configurations of (P e + , P e − ) = (+1 , −1), and etc. The ILC will run at √ s = 500 GeV with an integrated luminosity of 4 ab −1 , which will be equally shared by two beam polarization configurations of (P e + , P e − ) = (±0.3 , ∓0.8) [86,99] Prospects of measuring the cross sections of the SM Higgs boson pair production via the hh → bbbb and hh → bbW + W − final states have been investigated in Refs. [100,101]. The corresponding significance are listed in Table. V. Combining these measurements for (P e + , P e − ) = (±0.3 , ∓0.8), we obtain the ILC precision on the measurement of (σ RL , σ LR ) in Fig. 8, with two contours for 1and 2-σ regions, respectively. In Fig. 8, we also show the best-fit points for Type-I (Black) and Type-II (Red) 2HDM in the same plane. Each point represents different value of m 12 given by the corresponding label. Note that the 2HDM cross sections are calculated by using the best-fit points in the (α , β) plane from the LHC global fitting in previous section. From Fig. 8, we find that for the h/H mass-degenerate Type-I 2HDM case, the cross sections for all allowed m 12 are still within the precision of the ILC measurement. However, the h/H mass-degenerate Type-II 2HDM points can be excluded from the Higgs pair production measurements at the ILC. More details can be found in Fig. 9, where we combine hh → bbbb and hh → bbW + W − channels and use the log-likelihood ratio method to perform the hypothesis test against the SM predictions. panel for Type-I and Type-II cases, respectively. Dashed vertical lines with arrows indicate the allowed region from 2HDM potential unitarity/stability constraints. We find that for Type-II case, the theoretical allowed region of m 12 has already been excluded by the ILC measurement at 4 σ or more, with 2 ab −1 luminosity for each beam polarization configuration. In contrast, the allowed region for the Type-I case is still safely sitting within 1 σ region of the ILC measurement. In the right panel, we present the exclusion level as a function of the luminosity cumulated at the ILC for beam polarization configurations of P (e − , e + ) = (±0.8, ∓0.3). The black and red line are for Type-I and Type-II case, respectively. For illustration, we choose the most sensitive value of m 12 in the stability allowed region for each type: m 12 = 50 (0) GeV for Type-I (Type-II). We can see that, the ILC has sensitivity to Type-II in this case when an integrated luminosity of 400 fb −1 (equally shared by two beam polarization configurations) is cumulated. However, at least an integrated luminosity of 80 ab −1 is required for the ILC to be sensitive to the Type-I in this situation, which is unrealistic.

V. CONCLUSION
In this work, we study the future prospects of distinguishing the mass-degenerate Higgs boson scenario from the single resonance case at the LHC and future e + e − colliders. Our study is performed in the general CPC 2HDM framework, with two CP -even Higgs bosons of h/H to be mass-degenerate. The direct measurements of the Higgs boson signal rates at the 125 GeV only probe its/their gauge couplings and Yukawa couplings. Alternatively, we find this scenario can be further constrained by a series of theoretical bounds and direct experimental searches in such a framework. Moreover, we suggest that the study of the Higgs boson pair productions will be useful for this scenario. Specifically, there are four types of Higgs cubic self couplings involved in the Higgs pair productions, which are (λ hhh , λ hhH , λ hHH , λ HHH ). The physical processes to be considered are the ggF to Higgs boson pair productions at the LHC, and the Higgs boson pair productions associated with a Z-boson at the ILC.
By performing the global fit of the LHC measurements of the 125 GeV Higgs boson signal rates, we find the best-fit points with mass-degenerate h/H in Type-I and Type-II 2HDM in (α , β) plane. The best-fit point in the Type-I case deviate from the alignment limit as large as c β−α 0.21, while it approaches to the alignment limit as c β−α ∼ O(10 −2 ) in the Type-II case. Correspondingly, one Higgs boson H becomes almost gaugephobic in the Type-II case. The h/H mass-degenerate scenario also passes the current LHC run-II searches for the CP -odd Higgs boson A via the bb + + − final states. Meanwhile, it also means that the further LHC searches for the CP -odd Higgs boson A below the tt mass threshold will play a role to justify or falsify the h/H mass-degenerate scenario.
The relevant cubic Higgs self couplings, such as λ hHH and λ HHH , are not vanishing even when the 2HDM parameters approach to the alignment limit. This suggests that the Higgs pair production processes can crucial to justify or falsify the h/H mass-degenerate scenario. The signal predictions of the ggF to mass-degenerate Higgs boson pairs are made at the LHC 14 TeV runs, with the focus on two leading search channels of 4b and 2b 2γ. Moderate signal enhancements with respect to the SM predictions are expected, while the enhancements are at most ∼ 10 %. Therefore, the Higgs boson pair productions at the LHC are less likely to probe the h/H mass-degenerate scenario. At the ILC 500 GeV run, we find that the h/H mass-degenerate samples in the Type-I case are within the precision of the ILC measurement, while the h/H mass-degenerate samples in the Type-II case can be probed or excluded with an integrated luminosity of 400 fb −1 . It means that the ILC 500 GeV run offers an opportunity to fully probe the h/H mass-degenerate scenario in the Type-II case.
Though our predictions are model-dependent, the Higgs pair productions with two massdegenerate Higgs bosons can be generalized to any other NP models with multiple Higgs bosons. There are multiple Higgs self couplings involved in the Higgs pair productions in general. Depending on the model setup, these self-couplings may be bounded by the constraints mentioned in the current study. Our discussion through the context focus on the ggF process at the LHC and the Higgs pair strahlung at the ILC. This discussion can be extended to other Higgs pair production channels, including the vector boson fusion (VBF) and tt associated processes at the future high-energy e + e − and pp colliders.