Higgs boson pair production through gauge boson fusion at linear colliders within the general 2HDM

Inclusive Higgs boson pair production through the mechanism of gauge boson fusion e^{+} e^{-} ->V* V* ->h h + X (V=W,Z) in the general Two-Higgs-Doublet Model (2HDM), with h=h^0,H^0,A^0,H^{\pm}, is analyzed at order \alpha^4_{ew} in the linear colliders ILC and CLIC. This kind of processes is highly sensitive to the trilinear Higgs (3H) boson self-interactions and hence can be a true keystone in the reconstruction of the Higgs potential. For example, in the ILC at 1 TeV, the most favorable scenarios yield cross-sections up to roughly 1 pb, thus entailing 10^5 events per 100 fb^{-1} of integrated luminosity, whilst remaining fully consistent with the perturbativity and unitarity bounds on the 3H couplings, the electroweak precision data and the constraints from BR(b->s\gamma). Comparing with other competing mechanisms, we conclude that the Higgs boson-pair events could be the dominant signature for Higgs-boson production in the TeV-class linear colliders for a wide region of the 2HDM parameter space, with no counterpart in the Minimal Supersymmetric Standard Model. Owing to the extremely clean environment of these colliders, inclusive 2H events should allow a comfortable tagging and might therefore open privileged new vistas into the structure of the Higgs potential.


Introduction
For more than 40 years, the Standard Model (SM) of Strong and Electroweak interactions has furnished a successful arena in which to describe the physics of Elementary Particles. In spite of its formidable achievements, a number of longstanding challenges are still to be resolved. Perhaps the most paradigmatic one concerns the ultimate nature of Electroweak Symmetry Breaking (EWSB). Even though the Higgs mechanism provides an elegant description of EWSB within a perturbative quantum field theory framework, one is forced to postulate the existence of (at least) one scalar (J P = 0 + ) fundamental building-block of Nature, whose experimental confirmation is conspicuously missing for the time being. Nevertheless, with the recent start-up of the LHC operations at CERN, the prospects for the discovery of the Higgs boson are extremely encouraging and the curtains may soon be drawn back to reveal this final player in the story of the SM. In fact, the LHC will be able to amply sweep the natural SM Higgs mass range (spanning from a few hundred GeV up to ∼ 1 TeV) and, in the most favorable scenarios, it could produce a copious number of Higgs boson events [1].
Despite these optimistic prospects for the discovery of the SM Higgs boson, the upcoming LHC data might well reveal that the ultimate origin of EWSB is grounded somewhere beyond the minimal assumptions of the SM. For instance, a number of well-motivated (perturbative) extensions of the SM entail a two-Higgs-doublet structure, such as the Minimal Supersymmetric Standard Model (MSSM) [2] or the general (unconstrained) Two-Higgs-Doublet Model (2HDM) [3]. From the beginning, however, the quest for experimental signatures of Higgs boson physics beyond the SM concentrated its efforts mainly on the search for supersymmetric (SUSY) Higgs bosons [2,3]see [4] for more recent developments, and [5] for fresh reviews of this subject. At the same time, very detailed investigations were undertaken in the past on the SUSY quantum effects on the gauge boson masses [6] (cf. [7] for the current state of the art) and on the Z-boson and top quark widths [8], in which virtual Higgs bosons may also play a significant role. Later on, outstanding new sources of quantum corrections were identified in processes involving Higgs bosons, quarks and/or squarks as external particles. The possible new effects stemmed from the enhancement capabilities associated to the supersymmetric Yukawa couplings between Higgs bosons and quarks or between quarks, squarks and chargino-neutralinos. In a variety of quite different processes, the potential effects were shown to be truly spectacular, see e.g. [9][10][11][12]. More recently, it has been shown that such an enhancement also applies to the loop-generated Yukawa couplings of additional gauge-singlet Higgs bosons in minimal extensions of the MSSM [13].
So far so good, but what about the non-SUSY Higgs boson extensions of the SM? Remarkably, the crucial novelty here is that the most relevant effects could have an origin fundamentally different from the Yukawa couplings linked to the traditional SUSY sources of enhancement. This possibility has recently been illustrated for the general 2HDM in [14] (see also [15] in a different domain). In this letter, we dwell on another compelling (and complementary) example of this fact, which we encounter once more in the physics of the Higgs boson self-couplings at linear colliders. Assuming that the LHC uncovers a neutral Higgs boson, a critical issue will be to discern whether this particle is compatible with the SM and/or any of its extensions and, in the latter case, to which of these extensions it belongs. The next generation of TeV-class linear colliders (based on both e + e − and γγ collisions), such as the ILC and CLIC projects [16], will be invaluable in answering this fundamental question. Owing to its extremely clean environment, a linear collider should allow for precise measurements of: i) the Higgs boson mass (or masses, if more than one); ii) the couplings of the Higgs bosons to quarks, leptons and gauge bosons; and iii) the Higgs boson self-couplings mentioned above, i.e. the trilinear (3H) and quartic (4H) Higgs boson self-interactions. In a nutshell, it should allow us a keen insight into the physics lying beneath the EWSB mechanism and, ultimately, to reconstruct the Higgs potential itself.
In this work, we provide some complementary clues to this reconstruction process: specifically, we focus on the effects of the trilinear Higgs boson couplings on the inclusive production of Higgs boson pairs induced by weak gauge boson fusion, i.e. e + e − → V * V * → h h + X, where V = W ± , Z and h = h 0 , H 0 , A 0 , H ± . We show that this mechanism could be the leading Higgs boson production channel at TeV energies and, if so, it should furnish an intrinsic and totally unambiguous signal of non-supersymmetric new physics. Therefore, while the physics of the top (and bottom) quark [17] is the natural playground for the study of the SM and MSSM Higgs boson interactions, we find that in the case of non-SUSY theories there are comparable (if not better) opportunities in sectors of the theory not necessarily related to heavy quark flavors but to the very structure of the Higgs potential.

