Double Higgs boson production at $e^+e^-$ colliders in the two-Higgs-doublet model

We study the double Higgs boson production processes $e^+e^- \to hh f\bar{f}$ ($f\neq t$) with $h$ being the 125 GeV Higgs boson in the two-Higgs-doublet model with a softly-broken $Z_2$ symmetry. The cross section can be significantly enhanced, typically a few hundreds percent, as compared to the standard model prediction due to resonant effects of heavy neutral Higgs bosons, which becomes important in the case without the alignment limit. We find a strong correlation between the enhancement factor of the cross section and the scaling factor of the $hf\bar{f}$ couplings under constraints from perturbative unitarity, vacuum stability and current experimental data at the LHC as well as the electroweak precision data.


I. INTRODUCTION
Various signatures of the discovered Higgs boson at the LHC [1,2] show that its properties such as observed production cross sections and decay branching ratios are consistent with those predicted in the Standard Model (SM) [3]. Although the SM assumes the minimal form of the Higgs sector composed of one isospin scalar doublet, one may ask him-or herself a natural question: whether the discovered Higgs boson comes from just one doublet or not?
In fact, the discovered Higgs boson can be regarded as one of Higgs bosons arising from an extended structure of the Higgs sector, and there is no strong reason to restrict the Higgs sector to be minimal. On the other hand, extended Higgs sectors often appear as a low energy effective theory of physics beyond the SM based upon various physics motivations.
The important thing is that phenomenological features in the extended Higgs sectors strongly depend on a specific scenario of the underlying theory. Therefore, as a bottom-up approach, it is important to study the structure of the Higgs sector in order to narrow down new physics models.
Among various possibilities of the extended Higgs sector, a two-Higgs-doublet model (THDM) [4] is one of the simplest but important examples, as it appears in several new physics models. For example, models proposed to solve the gauge hierarchy problem predict THDMs as their low energy effective theories such as the minimal supersymmetric SM [5] and composite Higgs models [6][7][8]. In addition, extra CP violating phases can arise from the multi-doublet structure of the scalar potential, which are needed to realize the successful scenario for the electroweak baryogenesis [9,10]. Furthermore, the second Higgs doublet is often introduced in models beyond the SM to explain tiny neutrino masses [11][12][13] and dark matter in the Universe [14]. For these reasons, we consider the THDM as a reference model, in which we impose a softly-broken discrete Z 2 symmetry to avoid flavor changing neutral currents at tree level [15].
There are basically two ways to test the THDM at collider experiments, namely, the direct searches for additional Higgs bosons such as a charged Higgs boson and the indirect searches finding deviations in properties of the discovered Higgs boson (h) from the SM prediction.
Concerning the former way, if we discover the additional Higgs bosons, it turns out to be a direct evidence for the THDM, but no report has hitherto been provided for the discovery of such new particles at the LHC. In contrast, the latter way is getting much attention, since h has already been discovered and its properties will be precisely measured in near future. For example, at the high-luminosity LHC, the h couplings to weak gauge bosons (hV V , V = W, Z) and fermions (e.g., hτ τ , hbb) are expected to be measured with a few percent level [16], while at the International Linear Collider (ILC), hV V and hff (f = τ, c and b) couplings can be measured to be sub-percent and one percent level, respectively [17].
On the other hand, in the THDM the h couplings can deviate from the SM prediction with various patterns depending on a particular scenario as it has been clarified in Ref. [18] at tree level and in Refs. [19][20][21][22][23] at one loop level. Therefore, by looking at the deviation in the h couplings, which could be found in the future collider experiments, we can determine the scenario of the THDM.
In this paper, we focus on the cross section of double Higgs boson production processes at future e + e − colliders as another important observable regarding the indirect search of the THDM. The double Higgs boson production has been discussed in the pp [24], e + e − [24,25] and γγ [24] collision in the THDM to extract the Higgs boson self-coupling constant hhh, particularly the case with the alignment limit, where all the SM-like Higgs boson coupling become the same as those in the SM prediction. In Ref. [26], the Higgs boson pair production has also been discussed at the LHC in the case with and without the alignment limit. Our motivation to discuss the double Higgs production is to find the correlation between the deviation in the Higgs boson couplings and the modification of the cross section. These two variables are expected to be strongly correlated with each other, because the deviation in the h coupling appears in the case without the alignment limit at tree level, in which case the additional neutral Higgs bosons can mediate the double Higgs boson production process and can provide sizable enhancement of the cross section. We clarify how large enhancement can be obtained in the case without the alignment limit. In the numerical analysis, we focus on a special case of the THDM often referred to as Type-I, because scenarios without the alignment limit are highly constrained by the Higgs boson signal strengths in the other types of Yukawa interactions such as Type-II. Under the current theoretical and experimental constraints on the parameter space, we find a strong correlation between the enhancement of the cross section and the scaling factor of the hff couplings. We note that Higgs boson pair production at e + e − colliders requires the collision energy √ s being larger than 250 GeV while the scaling factor can be precisely determined with √ s = 250 GeV. Hence, the Higgs boson couplings will have been known when experiments reach to the collision energy enough for Higgs pair boson production measurement.
This paper is organized as follows. In Sec. II, we define the model and present the relevant Higgs boson interactions. In Sec. III, we discuss general properties of the Higgs boson pair production in e + e − colliders. Then, we explain a mechanism that enhances the cross section in the THDM without the alignment limit and give a rough estimate of the enhancement factor. Detailed numerical analysis on the aforementioned processes at tree level is performed in Sec. IV. Finally, Sec. V is devoted to our conclusions.