Higgs boson production induced by trilinear Higgs interactions
Of cardinal importance is to understand in detail the phenomenology of the Higgs sector in the context of linear colliders (both within the SM and its renormalizable generalizations). Quite an effort has already been devoted to this goal in the literature. Exclusive double (2H) and multiple Higgs boson production, for instance, has been comprehensively investigated, although mainly within the MSSM [18][19][20]. The exclusive 2H case consists of the simplest Higgs boson production processes: Similarly, the two-body final states e + e − → Zh and e + e − → Ah (with h = h 0 , H 0 ) have been long known to be complementary to each other in the MSSM [18]. Notice that processes with exclusively two identical neutral Higgs bosons in the final state cannot proceed at the tree-level (neither in the SM nor in the MSSM), and at one-loop they have very tiny cross-sections of order 10 −5 pb [14]. However, sizeable rates of two Higgs bosons in the final state in an e + e − collider may arise from the processes like (2.1), whose detection would signify an unmistakable observation of physics beyond the SM. Nevertheless, let us highlight that none of the exclusive 2H channels (2.1) is sensible to the trilinear Higgs boson couplings at the leading order. Within such pairwise Higgs events, therefore, probing the structure of the 3H self-interactions would only be possible through the analysis of the radiative corrections. Indeed, one needs to include such quantum effects in the exclusive 2H boson processes (2.1) so as to disentangle the genuine imprints of a SUSY Higgs sector from a non-SUSY one. In this context, a rich literature exists on the one-loop calculation of the cross-sections for the two-particle final states within the MSSM 1 . There are also studies considering radiative corrections to charged Higgs production in e + e − collisions within the 2HDM [25], and also on single, double and multiple Higgs production at the LHC but mostly for the MSSM [26]. A complete 1-loop analysis of the exclusive 2H channels (2.1) in the general 2HDM is missing, however, and will be presented elsewhere [27].
Our main endeavor in this letter is to further investigate the foreseeable impact of these 3H self-couplings in a class of processes which critically depend on them already at the tree-level. This was for instance the case of Ref. [14], in which the triple Higgs boson couplings were probed by performing a systematic analysis of the exclusive production processes with three Higgs bosons in the final state. There are three classes of processes of this kind compatible with CP-conservation, to wit: where, in class 2), we assume that the two neutral Higgs bosons h must be the same, i.e. the allowed final states are ( . The cross-sections in the 2HDM were shown in [14] to reach up to O(0.1) pb, i.e. several orders of magnitude over the corresponding MSSM predictions [18]. Similar effects have also been recently reported in 2H strahlung processes of the guise e + e − → Z 0 h h [28], and also in the loop-induced 2H production through γγ interactions [29].
It is important to note that, in a CP-conserving theory, all the 3H final states in (2.2) contain at least one charged or pseudoscalar Higgs boson, and this has practical implications. In fact, let us recall that there is an important distinction between the two basic types of generic 2HDM models [3]; namely, whilst light charged Higgs bosons are possible within a type-I 2HDM -in which only the Φ 2 field couples to fermions 2 -, in type-II models there are important constraints on the charged Higgs mass due to contributions to the flavor-changing neutral-current (FCNC) process b → sγ [30,31]. The Higgs pseudoscalar A 0 is then also constrained to be relatively heavy M A ∼ M H ± , due to the limits on δρ -see e.g. Eq. (2.6) of Ref. [14].
Triple Higgs boson couplings may drive different kinds of processes. The phenomenological impact, however, may seriously depend on whether the underlying Higgs sector belongs to the MSSM or to the general 2HDM. In practice, this means that we should be ready to identify highly distinctive signatures of the underlying model. We have mentioned above that the trilinear Higgs boson couplings are involved at the tree-level in processes with three Higgs bosons (3H) in the final state, see (2.2). But, in fact, they are also involved in inclusive processes with two Higgs bosons (indicated as 2HX) in the final state. For instance, the 3H-coupling has been investigated phenomenologically in TeV-class linear colliders in [19,22,23] through the double-Higgs strahlung process e + e − → HHZ or the WW double-Higgs fusion mechanism e + e − → H + H − ν e ν e . These processes, which include vertices like ZZH, WWH, ZZHH, WWHH and HHH, are possible both in the SM and its extensions, such as the MSSM and the general 2HDM. Unfortunately, the crosssection turns out to be rather small both in the SM and in the MSSM, being of order 10 −3 pb at most, i.e. equal or less than 1 fb [22]. Even worse is the situation regarding the 3H boson production in the MSSM, in which -except in the case of some particular resonant configurationthe typical cross-sections just border the value ∼ 0.01 fb or less [22]. In the latter reference it has been shown that, if the double and triple Higgs production cross-sections would yield sufficiently high signal rates, the system of couplings and corresponding double/triple Higgs production crosssections could in principle be solved for all trilinear Higgs self-couplings up to discrete ambiguities using only these processes. But in practice these cross-sections are manifestly too small to be all measurable.
In stark contrast with the pessimistic situation in the MSSM, the unconstrained 2HDM may exhibit quite promising signatures. A key point here is the potentially large enhancements that the 3H couplings may accommodate -unlike the MSSM case, where such self-couplings are constrained by SUSY invariance and are all predicted to be purely gauge [3,5]. To make it transparent with a single explicit example, the analytical expression for the particular h 0 h 0 H 0 self-coupling in both models reads 3 It is patent from these expressions that the 2HDM coupling can be enhanced both at low and high tan β, and also through the Higgs boson mass splittings, whereas the MSSM coupling cannot be enhanced in any way. In practice, the enhancement possibilities of the 2HDM couplings will be partially tamed by the aforementioned set of phenomenological and theoretical restrictions.
In the next section, we analyze the inclusive Higgs boson-pair production at linear colliders within the general 2HDM, where X = (e + , e − ;ν e , ν e ). At high energies (∼TeV) the vector boson fusion diagrams of the kind displayed in Fig. 1 constitute the dominant mechanism. Therefore, in practice the bulk mechanism of (2 Notice that the fusion mechanism may trigger either one or no Higgs boson as a virtual intermediate state. There is also the possibility that a single real Higgs boson is produced on-shell and subsequently decays into two Higgs bosons of smaller mass. In this respect, let us recall that, in the SM, the mechanism of single Higgs boson production via gauge boson fusion in e + e − collisions was already considered long ago [33] and was shown to be dominant with respect to the annihilation channels at high energy. Still, the cross-sections are very small in the SM (of order 1 − 10 fb) for masses M H ∼ 100 GeV and they were computed at that time for the "future" LEP energies. Some enhancement can be achieved for charged Higgs production in extended Higgs models [34]. However, none of the last two sorts of processes involve the Higgs boson self-couplings. In the following, we concentrate on the computation of the cross-section for the processes (2.4), which do involve the 3H-couplings in some of their Feynman diagrams, see Fig. 1. It turns out that these specific diagrams are responsible for the bulk of the cross-sections under the most favorable conditions for these processes. We will perform the calculation in the general 2HDM model and shall take the same set of phenomenological restrictions used in Ref. [14], to which we also refer the reader for more details on the relevant pieces of the interaction Lagrangian.

Double Higgs boson production from weak gauge boson fusion
In contrast to the simple 2H channels (2.1), the class of triple Higgs boson processes (2.2) and the inclusive 2H processes (2.4) are both directly sensitive to the trilinear self-interactions, which implies both a source of potential enhancement and a possible strategy to measure such couplings. In Ref. [14] it was shown that, for Higgs boson masses 100 GeV , the production cross-sections corresponding to the 3H processes (2.2) could be remarkably high in the general 2HDM, lying typically above the inclusive 2H cross-sections (2.1) at center-of-mass (CM) energies 1 TeV. This feature, which is impossible in the MSSM, is brought about by the enhanced Higgs boson self-couplings involved in the triple Higgs production mechanism (2.2).
In a similar way, when we move to the 2HX processes (2.4) and consider higher and higher energy (typically at the ∼ TeV range of linear colliders), several opportunities for significant enhancement may appear. The leading mechanism behind these processes is the weak gauge-boson fusion mechanism Fig. 1 for the case h = h 0 ). As a result, the three main sources of enhancement here are the following: i) first, the t-channel gauge boson fusion is not so severely damped at high energies as compared to s-channel annihilation; ii) second, the triple Higgs vertex is involved in some of the gauge boson fusion channels (see Fig. 1, first and second diagrams from the left in either row); and, finally, iii) in some cases, the virtual Higgs boson produced in the last sort of diagrams may not be too far away from the resonance and, in such circumstance, the 2HX final state can be more copiously produced. In particular, if the Higgs boson H 0 is heavy enough, the resonant decay H 0 → h 0 h 0 will be kinematically allowed. Although this decay mode is also possible within the MSSM, the production of the initial H 0 is mostly suppressed in this model, since its gauge couplings are known to be complementary to those of the SM-like h 0 . The rate for this process in the MSSM is therefore not competitive with the 2HDM one. The upshot is that some of the gauge boson fusion processes (3.1) can become fully competitive, if not the leading mechanism of Higgs boson production, at high energy in the general 2HDM. Let us emphasize that the cross-sections for the processes (3.1) grow up to very high values of the CM energy √ S. This is possible because, for the fusion diagrams, the weak gauge bosons V can be quasi-real and hence have virtual momenta well-below the CM energy of the process, which may satisfy S ≫ M 2 V -the rest of the energy being carried away by the concomitant lepton final states X = (e + , e − ;ν e , ν e ). Therefore, while the exclusive 2H and 3H production cross-sections are expected to decay irremissibly as ∼ 1/S owing to the Z-boson propagator that mediates the schannel diagrams, the energy behavior of the gauge fusion mechanism for 2HX production is quite different. It is actually reminiscent of the cross-section for two-photon processes e + e − → γ * γ * → Y + e + e − , which in the asymptotic limit goes roughly as ∼ (α 4 /M 2 ) ln 2 (S/m 2 e ) ln n (S/M 2 ), where M is the threshold mass of the produced final state Y and the number n ≥ 1 depends on the high energy behavior of σ(γγ → Y ) 4 . In our case, the situation is a bit more complicated because we have massive gauge bosons V (rather than photons), and moreover the diagrams with triple Higgs vertex contain an intermediate virtual Higgs state emerging from the V V -fusion. Thus, the dependence of the cross-sections (2.4) with the threshold mass of the Higgs boson pair is not so simple, but as in the two-photon case the cross-sections are expected to raise logarithmically with the energy, rather than decaying quadratically with it. We shall confirm these expectations with the numerical analysis in the next section.
The following comment is also in order. Even though the results on the inclusive Higgs bosonpair production (2.4) are overwhelmingly dominated by the gauge boson fusion mechanism (3.1), we point out that we have computed the cross-section of the processes (2.4) with full generality, i.e. by including all the diagrams at order O(α 4 ew ). This is actually necessary in order to insure the gauge invariance of the overall result. Some subsets of diagrams are nonetheless completely irrelevant, such as e.g. those Z-mediated amplitudes where the Higgs bosons are radiated off the lepton legs. However, there are other more subtle subsets. In particular, there are annihilation diagrams leading to the same final state (2.4) in which a pair (Z * , h * ) -consisting of a virtual gauge boson and a virtual Higgs boson -is produced and, subsequently, Z * decays into lepton pairs and h * radiates a real Higgs boson. Although none of these topologies gives a dominant effect in front of the primary gauge boson fusion mechanisms, a careful treatment is nevertheless required in this case so as to preserve the consistency of the overall procedure. In particular, the aforementioned Higgs strahlung process demands the introduction of a Breit-Wigner propagator for the virtual bosons.
The above considerations suggest that large production cross-sections for the inclusive 2HX processes (2.4) should be possible at high energy regimes when both the exclusive 2H and 3H channels (2.1) and (2.2) -all of them mediated by Z-boson exchange -are kinematically suppressed as ∼ 1/S. We have undertaken this calculation using the computational tool CompHEP [36]. In some steps of the computation, we have also made use of FeynArts and FormCalc [32], which served also for cross-checking purposes. Furthermore, to ease the comparison with the existing analysis of the 3H processes (2.2), we present our numerical calculation on the basis of the same set of free parameters as in [14], i.e.
This means that the general 2HDM Higgs potential is treated under the assumption that λ 5 = λ 6 (in order not to depart too much from the the MSSM structure of the Higgs sector, see [14] for details). For a realistic computation, one has to include all known existing constraints. For example, there are strong limitations on the parameters of a type-II 2HDM coming from FCNC processes, mainly from the charged Higgs boson contributions to B(b → sγ), which would overshoot the allowed experimental range unless M H ± 350 GeV [30]; and there are also the radiative corrections to δρ from the 2HDM sector, which must not exceed the limit |δρ 2HDM | ≤ 10 −3 . Finally, there are of course the direct search limits from LEP and the Tevatron [31]. All these constraints are integrated in our codes and, thus, the entire numerical analysis is consistent with these bounds, along with the general unitarity constraints which apply to both type-I and type-II 2HDM. More details on these constraints are discussed in [14], to which we refer the reader for additional information. Let us also note that, in order to obtain more accurate results, a running value for the electromagnetic coupling constant α(M Z ) = 1/127.9 has been used.