II. MODEL
We consider the THDM whose Higgs sector is constructed by two isospin scalar doublet fields Φ 1 and Φ 2 . For simplicity, we consider the CP-conserving case throughout the paper.
The vacuum expectation values (VEVs) of the two doublets, i.e., It is convenient to introduce the so-called Higgs basis [27,28] as follows: where In this basis, the Nambu-Goldstone bosons G ± and G 0 , which are absorbed into the longitudinal components of the W and Z bosons, are separated from the physical singly-charged Higgs boson H ± and the CP-odd Higgs boson A, while two CP-even Higgs states h ′ 1 and h ′ 2 are generally not mass eigenstates at this stage. By introducing another mixing angle α, the mass eigenstates are given by where we have abbreviated sin X and cos X as s X and c X , respectively. We identify h as the discovered Higgs boson with a mass of about 125 GeV.
The kinetic terms for the Higgs doublets are given as where D µ is the covariant derivative for the Higgs doublets. We note that from the first term on the right-hand side, the gauge-gauge-Higgs type interactions are obtained as From the second term in the right-hand side of Eq. (4), we obtain the Higgs-Higgs-gauge type interactions as with and θ W being the weak mixing angle.
In order to avoid Higgs boson mediating flavour changing neutral currents at tree level, we impose a discrete Z 2 symmetry [15] into the Higgs sector, where the two doublets are transformed as Φ 1 → +Φ 1 and Φ 2 → −Φ 2 . Under the Z 2 symmetry, only one of the two Higgs doublets can couple to each up-type, down-type quarks and charged leptons, by which the interaction matrices in the flavour space between neutral Higgs bosons and fermions are diagonalized in the fermion mass eigenbasis. It has been known that there are four independent types of Yukawa interactions so-called Type-I, -II, -X and -Y [29] depending on the way to assign the Z 2 charge for fermions [30,31].
The Yukawa Lagrangian for the third generation fermions is then given in the Higgs basis by whereΦ = iτ 2 Φ * andΦ ′ = iτ 2 Φ ′ * . The factors ξ b and ξ τ depend on the choice of the types of Yukawa interaction as while ξ t = cot β for all the types of Yukawa interaction. The interaction terms of the Higgs bosons and fermions are extracted as where I t (I b,τ ) = 1/2 (−1/2) and P L (P R ) is the projection operator for the left (right) hand chirality.
The Higgs potential is generally written by eight independent real parameters when we include the soft-breaking term of the Z 2 symmetry: As we already mentioned in the above, we assume the CP-invariance of the Higgs sector, so that m 2 3 and λ 5 parameters are taken to real. After imposing the tadpole conditions, i.e., the requirement of vanishing the linear terms of h and H in the Higgs potential, we can eliminate the m 2 1 and m 2 2 parameters. Then, these eight parameters, i.e., six parameters in the potential and two VEVs are expressed as follows: where M 2 ≡ m 2 3 /(s β c β ). Among these eight parameters, v and m h are fixed to about 246 GeV and 125 GeV, respectively.
From Eq. (10), we can extract the scalar three-point couplings. In particular, the hhh and Hhh couplings will be important in the later analysis, which are expressed in terms of the parameters shown in (11) as follows: where the above quantities are defined by the coefficient in front of the hhh and Hhh vertices in the Lagrangian. Before closing this section, it would be worth to mention that taking s β−α → 1, all the hV V (V = W, Z), hff and hhh couplings become the same values as in the SM at tree level. This limit has been known as the alignment limit, where the SM-like Higgs boson h completely comes from the doublet Φ in the Higgs bases in Eq. (2).