Non-resonant double Higgs boson production
For the numerical analysis, we will consider five different sets (I-V) of parameters of the general 2HDM Higgs sector, see Table 1. This should suffice to illustrate the enhancement possibilities of the inclusive 2HX cross-section (3.1). Notice that sets I-II and V are compatible with the type-II 2HDM (because the charged Higgs boson mass is sufficiently heavy to satisfy the aforementioned constraints). These sets are actually possible for type-I models too. On the other hand, the relatively light Higgs boson mass sets III-IV will be used (exclusively) for the less constrained type I models. Let us also remark that in the case of sets I-III the resonant decay H 0 → h 0 h 0 is not possible, although for Sets I and III the mass threshold for such resonant mode is closer than for Set II. Finally, Sets IV and V explore the possibility of having some resonant Higgs boson decays and are mainly intended for type-I and type-II models, respectively. (Set V is also compatible with type-I models, as remarked before, but we mainly aim at type-II ones for that set We consider first the numerical results obtained from the non-resonant scenarios, which are more general. In doing so, we wish to compare the cross-sections for both the 3H and 2HX channels (2.2) and (2.4), respectively. In Fig. 2a, we show these production cross-sections as a function of the CM energy for Set I of Higgs boson masses (full lines in that figure). The fixed value of the CPeven mixing angle in this figure (sin α = −0.7) has been determined in combination with tan β after maximizing the cross-sections for the mass Set I under the various constraints discussed previously. The relevant trilinear Higgs self-couplings h 0 h 0 h 0 and H 0 h 0 h 0 , and hence the overall production cross-section, become maximally enhanced at tan β ≃ 25. The main constraint that fixes the aforementioned tan β value is the unitarity bound of the trilinear Higgs couplings. In order to explore the effect on the cross-sections after significantly increasing some masses (while respecting all the necessary bounds mentioned above), we also show the corresponding 2HX and 3H crosssections for M A 0 = M H ± = 800 GeV, with all other parameters fixed as in Set I (dashed curves in that figure). In this case, the unitarity constraints pull the maximum value of tan β down to S ≡ Ecm (in GeV) for the particular 2HX production process e + e − → h 0 h 0 + X, together with the sum of all the exclusive 3H channels, e + e − → 3H for a choice of Higgs masses as in Set I (in full lines). We fix sin α = −0.7 and tan β = 25, in which case the relevant 3H couplings, and hence the production cross-section, are maximally enhanced. The 2HX and 3H cross-sections are also depicted (in dashed lines) as in the previous case but for M A 0 = M H ± = 800 GeV; here, however, tan β ≃ 4 in order to preserve the consistency with the unitarity bounds. In the right axis, we track the predicted number of events per 100fb -1 of integrated luminosity; b) As before, but for a lower enhancement of the 3H couplings: tan β ≃ 14 with Set I masses (in full lines) and tan β ≃ 2 with M A 0 = M H ± = 800 GeV (in dashed lines). Fig. 2 is that the kinematical threshold for the overall 3H contribution shifts to higher energies when we boost M A 0 up to 800 GeV. Let us clarify that the rise of further thresholds does not leave a visible footprint on σ( √ S), provided one of the 3H channels (in the present case h 0 h 0 A 0 ) dominates over the remaining ones. The subdominant channels, although have less enhanced 3H couplings and a larger phase-space suppression, contribute to smooth out significantly the damping of the total 3H cross-section as a function of √ S as compared to the individual channels.

∼ 4. Noteworthy in
In Fig. 3, we show the numerical analysis of the 2HX cross-section for various values of tan β, with the Higgs mass parameters taken as in Set I (Fig. 3a) and Set II (Fig. 3b). The corresponding 3H production cross-sections are shown for comparison in the lower panels (Figs. 3c and 3d). Likewise, in Fig. 2b we display the corresponding results obtained for a lower enhancement of the 3H couplings (viz. up to one half of the standard unitarity bound used in [14]). Notice that our energy scan actually sweeps a wide range which reaches up to 5 TeV in order to better show the asymptotic trend of the cross-sections. In practice, the ILC will cover only the approximate range 0.5 − 1.5 TeV, whereas CLIC may reach 3 TeV [16].
The dominant effect on the inclusive 2HX amplitudes is proportional to the Higgs trilinear couplings H 0 h 0 h 0 and h 0 h 0 h 0 . These couplings are enhanced at large values of tan β, as can be seen from the first Eq. (2.3) and in general from Table 1 of Ref. [14]. Let us also emphasize that, in order to better appreciate the tan β-dependence, the maximum 2HX and 3H cases corresponding to tan β = 25 are again included in Figs. 3a,c along with the other (smaller) tan β values. Although the overall production rates decrease with decreasing tan β, the 2HX channel remains dominant at (and above) √ S = 1 TeV for Set I. In this same energy range, but for Higgs boson masses as in Set II, the 2HX channel remains comparable to the 3H channel only for the largest allowed values of tan β. At higher energies, however, such as those planned for CLIC, the 2HX channels become the S (in GeV) for the inclusive Higgs boson-pair production process e + e − → h 0 h 0 + X (upper panels) and for the sum of all the exclusive 3H channels, e + e − → 3H (lower panels). In Figures a) and c), Set I of Higgs boson masses is employed, while Set II is used for b) and d) -See Table 1. As in the previous figure, we take sin α = −0.7 and we study the behavior of the cross-section over different values of tan β. In the right axis of each plot, we track the predicted number of events per 100fb -1 of integrated luminosity.
leading ones even at small values of tan β. To be sure, some of these gauge fusion processes furnish a very competitive strategy to probe the 3H self-couplings. This strategy nicely complements the prospects that were singled out for the exclusive 3H channels in Ref. [14], particularly at very large center-of-mass energies, where the 3H signal is greatly suppressed whilst the 2HX one remains sustained and even logarithmically enhanced. Cross-sections for the remaining 2HX final states containing heavier Higgs bosons are not shown in Fig. 3 as they are found to have negligible production rates when compared to the corresponding 3H final states in this scenario. Explicitly, the maximum production rates are found to be of the order 10 −2 pb for the h 0 H 0 case, and below 10 −3 pb for the rest of the 2HX final states.
Notwithstanding, a starkly distinct panorama shows up for lighter H ± and A 0 bosons. In Figure 4a, we display σ( √ S) corresponding to Set III of Higgs boson masses (see Table 1), in which case sizable cross-sections are attained for a number of 2HX channels: H + H − yields ∼ 1 pb; h 0 h 0 and A 0 A 0 give ∼ 0.1 pb, and h 0 H 0 renders ∼ 0.01 pb. Owing to the relatively light mass of the charged Higgs boson, the latter scenario is only allowed for type-I, not for type-II, 2HDMotherwise it would yield an exceedingly large FCNC contribution to B(b → sγ) [31]. Similarly, in S (in GeV) for the 2HX processes e + e − → hh + X for Sets III and IV of Higgs boson masses, which are suitable for type-I 2HDM -see Table 1. Displayed are the Higgs boson-pair channels whose cross-sections lie above 0.01 pb, together with the corresponding triple-Higgs production rate, e + e − → 3H, for each scenario. The values chosen for sin α and tan β are quoted in the figures. Fig. 4b we present the corresponding results for Set IV of Higgs boson masses, which also applies exclusively for type-I models, but refers to the resonant situation in which the on-shell decays H 0 → h 0 h 0 , A 0 A 0 , H + H − are all allowed. In the next section, we further dwell on the resonant case, but we consider a scenario (valid for both type-I and type-II models) where only the final state with two light CP-even Higgs bosons is permitted, i.e. H 0 → h 0 h 0 , and we study it in more detail.