III. DOUBLE HIGGS BOSON PRODUCTION
In this section, we discuss the general property of double Higgs boson production processes at e + e − colliders, and consider how much the production cross section can be different in the THDM with respect to that in the SM prediction. being the total width of H. 2 Therefore, the size of this cross section is typically obtained by multiplying the factor of 16π 2 × c 2 β−α × BR(H → hh) with respect to the cross section of the diagram (a) in the SM, where the factor 16π 2 appears due to the typical ratio of the two-body and three-body phase-space factors, while c 2 β−α comes from the HZZ coupling normalized to the hZZ one in the SM. For example, when s β−α is fixed to be 0.99 (0.995), the above factor becomes 3.14 (1.56) assuming BR(H → hh) = 1. The similar enhancement can also be obtained from the diagram (b') as long as A is produced with on-shell. We thus can expect that the total cross section of the double Higgs boson production can be several times larger than the SM prediction. 1 For the case of √ s larger than 2(m t + m h ) ∼ 600 GeV, the tthh production is allowed and its contribution could be comparable as compared with that shown in Fig. 1. In the numerical analysis of this paper given in the next section, we focus on the collision energy to be 500 GeV, so that this process is not needed to be considered. 2 In our scenario, the typical size of Γ H /m H is 1% level or smaller.
The important thing here is that such enhancement of the cross section happens when departure of the alignment limit, i.e., s β−α = 1, is realized, because the both HZZ and AZh coupling are proportional to c β−α . Therefore, the enhancement is strongly correlated with the deviation in the Higgs boson couplings from the SM prediction. On the other hand, the hV V (V = W, Z) and hff couplings are expected to be precisely measured at future e + e − colliders. For example, at the ILC the hZZ, hW W , hbb, hτ + τ − and hcc couplings may be able to be measured with 0.38, 1.8, 1.8, 1.9 and 2.4% at 1σ level assuming 250 GeV of the collision energy and 2 ab −1 of the integrated luminosity [17]. Therefore, if deviations in the Higgs boson couplings are detected in future, we expect the sizable enhancement of the double Higgs boson production as well.
In the next section, we numerically evaluate how large enhancement can be obtained in the scenario without the alignment limit and show the correlation between the double Higgs boson production cross section and the deviation in the Higgs boson couplings.

IV. NUMERICAL RESULTS
In this section, we compute the cross section of the e + e − → ff hh process, where f is a fermion except for the top quark. These four-body final states are obtained through the Zhh production with the decay of Z into a fermion pair and the vector boson fusion process, see Figs. 1 and 2. We note that when the Z boson from the Zhh process decays into e + e − (ν eνe ), this process interferes with the e + e − hh (ν eνe hh) final states from the Z (W ) boson fusion process. We take into account such interference effects in the numerical analysis.
We also note that the dependence of the type of Yukawa interactions slightly appears in the cross section through the decay width of H and A. In contrast, the type dependence significantly appears in the region of the parameter space allowed from experimental bounds.
In particular, except for the Type-I THDM, the scenario with s β−α = 1 which is considered in this section is highly constrained by the Higgs boson signal strengths, see e.g., [32]. We thus consider the Type-I THDM in what follows 3 . 3 In the Type-I THDM, constraints from flavor experiments such as those from B → X s γ data on the charged Higgs boson mass are also quite milder than those in the other types especially in the Type-II. For example, O(100) GeV of the charged Higgs boson mass is allowed by the B → X s γ data when As seen in (11), there are six free parameters in the THDM, i.e., three masses of extra Higgs bosons, M 2 , tan β and α. Instead of inputting α, we take s β−α and the sign of c β−α as inputs. To avoid large contributions to the electroweak oblique T parameter [34,35], we take the masses of H ± and A to be the same, i.e., m A = m H ± [36][37][38][39]. In this case, the quadratic dependence of the extra Higgs boson masses completely vanishes, and only the small logarithmic mass dependence remains, which is proportional to c 2 β−α . Before going to show the numerical results for the cross section, let us summarize the constraints on the parameter space what we take into account in the analysis. We impose the perturbative unitarity bound [40][41][42][43] and the vacuum stability bound [36][37][38][39] as the constraints from theoretical consistency. These bounds restrict the size of scalar quartic couplings in the potential, which can be translated into the bound on the Higgs boson masses and the mixing angles. As the experimental constraints, we take into account the electroweak oblique S and T parameters [44], direct searches for additional Higgs bosons at the LEP, Tevatron and LHC experiments as well as the compatibility of the signal strengths for the discovered Higgs boson with a mass of 125 GeV. For the direct searches, we use the HiggsBounds-5.3.0beta [45]. For the signal strengths of the discovered Higgs boson, we use the combined data from ATLAS and CMS at the LHC Run-I experiments [3], and require that the prediction of the signal strengths for the γγ, ZZ * , W W * and τ τ modes of the h state does not exceed the given width of the error bar. We impose these experimental bounds at the 95% confidence level.
In Fig. 3 in the THDM normalized to the SM value g SM hf f , see Eq. (7): We note that due to the choice of the Type-I THDM, the scaling factor κ f does not depend on the choice of a fermion f , see Sec. II. From Eq. (14), we see that κ f < s β−α (κ f > s β−α ) tan β > ∼ 2 [33] in the Type-I THDM, while m H ± < ∼ 600 GeV is excluded with 95% confidence level in the Type-II THDM.
is obtained by taking the sign of c β−α to be negative (positive), and κ f = s β−α is given in the limit of tan β → ∞.
We Now we present numerical results on the cross section. We consider the ratio of the di-Higgs production cross section in the Type-I THDM to that in the SM: where the summation for f is done over all the fermions except for f = t. For the numerical computation of the cross section, we used the public version of GRACE [46][47][48] with some modifications. All the calculations were performed at tree level with √ s = 500 GeV and 1 ≤ tan β ≤ 30.
The scatter plot in Fig. 4 shows the ratio R as functions of the scaling factor κ f with using sampling points passing all the constraints, namely, points in the white region of Fig. 3.   We find that the value of R can be maximally around 5 (3) for s β−α = 0.99 (0.995).
We also perform the similar calculation shown in Figs. 5 and 6, but for the case with m H = 300 GeV and m A = m H ± = 500 GeV. The region allowed by the constraints is almost the same as the previous plots in Fig. 3. On the other hand, by looking at Fig. 6, we find that the shape of the scattering points are shifted to below due to the A mediation being off-shell, and the maximally allowed value of R becomes around 4 and 2.5 for s β−α = 0.99 and 0.995, respectively.
Finally, we show the mass dependence of the value of R in Fig. 7 with a fixed value of tan β to be 2 (top), 3 (middle) and 5 (bottom). In these plots, we take the degenerate mass of the extra Higgs bosons, i.e., m H = m A = m H ± , and varying the mass range to be from 200 GeV to 500 GeV. Again, all the points are passed all the constraints explained in the above. For the low tan β case such as tan β = 2 shown in the top panels, the constraint from direct searches excludes the lower mass region, so that the points appear only in the larger mass region. For the larger tan β case particularly with c β−α > 0 (left panels), the lower mass region is allowed, and it is clearly seen that the value of R drastically increases at the mass of extra Higgs boson just above 215 and 250 GeV, because of the on-shell A → Zh and H → hh decays open. We also find that for the larger mass region, e.g., m H > ∼ 400 GeV the value of R becomes smaller than 1, in which the H and A appearing in the Zhh production are getting off-shell, so that these contributions become unimportant. The value R < 1 is simply explained by the fact that the contribution containing h * → hh as the diagram (a) in Fig. 1 becomes smaller than the SM prediction because of the smaller λ hhh coupling with respect to the SM value, see Eq. (12). We note that in these plots, the value of tan β is fixed, so that the value of κ f is determined to be (1.06, 1.04, 1.02; 0.92, 0.94, 0.96) for s β−α = 0.99 Although expected not to give so drastic change, radiative corrections to the cross-section enhancement studied in this work may be investigated by using, for example, H-COUP [49] to incorporate one-loop electroweak-corrected vertices or GRACE-loop [50,51] for full one-loop computation. We leave it for future works.