Resonant double Higgs boson production
We note that the scenarios considered above are proper of the general 2HDM. In the MSSM case, where the masses of the CP-even light and heavy Higgs bosons h 0 , H 0 are not independent parameters once tan β and M A 0 are given [3], the mass splittings indicated in Table 1 are not possible. For larger enough values of M H 0 , there is a drastic change in the behaviour for the production cross-section of the inclusive channel e + e − → h 0 h 0 + X since the on-shell decay H 0 → h 0 h 0 becomes kinematically available. Indeed, with M H 0 > 2M h 0 the cross-sections are somewhat less dependent upon the enhancement of the Higgs trilinear couplings, the reason being that highly enhanced trilinears now lead to a dramatic broadening of the H 0 resonance with a branching ratio essentially of order one.
In Figures 5a-5d, we exhibit a panoply of 2HX and 3H production cross-sections for the Higgs boson masses as in Set V of Table 1. We explore the dependence with the CM energy and the mixing angles of the Higgs sector. In all cases, it corresponds to a situation where the on-shell decay H 0 → h 0 h 0 is possible but, in contradistinction to Set IV considered in Fig. 4b, the alternate decays H 0 → A 0 A 0 , H + H − are not allowed. The case under consideration is more along the line of type-II models, which are the the closest ones to the SUSY case. From these figures, it is patent that the inclusive cross-section exceeds the 3H channel for all energies above the TeV scale. Furthermore, there is a significant dependence of the 2HX cross-section on sin α, which enters through the vertex factor cos 2 (β − α) associated to the on-shell production of H 0 from W + W − and Z 0 Z 0 fusion (cf. Fig. 1). At the same time, the 2HX production cross-sections are now largely independent of tan β, because changing this parameter just leads to small fluctuations of the H 0 branching ratio around one. At variance with this mild tan β dependence of the 2HX processes, the main enhancement source of the 3H final states still resides in the tan β effective dependence of the Higgs trilinears, as can be seen in Figure 5d.
Interestingly enough, let us emphasize that the potentially large values of the 2HX cross-sections studied up to now, either with resonant or non-resonant configurations, are a kind of trademark prediction of the non-supersymmetric two-Higgs-doublet models. Of course, the collection of diagrams shown in Fig. 1 also accounts for the corresponding 2HX processes within the MSSM. Nonetheless, the 3H couplings are constrained by supersymmetry and are directly related to the electroweak gauge couplings; there is therefore no possibility of enhancement. Furthermore, in the MSSM the region of parameter space where the relation M H 0 > 2 M h 0 holds is dominant. As a consequence, the inclusive 2HX production will be brought about predominantly by the production of on-shell H 0 bosons via W + W − fusion (cf. Fig. 1) and followed by their subsequent on-shell decay H 0 → h 0 h 0 . As already mentioned, it is precisely the W + W − H 0 ∼ cos(β − α) coupling that limits the resonant 2HX production in the MSSM because supersymmetry constrains the lightest CP-even Higgs boson h 0 to have SM-like couplings when the remaining partners are heavy. Consequently, the resulting cross-sections are expected to lie far below the optimal 2HDM yields. To illustrate these features in a concrete way, in Table 2 we compute the predicted MSSM cross-sections for the inclusive production of a h 0 h 0 pair, σ(e + e − → h 0 h 0 + X). Concerning the CM energies, we assume the operation range that is scheduled for either the ILC ( √ S = 0.5 − 1.5 TeV) and CLIC (up to √ S = 3 TeV), and take three different CP-odd Higgs boson masses M A 0 = 200 , 300 and 500 GeV. The value of tan β is not shown, but it is determined such that to (approximately) optimize the corresponding production rate. In calculating these values, we have taken all soft SUSY-breaking masses equal to M SUSY = 1 TeV, along with µ = 200 GeV and A t = A b = A τ = 1 TeV 5 . From Table 2, we see that the 2HX rates render a contribution of order of a few fb , at most, for M A 0 > 300 GeV. Nevertheless, there is a narrow corner in the region of low M A 0 /low tan β (which is severely restricted by the LEP mass bounds [31]) wherein the W + W − H 0 coupling is not so much hampered and hence the cross-section climbs up to ∼ 10 fb. By comparing the (optimal) values in Table 2 to the most favorable scenarios displayed in Figures 2-5, we conclude that the 2HX signal in the 2HDM is typically a factor 10 − 100 larger than its MSSM counterpart.
Finally, we have also evaluated the "SM background", i.e. the cross-section for double Higgs production from gauge fusion in the SM for the same ∼TeV energies under study. We find e.g. that σ(e + e − → V * V * → H H + X) = 10 −3 − 10 −5 pb for SM Higgs boson masses in the range M H = 115 − 300 GeV (see Table 3), which is relatively quite small as compared to our favorite 2HDM scenarios. Furthermore, we note that the produced SM Higgs boson will predominantly decay into gauge boson pairs W + W − , Z 0 Z 0 and (to a lesser extent) to bottom quark pairs bb and top quark pairs tt (if kinematically possible). Subsequently, the gauge bosons will decay both into leptonic and light quark channels. In contrast, in the 2HDM case, the very same conditions that favor the fusion production of Higgs bosons do also favor the decay of the produced Higgs boson into other Higgs bosons and these into heavy quarks. So the kind of signature is very distinct. Another source of background to the 2HDM signal could come from gauge boson fusion into gauge bosons, essentially Z 0 pairs. If these gauge bosons subsequently decay into quarks, this would mimic the Higgs boson themselves decaying into quark pairs. However, explicit calculation shows that the cross-section σ(e + e − → V * V * → Z 0 Z 0 X) reaches up to 0.1pb at most (Cf.  Table 3: Cross-sections for the leading background processes within the Standard Model, these are the SM Higgs boson pair-production, e + e − → H SM H SM + X, and the Z 0 boson pair production, e + e − → Z 0 Z 0 + X. The corresponding production rates are computed for MH SM = (115, 300) GeV in order to account for phenomenologically relevant scenarios. For a sufficiently heavy SM Higgs, resonant production of Z 0 -boson pairs occurs.

Conclusions
We have devoted this work to the study of the inclusive production of Higgs-boson pairs in e + e − collisions (2.4), mainly triggered by the weak gauge boson fusion mechanism (3.1) within the general 2HDM. We have found that the cross-sections for some of these processes can be several orders of magnitude larger than their MSSM counterparts. Moreover, in contrast to the two-body final states (2.1), a tree-level analysis of the gauge fusion mechanism could reveal the general 2HDM nature of the Higgs bosons involved, if the enhancement properties that we have explored effectively apply in the physical region of the parameter space. This was shown to be the case also for the previously considered processes (2.2) with three Higgs bosons in the final state [14]. However, the inclusive 2HX channels (3.1) could be by far the leading mechanism for Higgs boson production at the characteristic TeV energies of the planned ILC and CLIC colliders [16]. We find remarkable opportunities whose threefold origin stems from: i) the sustained (logarithmic) growing of the weak gauge boson fusion channels with √ S at TeV energies; ii) the enhanced regime of the trilinear (3H) Higgs boson couplings in the 2HDM; and iii) the possible resonant (or near resonant) decay of an intermediate Higgs boson into the final Higgs boson-pairs.
Phenomenologically interesting results are attained in large regions of the 2HDM parameter space. In these domains, the maximum cross-section of the inclusive production of Higgs boson pairs (2.4) can be far above (one to two orders of magnitude) the exclusive triple-Higgs boson events (2.2). The possible sources of background (coming from SM Higgs-and gauge-boson pair production) are mostly negligible. In the case of type-II models, which are closer to the MSSM Higgs sector, two simultaneously of the inclusive processes (2.4) can have cross-sections above 0.01 pb (and one of them at the level of 1 pb), therefore with production rates at the level of 10 3 − 10 5 events per 100 fb -1 of integrated luminosity. Let us note that type-I models may also lead to crosssections of order of 1 pb in particular channels, but here we may have up to four channels (2.4) simultaneously having sizeable cross-sections of 0.01 − 0.1 pb, thus amounting to production rates of 10 3 − 10 4 events for the same segment of integrated luminosity. In contrast, the corresponding MSSM maximum cross-sections are typically of the order of (0.1−1) fb, and hence some 100−1000 times smaller (similar to the fusion production of Higgs boson pairs in the SM).
Let us also mention that for type-II models (characterized by a heavier spectrum of Higgs boson masses), the leading decay modes for each of the Higgs bosons in a typical final state will be into heavy quarks. For example, if we take the channel (2.4) corresponding to the Higgs boson pair h 0 h 0 , each of the neutral Higgs bosons will mainly decay as h 0 → bb. In this region of parameter space, the alternate Higgs boson decays into gauge bosons (such as h 0 → W + W − , ZZ) are either kinematically forbidden or simply not dominant. In practice we would therefore expect to see a sizable number of 4-prong final states consisting of highly energetic (M h 0 /2 > 50 GeV ) b b-jets (or t b-jets from charged Higgs decays from the H + H − final state), which should be clearly distinguishable in a linear e + e − -collider.
We have also demonstrated that by exploring e + e − collisions at even higher energies (say, up to the characteristic 3 TeV range expected for CLIC) the opportunities could still be better. The most favored channel may then reach a cross-section at the level of 5 pb, and thus producing around half million events in the same range of integrated luminosity. In general, this upgrading should enable to perform accurate measurements of the cross-sections of several channels (2.4) and disentangle enough correlations among the parameters of the model so as to be able to pin down the corresponding Higgs boson masses and couplings at high precision. Needless to say, a truly accurate analysis demands the incorporation of quantum effects in the trilinear interactions. Such study, however, goes far beyond the scope of the present letter, whose main aim is only to show that clear signs of new physics can be highlighted already from an analysis of Higgs boson production through gauge boson fusion at leading order. Further investigation on these 2HX channels may also be of interest in the context of the LHC. However, appropriate studies of the distribution of the signal versus the background should be undertaken in order to ascertain whether the dominant signatures of the 2HX events (in the form of heavy-quark jets) could be disentangled from the important QCD background inherent in the physics of that collider.
In summary, experiments at linear colliders such as the ILC and CLIC can provide the cleanest signals of new physics and can be of paramount importance as a high precision tool to underpine the most sensitive building blocks of the gauge theories, and most significantly the structure of the Higgs potential.