Neutral Higgs-pair production at Linear Colliders within the general 2HDM: quantum effects and triple Higgs boson self-interactions

The pairwise production of neutral Higgs bosons is analyzed in the context of the future linear colliders, such as the ILC and CLIC, within the general Two-Higgs-Doublet Model (2HDM). The corresponding cross-sections are computed at the one-loop level in full compliance with the current phenomenological bounds and the stringent theoretical constraints inherent to the consistency of the model. We uncover regions across the 2HDM parameter space, mainly for low tan\beta near 1 and moderate values of the relevant lambda_5 parameter, wherein the radiative corrections to the Higgs-pair production cross sections can comfortably reach 50% This behavior can be traced back to the enhancement capabilities of the trilinear Higgs self-interactions -- a trademark feature of the 2HDM, with no counterpart in the Minimal Supersymmetric Standard Model. Interestingly enough, the quantum effects are positive for energies around 500 GeV, thereby producing a significant enhancement in the expected number of events precisely around the fiducial startup energy of the ILC. The Higgs-pair production rates can be substantial, typically amounting to a few thousand events per 500 inverse femtobarn of integrated luminosity. In contrast, the corrections are negative in the highest energy range (1 TeV). We also compare the exclusive pairwise production of Higgs bosons with the inclusive gauge boson fusion channels leading to 2H+X finals states, and also with the exclusive triple Higgs boson production. We find that these multiparticle final states can be highly complementary in the overall Higgs bosons search strategy.


I. INTRODUCTION
It was already in the 60's when the pioneering works of Higgs, Kibble, Englert and others suggested the existence of a fundamental spinless building-block of Nature [1], whose non-vanishing vacuum expectation value could explain the spontaneous breaking of the SU (2) L ⊗ U (1) Y gauge group of the Electroweak (EW) interactions down to the U (1) em abelian symmetry group of the Electromagnetism. Indeed, this so-called Higgs mechanism is the most fundamental longstanding issue that remains experimentally unsettled in Particle Physics. For one thing it is the only known strategy capable of building up a (perturbative) renormalizable quantum field theoretical description of the Electroweak Symmetry Breaking (EWSB) phenomenon. It is difficult to overemphasize that, to a great extent, the Higgs mechanism embodies the backbone of the SM structure and that in its absence we would have to cope with an entirely different conception of the inner functioning of the SM as a quantum field theory (QFT) .
The central assumption here is the existence of (at least) one elementary spinless field -the so-called Higgs boson. Its relevance is particularly evident if we take into account that the presence of one or more such fields in the structure of the model is indispensable for the unitarity of the theory. If no Higgs bosons exist below the TeV scale, weak interactions would actually become strong at that scale and e.g. the W + W − cross-section would violate the unitarity of the scattering matrix. Moreover, in the absence of Higgs bosons the various particle masses could not be generated consistently (i.e. without spoiling the ultraviolet behavior of the theory at higher orders of perturbation theory). In short, if no Higgs bosons are there to protect the (presumed) perturbative structure of the weak interactions, the latter would enter a perilous runaway regime at high energies. It is, thus, essential to confirm experimentally the existence of one or more Higgs boson particles through their explicit production in the colliders.
Amazingly enough, more than 40 years after the Higgs mechanism was naturally blended into the conventionally accepted SM landscape of the strong and electroweak interactions [2], it remains still thickly curtained and no direct evidence has been found yet of its main testable prediction: the existence of fundamental spinless particles. Notwithstanding the many efforts devoted at LEP and at the Tevatron in the last decades, all experimental searches for Higgs boson signatures have come up empty-handed, and one has to admit that no conclusive physical signature for the existence of these elementary scalar fields -not even the single one predicted by the SM -has been found for the time being. Thus, we do not know whether the conventional Higgs boson exists at all, or if there are no Higgs particles of any sort. On the contrary, it could be that there are extensions of the SM spanning a richer and elusive spectrum of Higgs bosons which have escaped all our dedicated experimental searches. Be that as it may, this pressing and highly intriguing enigma will hopefully be unveiled soon, specially with the advent of the new generation of supercolliders, like the brand new LHC, the future International Linear Collider (ILC) [3] and, hopefully, the Compact Linear Collider (CLIC) too [4].
Therefore, we have to (and we can) be well prepared to recognize all kind of hints from Higgs boson physics. Let us recall that apart from the single neutral (CP-even) Higgs boson H 0 of the SM, other opportunities arise in model building that could serve equally well the aforementioned purposes. Perhaps the most paradigmatic extension of the SM is the Minimal Supersymmetric Standard Model (MSSM) [5], whose Higgs sector involves two doublets of complex scalar fields. The physical spectrum consists of two charged states, H ± , two neutral CPeven states h 0 , H 0 (with masses conventionally chosen as M h 0 < M H 0 ) and one CP-odd state A 0 . We shall not dwell here on the whys and wherefores of supersymmetry (SUSY) as a most sought-after realization of physics beyond the SM. It suffices to say that it provides a Higgs sector which is stable (to all orders in perturbation theory) under embedding of the low-energy structure into a Grand Unified Theory, and in particular it provides the most natural link with gravity [5].
But, how to test the structure of the Higgs sector? In the MSSM case, such structure is highly constrained by the underlying supersymmetry. A consequence of it, which is particularly relevant for our considerations, is the fact that the self-interactions of the SUSY Higgs bosons turn out to be largely immaterial from the phenomenological point of view, in the sense that they cannot be enhanced (at the tree-level) as compared to the ordinary gauge interactions and, therefore, cannot provide outstanding signatures of physics beyond the SM. The bulk of the enhancing capabilities of the MSSM Lagrangian is concentrated, instead, in the rich structure of Yukawa-like couplings between Higgs bosons and quarks or between Higgs boson and squarks, or even between squarks and chargino-neutralinos, as it was shown long ago in plentiful phenomenological scenarios involving quantum effects on the gauge boson masses [6,7] and gauge boson and top quark widths [8], as well as in many other processes and observables (see e.g. [9][10][11][12][13]. For the present state of the art, cf. [14,15]; and for recent comprehensive reviews, see [16,17], for example. The two-Higgs SU L (2)-doublet structure of the Higgs sector in the MSSM is a trademark prediction of SUSY invariance. Nonetheless, the very same doublet structure can appear in the form of a non-supersymmetric framework, the so-called general (unconstrained) Two-Higgs-Doublet Model(2HDM), which does also entail a rich phenomenology [18]. Again, the same spectrum of CP-even (h 0 , H 0 ) and CP-odd (A 0 ) scalar fields arise. However, we should keep in mind that the most distinctive aspects of the general 2HDM Higgs bosons could well be located in sectors of the model quite different from the MSSM case [19]. Indeed, even the very couplings involved in the structure of the Higgs potential (viz. the so-called trilinear and quartic Higgs boson self-interactions) could be the most significant ones as far as the phenomenological implications are concerned, rather than just the enhanced Yukawa couplings with the heavy quarks. In such circumstance, we should be able to detect a very different kind of experimental signatures. Assuming, for instance, that the LHC will find evidence of a neutral Higgs boson, a critical issue will be to discern whether the newcomer is compatible with the SM or any of its extensions and, in the latter case, to which of these extensions it most likely belongs. In order to carry out this "finer" hunting of the Higgs boson(s), a TeV-range accessible linear collider machine (such as the aforesaid ILC/CLIC) will be needed. In this work, we will show how useful it could be such collider to discriminate neutral supersymmetric Higgs bosons from generic (non-SUSY) 2HDM ones.
The paper is organized as follows. In the next section, we remind the reader of some relevant Higgs boson production processes in the linear colliders. In Section III, we describe the structure of the 2HDM and the relevant phenomenological restrictions, including the properties of the trilinear Higgs self-interactions and their interplay with the notion of unitarity and vacuum stability. In Section IV, we develop in detail the renormalization program for the Higgs sector of the general 2HDM. The theoretical setup for the one-loop computation of the Higgs boson cross-sections is elaborated in Section V, leaving Section VI to present a comprehensive numerical analysis of our results. Finally, Section VII is devoted to a thorough general discussion and to present the conclusions of our work.

II. HIGGS BOSON PRODUCTION AT LINEAR COLLIDERS
The leading Higgs boson production mechanisms at the LHC have been studied thoroughly for the last twenty years, they are well under control both in the SM [20] and in the MSSM [21,22] and there are also some studies in the 2HDM [23]. Hopefully, they may soon help revealing some clues on Higgs boson physics at the LHCfor a review see e.g. [17] and references therein. However, it is not obvious that it will be possible to easily disentangle the nature of the potentially produced Higgs boson(s). The next generation of TeV-class linear colliders (based on both e + e − and γγ collisions), such as the ILC and the CLIC projects [3,4], will be of paramount importance in finally settling the experimental basis of the "Higgs issue" as the most fundamental theoretical construct of the SM of electroweak interactions. Thanks to its extremely clean environment (in contrast to the large QCD background inherent to a hadronic machine such as the LHC), a linear collider (linac for short) should allow for a precise measurement of the Higgs boson parameters, such as: i) the Higgs boson mass (or masses, if more than one); ii) the couplings of the Higgs bosons to quarks, leptons and gauge bosons; iii) the Higgs boson self-couplings mentioned above, i.e. the trilinear (3H) and quartic (4H) Higgs boson self-interactions. At the end of the day, a linac should allow us to dig deeper than ever into the structure of the EWSB and, hopefully, even to reconstruct the Higgs potential itself.
A detailed road map of predictions for Higgs-boson observables at the linear colliders is called for, with a special emphasis on those signatures which may be characteristic of the different extensions of the SM. For example, as indicated above, triple Higgs boson self-interactions (3H) may play a cardinal role in this endeavor because, in favorable circumstances, they could easily distinguish between the MSSM and the general 2HDM Higgs sectors. Such 3H couplings can mediate a plethora of interesting processes. The key point here is the potentially large enhancements that the 3H couplings may accommodate in the general 2HDM case, in contrast to the supersymmetric extensions. Actually, even for the SUSY case there is a large number of works attempting to extract vestiges of non-standard dynamics in these couplings, mainly through the radiative corrections that they can undergo. For instance, the 3H-couplings have been investigated in [24][25][26], and in some cases considered for possible phenomenological applications in TeV-class linear colliders through the double-Higgs strahlung process e + e − → HHZ or the double-Higgs WW-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 [25]. Even worse is the situation regarding the triple Higgs boson production in the MSSM, in which -except in the case of some particular resonant configuration -the typical cross-sections just border the line of ∼ 0.01 fb or less [25]. In the latter reference, for instance, it has been shown that if the double and triple Higgs production cross-sections would yield sufficiently high signal rates, the system of couplings could in principle be solved for all trilinear Higgs self-couplings up to discrete ambiguities using only these processes. But this is perhaps a bit too optimistic since in practice these cross-sections are manifestly too small to be all measurable in a comfortable way.
In stark contrast to this meager panorama within the SM and the MSSM, 3H couplings have been tested for tree-level processes in the context of the unconstrained 2HDM. For example, in Ref. [27] the tree-level production of triple Higgs boson final states was considered in the ILC. There are three classes of processes of this kind compatible with CP-conservation, namely 1) e + e − → H + H − h , 2) e + e − → hhA 0 , 3) e + e − → h 0 H 0 A 0 , (h = h 0 , H 0 , A 0 ) (1) where, in class 2), we assume that the two neutral Higgs bosons h must be the same, i.e. the allowed final states are (h hA 0 ) = (h 0 h 0 A 0 ), (H 0 H 0 A 0 ) and (A 0 A 0 A 0 ). The cross-sections in the 2HDM were shown to reach up to O(0.1) pb [27], i.e. several orders of magnitude over the corresponding MSSM predictions [28]. Besides the exclusive production of three Higgs bosons in the final state, also the inclusive pairwise production of Higgs bosons may be critically sensitive to the 3H self-interactions. Sizable production rates, again in the range of 0.1 − 1pb have recently been reported in Ref. [29], whose focus was on the inclusive Higgs boson-pair production at order O(α 4 ew ) through the mechanism of gauge-boson fusion Similar effects have also been recently computed on 2Hstrahlung processes of the guise e + e − → Z 0 h h [30]. In a complementary way, 3H couplings can also be probed in loop-induced processes. For instance, Ref. [31] presents a computation of the single Higgs boson production rate through the scattering process of i) two real photons, using the γγ mode of a linear collider, i.e. γγ → h; and ii) the more traditional mechanism of virtual photon fusion in e + e − colliders, e + e − → γ * γ * → h + X. In either case, the obtained cross sections within the general 2HDM (up to 1 pb for the first mechanism, and 0.01 pb for the latter) are 1 to 2 orders of magnitude larger than the expected SM yields and also well above the MSSM results (see [31] and references therein). Promising signatures have also been reported within the general 2HDM using loop-induced production of two neutral Higgs bosons through real γγ collisions [32], although in this case the cross-sections are smaller than in the primary mechanism γγ → h mentioned before. A crucial observation concerning our aim here is the following: all the processes described above are directly sensitive to the 3H self-couplings already at the leading order. In this paper, we continue exploiting the properties of these couplings in the general 2HDM, but we concentrate now on their impact at the level of indirect effects through radiative corrections. The latter may significantly affect processes that are kinematically more favored (e.g. because they have a smaller number of particles in the final state), but which are nevertheless totally insensitive to the trilinear Higgs boson couplings at the lowest order. Specifically, we wish to concentrate on identifying the largest quantum effects that the 3H couplings can stamp on the cross-sections for two-body neutral Higgs boson final states: Surprisingly enough, very little attention has been paid to these basic production processes and to the calculation of the corresponding radiative corrections within the 2HDM. In contrast, a lot of work has been invested on them from the point of view of the MSSM [24,25,[33][34][35][36][37][38]] -see also [39] and references therein 1 . To the best of our knowledge, only Refs. [41][42][43] have addressed this topic in the general 2HDM, although they are restricted to the production of charged Higgs pairs H + H − . In the present paper, we shall be concerned exclusively on computing the quantum effects involved on the neutral Higgs boson channels (3) at order O(α 3 ew ) (and eventually also the leading O(α 4 ew ) terms). As we will see, very large quantum effects (of order 50%) may arise on the crosssections of the two-body Higgs boson channels (3) as a result of the enhancement capabilities of the triple (and, to a lesser extent, also the quartic) Higgs boson selfinteractions. These effects are completely unmatched within the MSSM and should therefore be highly characteristic of non-supersymmetric Higgs boson physics.

III. THE TWO-HIGGS-DOUBLET MODEL: GENERAL SETTINGS AND RELEVANT RESTRICTIONS
The Two-Higgs-Doublet Model (2HDM) is defined upon the canonical extension of the SM Higgs sector with a second SU (2) L doublet with weak hypercharge Y = 1, so that it contains 4 complex scalar fields arranged as follows: (4) With the help of these two weak-isospin doublets, we can write down the most general structure of the Higgs potential fulfilling the conditions of CP-conservation, gauge invariance and renormalizability. Moreover, a discrete Z 2 symmetry Φ i → (−1) i Φ i (i = 1, 2) -which will be exact up to soft-breaking terms of dimension 2 -is usually imposed as a sufficient condition to guarantee the proper suppression of the Flavor-Changing Neutral-Current (FCNC) effects that would otherwise arise within the quark Yukawa sector [44] 2 (we shall return to this (A 0 ); and a couple of charged Higgs bosons (H ± ). The free parameters of the general 2HDM are usually chosen to be: the masses of the physical Higgs particles; the aforementioned parameter tan β; the mixing angle α between the two CP-even states; and finally the coupling λ 5 , which cannot be absorbed in the previous quantities and becomes tied to the structure of the Higgs self-couplings. To summarize, the vector of free inputs reads Therefore, we are left with 7 free parameters, which indeed correspond to the original 6 couplings λ i and the two VEV's v 1 , v 2 -the latter being submitted to the constraint v 2 ≡ v 2 1 + v 2 2 = 4 M 2 W /g 2 , which is valid order by order in perturbation theory, where M W is the W ± mass and g is the SU L (2) gauge coupling constant. The dimension 2 term that softly breaks the Z 2 symmetry can be written as −m 2 12 Φ † 1 Φ 2 + h.c., with the second equality being valid at the tree-level. Given tan β one may trade the parameter λ 5 for m 12 through this formula, if desired. Since we are going to adopt the on-shell renormalization scheme [46][47][48] for the one-loop calculation of the cross-sections (cf. section IV), it is convenient to introduce the electromagnetic fine structure constant α em = e 2 /4π, which is one of the fundamental inputs in this scheme (together with the physical Z 0 mass, M Z , and the Higgs and fermion particle masses). The electron charge −e is connected to the original SU L (2) gauge coupling and weak mixing angle θ W through the well-known relation e = g sin θ W , which is also preserved order by order in perturbation theory. Then the W ± mass can be related to α em and the Fermi constant G F in the standard way in the on-shell scheme, namely where the parameter ∆r [48] vanishes at the tree-level, but is affected by the radiative corrections to µ-decay both from standard as well as from new physics, for instance from the MSSM [8,14] and also from the 2HDM (see below). Notice that, in the above formula, M W also enters implicitly through the relation sin 2 θ W = 1 − M 2 W /M 2 Z , which is valid order by order in the onshell scheme. Since G F can be accurately determined from the µ-decay, and M Z has been measured with high precision at LEP, it is natural to take them both as experimental inputs. Then, with the help of the non-trivial relation (11), the W ± mass can be accurately predicted in a modified on-shell scheme where G F , rather than M W , enters as a physical input.
We shall not dwell here on the detailed structure of the Yukawa couplings of the Higgs boson to fermions (see [18] for an exhaustive treatment), since their contribution to the processes e + e − → h 0 A 0 /H 0 A 0 under study at the one-loop level is largely subdominant in the main regions of parameter space. However, it is precisely the coupling pattern to fermions that motivates the distinction between the different types of 2HDM's, so let us briefly describe them. Assuming natural flavor conservation [44], we must impose that at most one Higgs doublet can couple to any particular fermion type. As a result, the coupling of the 2HDM Higgs bosons with fermions (say the quarks) can be implemented in essentially two manners which insure the absence of potentially dangerous (tree-level) FCNC's, to wit (cf . Table III) 4 : i) type-I 2HDM, wherein the Higgs doublet (Φ 2 ) couples to all of quarks, whereas the other one (Φ 1 ) does not couple to them at all; or ii) type-II 2HDM, in which the doublet Φ 1 (resp. Φ 2 ) couples only to down-like (resp. up-like) right-handed quarks. In the latter case, an additional discrete symmetry involving the chiral components of the fermion sector must be imposed in order to banish the tree-level FCNC processes, e.g.
for the down and up right-handed quarks in the three flavor families (i = 1, 2, 3). In turn, the Higgs self-couplings λ i in the Higgs potential can be rewritten in terms of tan β and the physical parameters of the on-shell scheme such as the masses and the electromagnetic fine structure constant α em . At the tree-level, we have where s W = sin θ W , c W = cos θ W . From these Lagrangian couplings in the Higgs potential, we can de-type I type II h 0 tt cos α/ sin β cos α/ sin β h 0 bb cos α/ sin β − sin α/ cos β H 0 tt sin α/ sin β sin α/ sin β H 0 bb sin α/ sin β cos α/ cos β A 0 tt cot β cot β A 0 bb − cot β tan β rive the "physical couplings", namely those affecting the physical Higgs bosons in the mass-eigenstate basis. We call these couplings the triple (3H) and quartic (4H) Higgs couplings. Their behavior and enhancement capabilities are at the very core of our discussion. Indeed, from our analysis it will become clear that they furnish the dominant source of quantum corrections to the Higgspair production processes (3) within the general 2HDM. The physical 3H and 4H couplings are not explicitly present in the 2HDM potential (5). They are derived from it after spontaneous breaking of the EW symmetry and corresponding diagonalization of the Higgs boson mass matrix using the rotation angles α and β in (7). In the particular case of the SM, the trilinear and quartic Higgs couplings have fixed values, which depend uniquely on the actual mass of the Higgs boson. In the MSSM, however, and due to the SUSY invariance, the Higgs boson self-couplings are of pure gauge nature, as we shall revise briefly below. This is in fact the primary reason for the tiny production rates obtained for the triple Higgs boson processes (1) within the framework of the MSSM [25]. In contrast, the general 2HDM accommodates trilinear Higgs couplings with great potential enhancement. The full list is displayed in Table II. As can be seen, the couplings in this table depend on the 7 free parameters (9). In the particular case where λ 5 = λ 6 , this table reduces to Table 1 of Ref. [27]. The equality of these couplings takes place automatically e.g. in the MSSM, where in addition other simplifications occur, as we discuss below. Let us note, for example, that in the limit α = β − π/2 the h 0 h 0 h 0 -trilinear coupling in Table II  Let us recall that the MSSM Higgs sector is a type-II one of a very restricted sort: it is enforced to be SUSY invariant. This is a very demanding requirement as, for example, it requires that the superpotential has to be a holomorphic function of the chiral superfields [5,18]). This implies that we cannot construct the MSSM with two Y = +1 Higgs superfields. Then, in order to be able to generate masses for both the top and bottom quarks through EWSB, the Φ 2 doublet can be kept as it is with Y = +1, although we must replace Φ 1 with the conjugate (Y = −1) SU L (2) doublet where ǫ = i σ 2 , with σ 2 the second Pauli matrix. Thus, for the first doublet the correspondence with the MSSM case reads Φ 1 = −ǫ H * 1 , and the second doublet (the one with no change) is just relabeled H 2 . Moreover, the treelevel relationships imposed by Supersymmetry between the λ i couplings of the potential (5) are the following: Substituting (13) and (14) in (5) one obtains the usual MSSM potential, which in practice must be supplemented with soft SUSY-breaking scalar mass terms m 2 i H 2 i , including a bilinear mixing term m 2 12 H 1 H 2 for the two Higgs doublets. This term is the analog of (10) for the 2HDM, but in the soft-SUSY breaking context is arbitrary. The result can be cast as follows: where µ is the higgsino mass term in the superpotential [5]. The obtained potential is one where all quartic couplings become proportional to α em , which is tantamount to say that in the SUSY case these self-couplings are purely gauge. Thus, after EWSB also the trilinear self-couplings will be purely gauge. As warned before, this is the primary reason for their phenomenological inconspicuousness. Under the supersymmetric constraints, the entries of Table II boil down to the MSSM form listed e.g. in [18]. Moreover, the five constraints (14) reduce the number of free parameters from 7 to 2 in the supersymmetric context, typically chosen to be tan β and M A 0 . For the general 2HDM case, however, we shall stick to the form presented in Table II, where the 7 free inputs are chosen as in (9).
In another vein, it is of paramount importance when studying possible sources of unconventional physics to make sure that the SM behavior is sufficiently well reproduced up to the energies explored so far. Such a requirement translates into a number of constraints over the parameter space of the given model. In particular, these constraints severely limit the a priori enhancement possibilities of the Higgs boson self-interactions in the 2HDM.
• To begin with, we obviously need to keep track of the exclusion bounds from direct searches at LEP and the Tevatron. These amount to M l 114GeV for a SM-like Higgs boson -although the bound is weakened down to 92.8 GeV for models with more than one Higgs doublet [51]. In the 2HDM all the mass bounds can be easily satisfied upon a proper choice of the different Higgs masses -which, unlike the SUSY case, are not related with each other.
• Second, an important requirement to be enforced is related to the (approximate) SU L+R (2) custodial symmetry satisfied by models with an arbitrary number of Higgs doublets [52]. In practice, this restriction is implemented in terms of the parameter ρ, which defines the ratio of the neutralto-charged current Fermi constants. In general it takes the form ρ = ρ 0 + δρ, where ρ 0 is the tree-level value. In any model containing an arbitrary number of doublets (in particular the 2HDM), we have ρ 0 = M 2 W /M 2 Z c 2 W = 1, and then δρ represents the deviations from 1 induced by pure quantum corrections. From the known SM contribution and the experimental constraints [51] we must demand that the additional quantum effects coming from 2HDM dynamics satisfy the approximate condition |δρ 2HDM | 10 −3 . It is thus important to stay in a region of parameter space where this bound is respected. In our calculation we include the dominant part of these corrections, which involves oneloop contributions to δρ mediated by the Higgsboson and yields [53] From this expression, it is clear that arbitrary mass splittings between the Higgs bosons could easily overshoot the limits on δρ. However, we note that if M A 0 → M H ± then δρ 2HDM → 0, and hence if the mass splitting between M A 0 and M H ± is not significant δρ 2HDM can be kept under control 5 . Let us also recall in passing that the δρ correction translates into a contribution to the parameter ∆r of (11) given by −(c 2 W /s 2 W ) δρ. Therefore, the tight bounds on ∆r, together with the direct bounds on the ratio of charged and neutral weak interactions, restrict this contribution as indicated above.
• Also remarkable are the restrictions over the charged Higgs masses coming from FCNC radiative B-meson decays, whose branching ratio B(b → sγ) ≃ 3 × 10 −4 [51] is measured with sufficient precision to be sensitive to new physics. Current lower bounds render M H ± 295 GeV for tan β 1 [55,56]. It must be recalled that these bounds apply for both Type-I and Type-II models in the region of tan β < 1. In contrast, they are only relevant for Type-I in the large tan β regime, since for them the charged Higgs couplings to fermions are proportional to cot β and hence the loop contributions are highly suppressed at large tan β -cf.
Ref. [57] for detailed analytical expressions of the contributions to B(b → sγ) from the 2HDM. Besides, data from B → lν l may supply additional restrictions on the charged Higgs mass at large tan β for the type-II 2HDM [58]. For a very recent updated analysis of the different B-physics constraints on the 2HDM parameter space, see Ref. [45].
• Further constraints apply to tan β coming from the following two sources: i) The Z 0 → bb and B − B mixing processes strongly disfavor tan β below 1 [56]; ii) The requirement that the Higgs couplings to heavy quarks remain perturbative translates into an (approximate) allowed range of O(0.1) < tan β < 60 [49].
• Besides the available experimental data, an additional set of requirements ensues from the theoretical consistency of the model. In particular, we shall introduce the following set of conditions ensuring the stability of the vacuum [59][60][61]: The unitarity bounds deserve a separate discussion. Indeed, they turn out to be a fundamental ingredient of our computation. As we have already been mentioning, the trademark behavior that we aim to explore shall depend critically on the enhancement capabilities of the Higgsboson self-interactions which, in turn, will be sharply constrained by unitarity. The basic idea here is that, within perturbative QFT, the scattering amplitudes are "asymptotically flat", meaning that they cannot grow indefinitely with the energy. This is tantamount to say that the unitarity of the S-matrix must be guaranteed at the perturbative level. In a pioneering work by Lee, Quigg and Thacker (LQT) [62], the above arguments were first applied to the SM. As a very relevant outcome, an upper bound for the Higgs boson was obtained: The basic idea beneath the LQT argument is to compute the scattering amplitudes in a variety of processes (namely scalar-scalar, scalar-vector and vector-vector scattering; most particularly, the longitudinal components of the vector bosons) and demand all of them to be in accordance with the generic tree-level unitarity condition. The latter is a pure quantum-mechanical requirement, which can be expressed in terms of the s-wave scattering lenght in the high-energy limit: Taking into account the equivalence between the longitudinal components of the gauge and the Goldstone bosons (which is exact to order O( M 2 W /s), s being the centerof-mass energy squared), we can convince ourselves that only the scalar-scalar processes will be relevant to this concern. Moreover, it can also be shown that the bulk of the contribution is brought about by the quartic scalar self-couplings -the triple ones are suppressed as 1/s. Thereby the quartic Higgs self-couplings (4H) turn out to be the most strongly constrained by the bounds of Eq. (19).
Several authors have exported the above ideas to the 2HDM [63][64][65][66][67][68][69][70][71][72][73][74][75][76] and also to its complex (CP-violating) extensions [77]. The strategy here is to compute the Smatrix elements S ij in a set of 2 → 2 scattering processes involving Higgs and Goldstone bosons. In this case there are more scalar-scalar channels than in the SM, which are coupled among themselves. At this point, we could perform a unitary transformation U relating the S-matrix expressed in the original weak-interaction eigenstate basis (S w ) with the S-matrix written in the physical mass-eigenstate basis (S m ). Since, however, the former is much simpler than the latter and both carry exactly the same information, we may generalize the original LQT procedure [62] by focusing on S w and requiring that its eigenvalues satisfy the tree-level unitarity condition |α i | < 1 (∀i). In practice, inspired by the scattering length analysis (19), we require |α i | < 1/2 (∀i). This leads to a set of conditions over several linear combinations of the quartic couplings in the original Higgs potential [68][69][70][71][72][73][74][75][76], as can be seen from the explicit form of the various eigenvalues α i at the tree-level: (20)

IV. RENORMALIZATION OF THE 2HDM HIGGS SECTOR
In this Section, we discuss in detail the renormalization of the Higgs sector. The renormalization of the SM fields and parameters is performed in the conventional on-shell scheme in the Feynman gauge, see e.g. Ref. [46-48, 78, 79]. At present, highly automatized procedures are available for loop calculations, especially in the MSSM, see e.g. [16] and [80][81][82][83]. However, in our case the calculation of the cross-sections for the processes (3) is performed within the general (non-supersymmetric) 2HDM and we must deal with the renormalization of the Higgs sector in this class of generic models. To this end, we attach a multiplicative wave-function (WF) renormalization constant to each of the SU L (2) Higgs doublets in the 2HDM, The renormalized fields are those on the r.h.s. of these expressions. At one-loop we decompose Z Φ1,2 = 1 + δZ Φ1,2 + O(α ew ). These WF renormalization constants in the weak-eigenstate basis can be used to construct the WF renormalization constants Z hihj = 1 + δZ hihj in the mass-eigenstate basis by means of the set of linear transformations (7). We are interested only in the production of neutral Higgs bosons in this paper, and therefore we shall only quote the relations referring to them. The corresponding δ Z hihj in the neutral Higgs sector (using the notation δ Z hihi ≡ δ Z hi ) are determined as follows: where in the last two equations we have also included the mixing of the two CP-even Higgs bosons and that of the CP-odd state with the neutral Goldstone boson in the Feynman gauge. The different counterterms δ Z hihj must be anchored by specifying a set of subtraction conditions on a finite number of 2-point functions. To that purpose, we need to specify the renormalized self-energies corresponding to the gauge bosons and each of the physical Higgs-boson fields in the model. First of all let us define the relation between the polarization tensor of a gauge boson and the transverse self-energies. We define it as follows, where the dots indicate that we neglect the longitudinal part (proportional to q µ q ν ) because our calculations for a linear collider involve extremely light external fermions (electrons and positrons) and as a result the longitudinal contributions are suppressed as m 2 e /M 2 V . The same applies to the q µ q ν terms in the transverse part, of course. We identify +iΣ µν (q) with the Feynman diagram associated to the structure of the polarization tensor. As a result, the unrenormalized self-energy part is defined in such a way that −i Σ V V ′ (q 2 ) (notice the minus sign) is equated to the blob of external lines V and V ′ . For simplicity, we define Σ V (q 2 ) ≡ Σ V V (q 2 ) and similarly for other fields. Thus e.g. the free propagator for the gauge boson V becomes "dressed" by the quantum effects as follows: −i g µν being the bare mass of V and where in the second equality we have introduced the renormalized form of the propagator. Similarly, for self-energies involving two scalars (cf. Fig. 1) we identify +i Σ φiφj (q 2 ) with the Feynman blob with external legs φ i and φ j . For instance, the dressed propagator of the scalar φ with bare mass m (0) φ is given by i .
The last equality refers to the corresponding renormalized scalar field propagator, m φ being the renormalized mass andΣ φ the renormalized self-energy. Without wanting to appear too pedagogical, we remind the reader that the renormalization transformation at one-loop is performed by canonically introducing renormalization constants relating every bare field φ (0) i to the renormalized ones as follows φ this case takes on the particular form (21)and at the same time decomposing the bare masses and couplings as the sum of the renormalized quantity plus a counterterm (m With these definitions, we can check that the relations (24) and (25) are consistent at one-loop if we define the renormalized self-energy of V and φ aŝ (26) Finally, for blobs involving a vector (V ) and a scalar (φ) field (cf. Fig. 1), there is a free vector index µ and , where ±q is the external momentum of φ flowing into (out of) the blob. As with the tensor iΣ µν (q), we equate the vector Σ µ φ V (q) (in this case, with no i-factor) to the corresponding blob with external lines φ and V 6 . With these conventions, we may write the various scalar-scalar self-energies needed for this calculation involving both Higgs bosons and Goldstone bosons in diagonal or mixing form. For the CP-even sector: For the CP-odd sector: The structure of the Higgs sector beyond the leading order is somewhat more involved. Similarly as in the SM, scalar-vector mixing terms arise in the Lagrangian of the general 2HDM. Such contributions originate from the gauged kinetic terms of the Higgs doublets, (29) where, in our conventions, the SU (2) L × U (1) Y covariant derivative reads (using pretty standard notations), with −e being the electron charge. Upon spontaneous EW symmetry breaking, the following scalar-vector mixing terms are generated: The subindex in L SV 0 means that this is a bare Lagrangian, and therefore all fields in it are bare fields φ i ) -although we did not bother to label them as such in order to avoid an increasing notational snarl-up.
Similarly for the VEV's v i of the two Higgs doublets. Here we have to include the corresponding WF renormalization constants from (21). Therefore, the decomposition of the bare VEV's into the renormalized ones plus counterterms reads as follows: As we shall explain below in more detail, we will adopt the following renormalization condition for the parameter tan β = v 2 /v 1 [33,84]: This condition insures that the ratio v 2 /v 1 is always expressed in terms of the true vacua (that is, the VEV's resulting from carrying out the renormalization of the Higgs potential). The corresponding counterterm tan β → tan β+δ tan β associated to this definition therefore reads Proceeding in this way, we may rewrite Eq. (31) as L SV and where all the fields and parameters here are the renormalized ones. Similarly, the overall counterterm Lagrangian reads where we have reabsorbed the parameter v = v 2 1 + v 2 2 into the gauge boson masses, M W = gv/2, M Z = M W /c W and their associated counterterms.
The specific structure of the various mass and WF counterterms becomes fixed after we choose the appropriate renormalization conditions in a given subtraction scheme. In particular, a set of prescriptions for the determination of the Higgs mass counterterms arising in Eq. (28) is called for. Following the on-shell prescription, we arrange that the four renormalized Higgs boson masses coincide with the four input (physical) Higgs masses, and hence we enforce them to be the pole masses of the corresponding renormalized propagators (25). Equivalently, we declare that the real part of the corresponding renormalized self-energies vanishes for onshell Higgs bosons. In a nutshell: Incidentally, let us point out that this is a different approach with respect to the MSSM (see e.g. [16]), in which case and owing to the supersymmetric invariance, less freedom is left over to settle independent renormalization conditions. However, additional conditions are needed due to the extended structure of the Higgs sector of the 2HDM with respect to the simpler SM case. In particular, the condition (33) has been applied, but we still have to fix the structure of the counterterm δ tan β in (36), which we will unravel later on below.
The scalar-vector mixing contributions at the tree-level in (36) are canceled by the addition of gauge-fixing terms. Using the 't Hooft-Feynman gauge choice ξ = 1 in the class of linear R ξ gauges, we have: Let us, however, point out that the scalar-vector mixing can always take place via quantum corrections -see the diagrams of Fig. 1, and thereby a related set of counterterms must be introduced in order to renormalize the corresponding mixed self-energies. In particular, from Eq. (36) it follows that the A 0 − Z 0 mixing counterterm reads where in the last equality we have taken δ Z A 0 G 0 from Eq. (22) and used (34). Let us also derive the explicit expression for the renormalized A 0 − Z 0 mixing self-energy, as displayed in Fig. 1a. By plugging the corresponding Renormalized scalar-vector and scalar-scalar selfenergies at all orders in perturbation theory.
Feynman rules, we may cast conveniently the amplitude in the form where the structure of the renormalized A 0 − Z 0 selfenergy encapsulates both the corresponding bare selfenergy contribution and the related counterterm as follows:Σ with δ Z A 0 Z 0 given by Eq. (39). In the same vein, a similar expression can be derived for the mixing with the Goldstone boson, G 0 : withΣ A 0 G 0 defined in (28). We are now ready to tackle the renormalization of the Higgs-boson fields. Two independent conditions are needed to fix δ Z Φ1 , δ Z Φ2 . Insofar as we will be dealing with the neutral Higgs sector, it is more convenient for us to impose such conditions on the A 0 boson: where the short-hand notation f ′ (q 2 ) ≡ d f (q 2 )/dq 2 has been employed. The above conditions, together with (37), ensure that the outgoing A 0 particle is on its mass shell and that the residue of the renormalized propagator (25) is precisely equal to 1/[1 +Σ ′ A 0 (M 2 A 0 )] = 1 -see section V. Besides, the second condition guarantees that no mixing between the A 0 and Z 0 bosons will take place at any order in perturbation theory as long as that A 0 Higgs boson is on-shell. Actually, also the scalar-scalar mixing A 0 − G 0 is generated through a similar class of one-loop interactions (recall Fig. 1b). However, one can see that both mixing phenomena turn out to be related by the Slavnov-Taylor (ST) identity which ultimately stems from the underlying BRS symmetry of the theory [85]. Notice that Eqs. (44) and (45) imply that the real part of the A 0 − G 0 mixed selfenergy, too, vanishes on-shell. Let us recall that a BRS transformation is a local Grassmann gauge transformation in which the infinitesimal gauge parameters read δω a (x) = η a (x)δλ, with η a the FP ghost fields (scalar fields with Fermi statistics) and where δλ is an infinitesimal Grassmann (i.e. anticommuting) parameter. In fact, using the notation X ≡ 0| T X |0 for a Green's function constructed out of a product of field operators X on which the time-ordering operator T acts, the BRS invariance is expressed as δ λ X = 0, where δ λ = δλ δ/δλ. Recalling the expression for the gauge-fixing Lagrangian Eq. (38), we consider the BRS invariance of the 2-point Green's function constructed with an anti-FP ghost field η 0 and a neutral CP-odd Higgs boson. We have Here we have used the standard BRS transformation of the anti-FP ghost field, which is simply related to the gauge fixing term F 0 as follows: δ λη 0 = −F 0 δλ. Moreover, the second term of the BRS variation vanishes because δ λ A ∝ A (τ a /2)η a (x)δλ and we cannot have Green's functions with external ghost legs.
Fourier transforming the expression (46), we can rewrite it in momentum space, with ∆ XY µ denoting the renormalized (non-amputated) 2-point Green's functions introduced in Eqs. (40)- (42). The ST-identity (45) is just a transcription of (47) in terms of the renormalized mixed self-energies. In words, the identity states that an on-shell A 0 does not mix neither with the neutral gauge boson Z 0 nor with the neutral Goldstone boson G 0 , provided that the conditions in Eqs. (43) and (44) hold. A similar ST identity can be derived in the charged Higgs sector (see Ref. [9]), and although we do not need it in the context of our calculation we quote it for completeness. In the present notation, it reads For an extended discussion of these identities, see Ref. [86], and in particular for the off-shell regime and its form in generalized non-linear gauges [83].
A comment on the physical definition of tan β is now in order. Recall that we have assumed that the relation (33) holds between the counterterms of the VEV's. Such relation says that the renormalized value of tan β is defined to be the ratio of the renormalized values of v 2 and v 1 . This is tantamount to enforce the Higgs tadpoles to vanish [33,84], Indeed, the tadpoles are absent at the tree-level and must be absent too at one-loop if v 1 and v 2 are to be the renormalized VEV's characterizing the true ground state of the potential at the given order of perturbation theory. A note of caution may be pertinent at this point: renormalized quantities are not always physical quantities, although sometimes they are; for example, on-shell renormalized masses are of course physical masses, and this is indeed what is stipulated by the fundamental on-shell renormalization conditions (37). However, the renormalized VEV's v i are not directly measurable quantities. So, even though the above definition of tan β may look quite natural, it is not entirely physical yet, in the sense that it is not tied to a physical observable at this point. In order to link tan β to a genuine physical quantity we have to specify the process from which tan β could be extracted from the experiment. For example, in [9] it is derived from the charged Higgs decay H ± → τ ± ν τ , whereas in other cases it is collected from other observables, e.g. [83]. The ultimate connection of the current definition of tan β with physics demands to take account of the contribution from the process-dependent quantum effects, as explained in detail e.g. in [9]. Notice that the fact that this reference focuses on the MSSM is not relevant for this feature. All that said, let us stress that these process-dependent effects are not crucially important for the present analysis because they do not involve the triple Higgs boson couplings under study. Therefore, since we wish to place ourselves in the region of parameter space where these couplings are most favored, the process-dependent corrections cannot noticeably alter the bulk supply of the significant renormalization effects on the cross-sections (3) computed under these conditions. Thus, it is not crucial for this study of the leading triple Higgs boson self-interactions to deal with the details of the particular process from which one could eventually extract the value of tan β, and for this reason we will not include these corrections in our analysis 7 . In short, in all practi- 7 The situation is, in the quantitative aspect, quite similar to the analysis of the branching ratio of the charged Higgs decay of the top quark (t → H + b) both in the 2HDM [87] and in the MSSM ("a window to virtual SUSY"-see Ref. [9]). The potential quantum effects related to the intrinsic dynamics of this decay can be so large that the process-dependent corrections associated to the definition of tan β pale in comparison and are, therefore, unessential to demonstrate the very origin of the quantum corrections.
cal respects the definition of tan β employed here, which is based on the condition (33), is well-suited to illustrate the leading renormalization effects produced by the Higgs boson self-interactions on the 2HDM cross-sections (3). Combining the renormalization conditions (43) and (44), and using (22), (28) and (39), we are lead to The WF renormalization constants for the Higgs doublets can now be solved explicitly as follows: From (34) and (51) the specific form of the counterterm associated to our definition of tan β finally ensues: Up to now, we have established a set of conditions that allows us to renormalize five out of the seven free parameters (9) of the 2HDM. These are the conditions (33), (37) and (44). Notice that the additional one (43) is a field renormalization condition and, therefore, it does not affect the definition of the parameters themselvesalthough it certainly affects the renormalization of the Green's functions and cross-sections. The two renormalized parameters that remain to be defined are the CPeven mixing angle α and the coupling λ 5 . As far as α is concerned, we can split its bare value in the usual form α (0) = α + δα, where we define the renormalized value α by setting to zero the renormalized h 0 − H 0 mixing self energy at the scale of the Higgs mass in the final state. If e.g. the final state is the lightest CP-even Higgs boson (as in Figs. 3-11), then whereas if it is the heavy CP-even Higgs boson then M 2 In what follows we presume that the final state is the lightest CP-even state. It is not difficult to see that the renormalization condition (53) fixes the mixing angle counterterm δα in terms of the mixing mass counterterm δm 2 h 0 H 0 in the expression of the renormalized mixed self-energyΣ h 0 H 0 -see Eq. (27). Obviously, (53) implies with δ Z h 0 H 0 = sin 2α(δ tan β/ tan β). On the other hand, one may show from the renormalized Higgs potential that Actually, the condition (53) is bound to the fact that the external Higgs bosons (h 0 , H 0 ) are to be on-shell to all orders in perturbation theory. To better assess the significance of our renormalization conditions in the CP-even Higgs sector, we note that the quantum corrections in this sector induce the following non-trivial inverse propagator matrix: The physical Higgs boson masses 8 ensue from the real part of the roots of det ∆ −1 h 0 H 0 = 0, i.e. of the algebraic equation The above equation makes clear that both conditions (37) and (53) are needed in order to settle the correct on-shell properties for the h 0 boson. Similarly, we have Eq. (37) and The set of renormalization conditions that we have introduced so far enable us to anchor the full collection of counterterms within the Higgs sector, in particular the mass counterterms δ m 2 G 0 , δm 2 A 0 G 0 that remain unsettled in Eq. (28). In terms of the tadpole counterterms and δ tan β (all of them already fixed in this renormalization scheme), we may write their explicit expressions as follows: The only parameter left is λ 5 . However, this parameter is not involved at the tree-level in any of the cross-sections under study, which are entirely dependent on just the electroweak gauge boson couplings at the lowest order. Therefore, there cannot be any one-loop UV divergent quantity that needs to be absorbed into the renormalization of the parameter λ 5 . In other words, this parameter is not renormalized at one-loop for the processes under consideration. It goes without saying that the situation would be different at higher orders, e.g. at 2-loop order and beyond. But we need not go that far to illustrate the main message of this paper and its potential implications. In Section VII, we shall briefly return to this point. 8 Let us recall that all Higgs boson masses are input parameters within the general 2HDM, unlike the MSSM case.

V. NEUTRAL HIGGS-BOSON PAIR PRODUCTION AT 1 LOOP IN THE LINEAR COLLIDER: THEORETICAL SETUP
Due to its vantage point as a high precision experimental machine, and hence owing to its complementarity with the LHC, it is of cardinal importance to understand in detail the phenomenology of the Higgs sector in the context of a linear collider (both within the SM and its renormalizable extensions). Quite an effort has already been devoted to this goal in the literature, specially within the MSSM, see section II and references therein. Analyses of Higgs boson production in the LHC within the 2HDM have been considered e.g. in [23,88], and the results are as encouraging as the corresponding calculation for the linear colliders [26,27,29] except that in the latter case the experimental environment is much cleaner and the real possibilities to discriminate the nature of the Higgs bosons are correspondingly higher.
Exclusive processes in the linear collider can be specially significant. The simplest neutral Higgs boson processes of this kind in a linac involve the production of two Higgs bosons in the final state, Eq. (3). At the treelevel, they are of order O(α 2 ew ) and are mediated by Z 0 exchange. The next-to-simplest ones are the triple Higgs boson channels (1). These have been studied in [27], and in section VII we will compare their cross-sections with the 2H ones under the same set of conditions.
Notice that due to CP-conservation -or, what is more, due to Bose-Einstein symmetry -processes with only two identical neutral Higgs bosons in the final state cannot proceed at the tree-level (neither in the SM nor in any of its extensions). Specifically, the combinations 2H = h 0 h 0 , H 0 H 0 , h 0 H 0 of CP-even Higgs bosons can only be generated through CP-conserving box-diagrams at O(α 4 ew ), see Fig. 2. Such diagrams render up to σ 2H ∼ 5 × 10 −6 pb at √ s = 1.5 TeV [27], which would translate into a rate of scarcely a few events per 500 fb −1 of integrated luminosity, thus a too handicapped rate to provide feasible detection signals. It is worth pointing out that, since all couplings appearing in the box-type diagrams of Fig. 2 are fixed by the gauge symmetry, the obtained results do not depend on whether we consider the SM, the MSSM or the general (unconstrained) version of the 2HDM. Being these couplings generated from the electroweak kinetic terms (29), they are fully determined by the gauge symmetry, and as a consequence they are formally the same both in the MSSM and in the 2HDM. It means that, if we aim at distinguishing between the two basic types (supersymmetric and non-supersymmetric) of Higgs boson scenarios, the study of the quantum corrections is mandatory here. Even at one-loop order, this implies to cope with a formidable calculation.
Even though the detection of the 2H processes (3) would signify an unmistakable sign of physics beyond the SM, none of them is sensible to the trilinear Higgs boson couplings at the leading order. It means that the order of magnitude of the cross-sections (1 − 30 fb) is not significantly enhanced with respect to, say, the MSSM ones. What can be remarkably affected in the general 2HDM, however, is the size of the radiative corrections, i.e. the virtual quantum effects on these processes. For the triple Higgs boson channels (1) and the gauge boson fusion ones (2), instead, already the tree-level cross-sections can be greatly modified by the triple Higgs boson selfinteractions. Indeed, one can achieve sizable production rates up to roughly 1pb, whereas the maximum MSSM payback lies around 10 −3 − 10 −2 pb (see e.g. Tables 5 and 2 of Refs. [27] and [29], respectively). These results can be traced back to the intrinsically different nature of the 3H self-coupling in the supersymmetric model (a pure gauge coupling) as compared to the non-supersymmetric ones.
The use of the two-body processes (3) as a tool to probe the nature of Higgs bosons in the final states is nevertheless possible through the analysis of the radiative corrections. Closely related to them are the two-body final states e + e − → Z 0 h (with h = h 0 , H 0 ), which have long been known to be complementary to the former in the MSSM [28]. We will not consider the latter in the present analysis [89].
The theoretical setup for the one-loop calculation starts from the basic bare interaction vertices extracted from Eq. (29): In Fourier space, they read: p and p ′ being the 4-momenta of the CP-odd and even Higgs boson pointing outwards, respectively. Indeed, one-loop quantum corrections to the e + e − → A 0 h 0 amplitude are driven by a large set of Feynman diagrams involving the following subsets of contributions: i) the self-energy corrections to the Z 0 propagator and the γ − Z 0 mixing propagator (cf. Figs. 3-4); ii) the vertex corrections to the A 0 h 0 Z 0 interaction (cf. Figs. 5-6); iii) the loop-induced γ h 0 A 0 interaction (cf. Fig. 7); iv) the vertex corrections to the e + e − Z 0 interaction (cf. Fig. 8); v) the box-type diagrams (cf. Fig. 9); vi) the finite wave-function renormalization of the external Higgsboson legs (cf. Fig. 10); and, finally, the counterterm diagrams (cf. Fig. 11). Of course an equivalent set of diagrams would be needed to describe the complementary process e + e − → A 0 H 0 , but we refrain from displaying them explicitly.
One might expect that the cross-sections under study could be sensitive to the model differences between type-I and type-II 2HDM's through the Yukawa coupling corrections (cf. e.g. the first two rows of diagrams in Fig. 5). However, the leading effects are concentrated on the triple Higgs boson self-interactions, and being the latter identical for type-I and type-II, the model differences are tiny in the relevant regions of the parameter space.
Let us also mention that for our current purposes we can safely discard the pure QED corrections to the e + e − Z 0 vertex. These QED corrections and the pure  weak ones factorize into two subsets which are separately UV finite and gauge-invariant. Moreover, these photonic contributions are fully insensitive, at the order under consideration, to the relevant 3H self-couplings on which we are focusing in this work. For the processes under consideration, where the Higgs bosons in the final state are electrically neutral, the one-loop QED effects are confined to the initial e + e − vertex. In practice, the net outcome of the accompanying initial-state real photon bremsstrahlung is to lower the effective center-of-mass energy, and hence to modify the shape of the cross sections as a function of √ s. These QED effects are the only source of infrared (IR) divergences in our calculation, and since we do not consider them the obtained scattering amplitude is IR finite. Finally, let us also emphasize that its inclusion would not change significantly the overall size of the radiative corrections, as in the regions of the parameter space where the 3H self-couplings are dominant the QED effects are comparatively negligible (as we have checked explicitly). In short, they are unessential at this stage to test the presence of the new dynamical features triggered by the 2HDM in the processes under consideration.
The UV divergent contributions arising from the set of Feynman diagrams that we have just described can be absorbed by means of the associated counterterm diagrams of Fig. 11. In particular, we need to include corresponding counterterms to renormalize the A 0 h 0 Z 0 and A 0 H 0 A 0 basic interaction vertices at the one-loop level. The analytical expressions ensue from Eqs. (60-61) after splitting parameters and fields in the usual way into renormalized ones plus counterterms. The final result reads as follows: When writing down the above counterterms, one has to keep track explicitly of the UV divergences stemming from the Higgs-Higgs (h 0 − H 0 ), Higgs-Goldstone (A 0 −G 0 ) and Higgs-vector (A 0 −Z 0 ) mixing effects which arise at the quantum level. By the same token, and despite there are no γh 0 A 0 / γH 0 A 0 tree-level couplings, we must also introduce appropriate counterterms in order to dispose of the associated UV-divergences that appear on account of the γ − Z 0 mixing at the quantum level. The relevant terms are and Likewise, we must introduce the corresponding counterterms for the e + e − Z 0 vertex, as well as for the Z 0 − Z 0 and Z 0 − γ self energies, for which we import the same expressions that hold in the standard on-shell scheme for the SM, see e.g. [47,78], but of course including in their calculation the loop contributions from the 2HDM Higgs bosons (cf. Figs. 3 and 4)  we can split the one-loop amplitudes into the following (UV-finite) subsets; and similarly for the case of the heavy Higgs boson, e + e − → A 0 H 0 . In the above equations, we understand that each amplitude is supplemented with the corresponding counterterm. The complete expression M one−loop = M (1) + δ M (1) is free of UV divergences, as we have explicitly checked. Moreover, at this point of the discussion we are keeping track explicitly of the λ 3H factors, in order to highlight those contributions which can be enhanced by the 3H self-couplings. For example, among the interference terms between the tree-level amplitude and the one-loop amplitude (67) we have the class of terms O(α 3 ew ) and also the class of terms O(α 2 ew λ 2 3H ), both classes of the same order in perturbation theory, although the latter are expected to be comparatively enhanced. We will check it numerically in the next section.
Taking into account the above expressions, the overall (one-loop corrected) scattering amplitude for the process e + e − → A 0 h 0 reads as follows:  CP-odd field on account of the field renormalization condition (44), this is not possible for h 0 . The reason is that we use a WF renormalization constant for each Higgs doublet, not for each field, see Eq. (21), and thus the relations (22) give no further room to normalize to one the residue of the propagator for the other fields. The nontrivial finite renormalization contributionẐ h 0 = 1 for the other external Higgs boson in this process, h 0 , must be computed explicitly, the result beinĝ Retaining from this expression only the O(α ew λ 2 3H ) contributions to the scattering amplitude, we are left with At the same time, our definition of the renormalized CP-even mixing angle implies that the h 0 − H 0 contribution to the one-loop amplitude (67) vanishes identically: In turn, the finite correction associated to the on-shell A 0 leg can be expressed in terms of the renormalized mixing self-energies in the guise where in the last step we have made use of the Slavnov-Taylor identity (45). The 4-momentum p µ stands for an outgoing A 0 boson. Finally, the expression for the renormalized mixing self-energy A 0 − G 0 appearing in Eq. (72) is obtained from (28) after substituting in it the explicit form of the counterterm (59).
Notice that since Eq. (44) decrees that the real part of the functionΣ A 0 Z 0 (M 2 A 0 ) must vanish, in the above amplitude (72) the only non-vanishing contribution shall stem from the imaginary part, ℑmΣ A 0 G 0 (M 2 A 0 ). A similar expression may be derived for the H 0 A 0 channel: The complete e + e − → A 0 h 0 amplitude at order O(α 3 ew ) can finally be rearranged in the following way: where the tree-level amplitude takes the explicit form: In the above expression, the notation reads as follows: e + , with 4-momentum p 1 and helicity η 1 ; e − , with 4momentum p 2 and helicity η 2 ; A 0 , with 4-momentum k 1 ; and h 0 , with 4-momentum k 2 . We have also introduced the left and right-handed weak couplings of the Z 0 boson to the electron, g L = (−1/2 + s 2 W )/c W s W , g R = s W /c W and the left and right-handed projectors P L,R = (1/2)(1 ∓ γ 5 ). Let us also mention that we have neglected the contribution of the Z 0 -width, as it is completely irrelevant at the working energies of the linear colliders (s ≫ M 2 Z 0 ). Finally, the total cross section is obtained after squaring the matrix element, performing an averaged sum over the polarizations of the colliding e + e − beams and integrating over the scattering angle: where we have introduced the standard definition of the Kählen function λ(x, y, x) = x 2 +y 2 +z 2 −2xy−2xz−2yz.
Before closing this section, let us stress once more that a similar set of formulae hold for the companion process e + e − → A 0 H 0 , which can be related to the equations above through the correspondence cos(β−α) → − sin(β− α) and M h 0 → M H 0 .    In this section, we describe in detail the results that are obtained from the numerical analysis of the processes e + e − → h 0 A 0 /H 0 A 0 at the one-loop level. The basic quantities of interest are two: i) the predicted cross section (76) at the Born-level σ (0) and at one-loop σ (0+1) ; and ii) the relative size of the one-loop radiative corrections, which we track through the parameter In practice, the computation has been performed with the help of the standard algebraic and numerical packages FeynArts, FormCalc and Looptools [80] for the generation of the Feynman diagrams, the analytical calculation and simplification of the scattering amplitudes as well as the numerical evaluation of the cross section. Aiming at a wider survey of the regions of phenomenological interest in the 2HDM parameter space, we will explore the relevant range of values for both cos(β − α) and sin(β − α) by picking up the representative configurations quoted in Table III. In particular, in the case α = β − π/2 the h 0 h 0 h 0 coupling takes on the SM form, as warned before. Furthermore, the case α = π/2 corresponds to one of the so-called fermiophobic scenarios; namely, for type I 2HDM, the lightest CP-even Higgs couples to all fermions as cos α/ sin β times the SM coupling (cf . Table  III); therefore, if α = ±π/2 the lightest Higgs decouples from all the fermions 9 . Full fermiophobia is of course not possible for type-II models (although it is partly possible). In particular, for the MSSM this phenomenon is impossible and its detection would be a direct signal of non-SUSY physics. Let us clarify, however, that per se the fermiophobic scenarios are not particularly relevant for the current analysis because the Yukawa couplings do not play a central role here (as we will discuss below). However, our cross-sections are indeed sensitive to the particular value of α that we select, and in this sense the two fiducial choices α = β and α = π/2 represent two distinctive regimes for our numerical analysis, as we will see.
For practical reasons, throughout our calculation we will be interested in projecting out the large radiative corrections from regions of the parameter space for which the predicted production rates at the tree-level are already sizable. This means that we shall mainly dwell on regimes nearby the optimal value of the tree-level coupling for the lA 0 channel (β ≃ α) and for the H 0 A 0 one (β ≃ π/2 + α).
In Table IV we quote the different Higgs-boson mass spectra that we shall consider in our numerical analysis. These mass sets have been designed to cover a wide span of phenomenologically motivated regimes. Because of their mass splittings, Sets I-III can only be realized within a general 2HDM, whereas Sets IV-VI reflect characteristic benchmark scenarios of Higgs boson masses within the MSSM; they have been generated through the MSSM parameter inputs specified in Table V. Note that Sets I to III (similarly, Sets IV to VI) are organized through increasing values of the Higgs boson masses. Furthermore, let us also highlight that, due to the presence of a light charged Higgs boson, Sets I-II and IV-V are adequate only for type-I 2HDM, whereas Sets III and VI are valid either for type-I or type-II 2HDM. I  100  150  140  120  Set II  130  150  200  160  Set III  150  200  260  300  Set IV  95  205  200  215  Set V  115  220  220  235  Set VI  130  285 285 300   Table IV, thus mimicking the characteristic one-loop mass splittings within the Higgs sector of the MSSM. The numerical value for these masses has been obtained with the aid of the program FeynHiggs by taking the full set of EW corrections at one-loop [81]. Universal trilinear couplings (At = A b = Aτ ) and GUT relations for the gaugino soft-SUSY breaking mass terms are assumed throughout.
In turn, the value of λ 5 is severely restrained by the theoretical constraints stemming from the tree-level unitarity and vacuum stability, i.e. from equations (17) and (19)- (20). The bottom line is that λ 5 > 0 is strongly disfavored, and the permissible values lie mostly in the range The value of tan β that is preferred so as to optimize the cross-sections in the above range is tan β ≃ 1. We will carefully derive these constraints in what follows and shall examine their implications on the numerical analysis of the cross-sections under study.
B. e + e − → h 0 A 0 A fairly generous overview of the numerical results for the Higgs-pair production through e + e − → h 0 A 0 at one-loop is provided systematically in Figs 12 -14 and in Table VI. In the latter, we display the predicted values for the total cross-section σ(e + e − → h 0 A 0 ), together with the relative size of the associated quantum corrections (77) for each of the Higgs-boson mass sets in Table IV. The results are obtained for fixed tan β = 1 and different values of α indicated in the table. The fiducial ILC center-of-mass energy is taken to be either √ s = 500 GeV or √ s = 1 TeV. For each mass set, λ 5 is allowed to take its largest (negative) attainable value indicated in Table VII. As we will confirm throughout the numerical analysis, this prescription insures maximally enhanced quantum corrections.
The corresponding loop-corrected cross sections σ (0+1) (e + e − → h 0 A 0 ) stay in the approximate range 2 − 30 fb for √ s = 500 GeV -this would entail up to barely 10 3 − 10 4 events in the standard segment 500 fb −1 of integrated luminosity. At √ s = 1 TeV, however, the predicted yields deplete down to O(1) fb in most cases, although this would still give a turnover of a few hundred events at the end of the luminosity shift.
As for the radiative corrections themselves, they are also explicitly quoted in Table VI. One can see that they can be remarkably large, being either positive (at "low" energies, viz. for √ s ∼ 500 GeV) or negative (for √ s ∼ 1 TeV and above). In the most favorable instances, such corrections can boost the cross-section value up to δ r ∼ 100 % at √ s = 500 GeV, whereas the corresponding suppression in the high energy range ( √ s = 1 TeV) can attain δ r ∼ −80% owing to a severe destructive interference between the tree-level and the one-loop amplitudes.
In Figs 12 -14 we explore in more detail the behavior of the radiative corrections and their interplay with the theoretical constraints associated to the perturbative unitarity and vacuum stability conditions. Although the range tan β 1 is usually the preferred one from the theoretical point of view, we entertain the possibility that tan β can be slightly below 1 in order to better assess the behavior around this value. In these plots, we depict δ r (in %) for two choices of energies, √ s = 500 GeV (left panels) and √ s = 1 TeV (right panels). For each of these figures we show independent contour plots corresponding respectively to the Sets II-IV of Higgs boson masses (cf. Table IV) in the tan β − λ 5 plane. The first of these sets involves an assorted spectrum of Higgs boson masses, which is impossible to reproduce in the MSSM, whereas the second set closely mimics a typical MSSM-like Higgs mass spectrum, in which we recognize the characteristic (approximate) degeneracy of the heavy Higgs boson masses: M H 0 ≃ M A 0 ≃ M H ± . We find that these scenarios are well representative of the phenomenological trends shown by all Higgs mass spectra under analysis. By setting α = β we optimize the tree-level A 0 h 0 Z 0 coupling and in this way we insure a sizable value of the lowest order cross-section, which is of course a good starting point for having a chance to eventually measure quantum corrections on it.
From these plots, it is eye-catching that the unitar- FIG. 13: As in Fig. 12, but for the Set III of Higgs boson masses (cf. ity constraints are the most restrictive ones, and tend to disfavor large values of tan β > 1 or values tan β ≪ 1. This is natural because the triple and quartic Higgs selfcouplings rocket fast with large and small tan β, see Eq. (12) and Table II. At the same time, these unitarity constraints place a lower bound on another highly sensitive parameter of this study, viz. λ 5 , whereby λ 5 cannot be smaller than −11 or thereabouts (in other words, |λ 5 | 11). Besides, this parameter is sharply stopped from above near λ 5 0 as a direct consequence of the vacuum stability condition (cf. the shaded band in the upper part of the plots). Thus, the combined set of bounds build up a characteristic physical domain in the parameter space of Figs 12 -14 which pretty much looks like a deep valley flanked with sharp cliffs on each side, and centered at tan β = 1. As a result, large departures from this central value are incompatible with large values of |λ 5 |. In particular, the range from tan β = 3 onwards is circumscribed to a narrow-edged neck around λ 5 ≃ 0. In addition, the heavier the Higgs mass spectrum, the more severe the unitarity bounds are. Let us also remark that all these constraints over |λ 5 | exhibit FIG. 14: As in Fig. 12, but for the Set IV of Higgs boson masses (cf.   a mild dependence on the actual value of the CP-even mixing angle α in the regimes considered here, whereas our cross-sections do vary significantly with it. In short, from the analysis of Figs. 12-14 two basic regimes of phenomenological interest can be sorted out: 1) On the one hand, we find scenarios in which λ 5 < 0 and where this coupling can stay moderately large in absolute value (viz. |λ 5 | ∼ 5 − 10) while tan β ≃ 1; this would correspond to the relatively narrow allowed stretch on the left-hand-side of these plots. In such configurations, a subset of 3H self-couplings are remarkably singled out -the more negative is λ 5 , the greater is the coupling strength. Such enhancement, which is transferred to the e + e − → A 0 h 0 amplitude through a set of Higgs boson-mediated one-loop vertex diagrams (see Fig. 5), translates into sizable quantum corrections which become more vigorous with growing |λ 5 |. For instance, from Fig. 12 we see that around tan β ≃ 1, δ r reaches ∼ +10 % at λ 5 = −6, whilst it becomes +30 % for λ 5 = −10. Furthermore, we encounter that positive radiative corrections as big as δ r ∼ +40% at √ s = 500 GeV may switch drastically into large negative effects of the same order ∼ −40% (hence suppressing significantly the crosssection) when the center-of-mass energy is increased up to √ s = 1TeV -see e.g. the right panel of Fig. 12. Obviously, positive radiative corrections are preferred in practice because they invigorate the physical signal. Therefore, the startup energy √ s = 500 GeV of the ILC seems to be the ideal regime to probe these particular quantum effects rather than moving to higher energies;   2) On the other hand, a very different behavior occurs in the neck-shaped region of the contour plots at moderate tan β > 1, as it is apparent in Figs. 12-14. Here the enhancement capabilities of the leading 3H self-couplings are rather meager with respect to the aforementioned large-|λ 5 | scenario. The reason is that for tan β 3 the coupling |λ 5 | is nailed down to take very low values.
In the absence of a significant yield from the Higgs boson-mediated loop corrections, we might expect alter-native sources of quantum effects to pop up. The first ones coming to mind are of course those originating in the Higgs-fermion Yukawa sector. However, it turns out that the Yukawa interactions are not particularly efficient in the regions explored in Figs. 12-14. Indeed, from Table III it is apparent that none of the Higgs boson couplings to fermions is spurred on for tan β 1 (in any of the two canonical 2HDM's). Only the regions with tan β < 1 -which are disfavored by the general theoretical expectations -and cos α ∼ 1 could contribute here mainly via the Higgs-top interaction, whose strength is h 0 tt ∼ m t cos α/ sin β regardless of considering type-I or type-II models. Additionally, the Higgs-bottom Yukawa coupling (∼ m b tan β within type-II 2HDM), could furnish a competitive source of enhancement, as it does e.g. in the MSSM. In the present case, however, no trace is left of such effect owing to our main focus on the range tan β 1, although this situation could change in the region tan β < 1 (see below).
Closely related with the previous comment, we should also stress the fact that there are no outstanding differences between the obtained results from type-I and type-II 2HDM's in the context of this analysis. Indeed, to start with the Higgs-quark Yukawa interactions do not  Table IV. Shown are the results obtained within three different values of λ5, at tan β = 1 and for α = β (left) and α = π/2 (right).  VIII: Relative contribution to the total cross section δr,i = (σ (0+i) − σ (0) )/σ (0) from each topology of one-loop diagrams normalized to the tree-level rate. The results are derived at fixed √ s = 500 GeV for tan β = 1, α = β and choosing the value of λ5 that optimizes the quantum corrections for each mass regime, see Table VII. enter at the tree-level in the processes (3) under consideration; second, the trilinear couplings are common to both types of models; and third, we have already emphasized that the one-loop differences that could manifest through the distinct form of these Yukawa couplings in type-I and type-II models are virtually obliterated in the physical region permitted by the constraints. All in all, perhaps the most distinctive feature between both models boils down, in this context, to just the requirement that M H ± 300 GeV for type-II 2HDM's as a result of the low-energy B-meson physics constraints, which tend to substantially raise the average mass of the Higgs spectrum for this class of models -cf. Sets III and VI of in Table IV. The basic differences in the physics of e + e − → A 0 h 0 in both models thereby narrow down to mass differences between their characteristic Higgs-boson mass spectra 10 .
A very similar pattern is met in the other explored scenarios, the only differences being ascribable to the distinct Higgs boson mass spectrum that each of these scenarios accommodate. In particular, heavier Higgs boson masses may have a twofold impact: i) dynamical, in that they soften the influence of Higgs-mediated one-loop diagrams -and hence the strength of the 3H self-couplings -inasmuch as these contributions become suppressed by inverse powers of the Higgs masses; and ii) kinematical, in that they shift the peak of the σ( √ s) curve, since the latter is correlated with the production threshold of the Higgs pair.
All the above mentioned phenomenological features are transparent in Figs. 15-17, in which we explore the evolution of the cross section as a function of the centerof-mass energy. We include, in each plot, the tree-level contribution σ (0) and the loop-corrected one σ (0+1) for different values of λ 5 . Worth noticing is that by testing the influence of λ 5 on the cross-sections, we are testing the bulk capability of the 3H self-interactions, and hence Therefore, a shift of M H ± slightly upwards restores the possibility of tan β ∼ 1 at, say, 2σ without altering significantly our results (as we have checked.) their potential fingerprint on the cross-sections at the quantum level. In the lower panels, we track the related behavior of the quantum correction δ r with √ s for the same set of λ 5 values. The plots are generated at fixed tan β = 1 and for two different values of the tree-level A 0 h 0 Z 0 coupling: namely, α = β (where the tree-level coupling is maximum) and α = π/2 (corresponding to the aforementioned fermiophobic limit, in which h 0 fully decouples from the fermionic sector for type I models). In Figs. 15,16 and 17, we display the numerical results obtained for Sets II, III and IV of Higgs boson masses, respectively. A similar numerical output is obtained for the remaining mass sets in Table IV, showing no relevant departure from the phenomenological trend that we have recorded as yet.
Let us focus e.g. on Fig. 15 for a while. The treelevel curve exhibits the expected behavior with √ s, as it scales with the s-channel Z 0 -boson propagator ∼ 1/(s − M 2 Z ). The maximum cross section is achieved at ∼ 500 GeV. The same pattern is followed by the full loopcorrected cross-section, although the dependence with √ s is sharper. In fact, the range where δ r > 0 is much briefer than the one where δ r < 0, although the former is characterized by very significant quantum effects; so much so, that they may skyrocket up to 100 % (at the largest allowed values of |λ 5 |) around the critical energy domain where √ s = 500 GeV. Such effects could hardly be missed in the first runs of the ILC. In contrast, when we move from √ s = 500 GeV to √ s = 1000 GeV the resulting cross section σ (0+1) drops by at least 50%. The impact of the enhanced 3H self-couplings can be easily read off from the dramatic differences in the one-loop production rates when the value of λ 5 is varied in the theoretically allowed range (78). Thus, when |λ 5 | is pulled down at fixed tan β, the quantum effects are rapidly tamed, in correspondence with the fact that the 3H self-couplings are no longer stirred up.
The drastic decline of the triple Higgs boson selfinteractions in this domain cannot be counterbalanced by the contribution of the Yukawa couplings, as they are not very significant in the region of the parameter space where tan β is of order one. At this point, it comes to mind our former observation about the potential role that the region tan β < 1 could play. Even though it is not favored on theoretical grounds, we can not exclude the possibility that O(0.1) < tan β < 1, and in fact the unitarity and vacuum stability bounds cannot dismiss in block this range when λ 5 is sufficiently small. Thus, for example, in the conditions of Fig. 12, and under the assumption that λ 5 = 0, we have checked that tan β ∼ 0.2 is allowed by the aforesaid constraints, and in this region we encounter quantum effects as large as +30%. They are the result of the combined effort from the diagrams involving the top quark Yukawa coupling and those involving the trilinear Higgs boson couplings -which, in this case, become enhanced by a pure tan β < 1 effect, rather than by an oversized λ 5 .
It is also interesting to track carefully the fractional payoff δ r = (σ (0+i) − σ (0) )/σ (0) from the various oneloop topologies (i = 1, 2, 3, ...) of Feynman diagrams contributing to the process e + e − → h 0 A 0 . The corresponding results are displayed in Table VIII under the given conditions. Included there are the gauge invariant subsets of diagrams of various types, namely Z 0 − Z 0 /γ vacuum polarization effects, box diagrams, the renormalized A 0 h 0 Z 0 and e + e − Z 0 vertex corrections, the loop-induced form factor γh 0 A 0 , and finally the finite WF corrections to the external h 0 leg. We can see that the gauge boson self-energies provide tiny contributions of order 1% at most. The box diagrams, in their stead, are nonnegligible (of order 10%), albeit negative. At the same time we detect mutually destructive effects between e.g. the finite WF correction and the e + e − Z 0 vertex.
Most conspicuously, and above all, we revalidate the dominant role of the Higgs boson mediated loopdiagrams (triggered by 3H self-interactions) as the truly leading source of quantum corrections to the cross- section. These dominant effects are concentrated on the fundamental A 0 h 0 Z 0 vertex. This vertex, which is purely gauge-like at the lowest order of perturbation theory, becomes drastically promoted at the one-loop level as a result of the dominant quantum effects sourced by the 3H self-interactions. We may cast the overall renormalization of such vertex as follows 11 : where Γ 1,3H denotes the leading subset of (Higgsmediated) one-loop corrections. The latter can be encapsulated in a UV-finite subset of one-loop Feynman diagrams and subsequently reabsorbed into a modified or effective A 0 h 0 Z 0 coupling. Simple power counting and educated guess provides the following estimate: 11 A similar discussion is of course also applicable to the complementary A 0 H 0 Z 0 coupling and the corresponding production process e + e − → H 0 A 0 , see the next subsection.
where f (M 2 h 0 /s, M 2 A 0 /s) is a dimensionless form factor that accounts for the complete one-loop amplitude. For all its cumbersome structure, we do know at least that f dies out with the center-of-mass energy of the process and encodes the sign flip of the radiative corrections when we increase the center-of-mass energy from √ s = 0.5 TeV to √ s = 1 TeV. The prefactor 1/16π 2 stands for the trademark numerical suppression attached to one-loop form factors, whilst λ 3H denotes a generic 3H self-coupling. Noteworthy is the fact that this coupling can reach values as high as λ 3H /M W ≃ |λ 5 |/e ≃ 30 and still under the stringent constraints imposed by unitarity.
It is enlightening to try to estimate the expected size of the maximum corrections. From Eq. (80) we can roughly gauge the maximum relative contribution, e.g. at the startup energy √ s = 500 GeV: where < ... > stands for the various operations of averaging and integration of the squared amplitudes. In the above equation, we have neglected the square of the one-  Fig. 19, but for Set IV of Higgs boson masses, cf. Table IV. loop amplitude (see, however, just below) and set f ∼ 1, which is a reasonable assumption since the current value of the center-of-mass energy, √ s = 500 GeV, is far from the range in which one-loop corrections are proved to be negligible. The latter result of Eq. (81) falls in the right ballpark of the optimal values for δ r that we have identified from our exact numerical analysis, for Higgs boson masses of the order of a few hundred GeV and maximum allowed |λ 5 | values (so as to push λ 3H close to the unitarity border). By means of the same formula we can also estimate the contribution to δ r arising from the squared of the one-loop amplitude, which we have neglected before. We find thus roughly 10% of the one-loop leading effect. It turns out that this is again in good agreement with the numerical calculation 12 . All in all we conclude that, as long as 3H couplings are sizable, quantum effects will be manifest -and potentially very large. Some more comments are in order concerning important features contained in Figs. 15-17. As we move from 12 We will further dwell on the potential relevance of such O(α 4 ew ) effects in Section VII. Fig. 15 to 16, hence as we transit from a relatively light to a heavier Higgs boson spectrum (cf . Table IV), we encounter that the resulting cross sections, both at the leading and at one-loop level, are pulled down by roughly a factor of 2. Unsurprisingly, this is signaling the exhaustion of the available phase space. By the same token, the peak-shaped maximum of the σ( √ s) curve gets shifted forward by about 100 GeV, in correspondence to the production threshold of the Higgs pair. Moreover, the heavier the Higgs bosons, the more suppressed the Higgs-mediated loops become, and hence the lower is the influence of 3H self-couplings. This is precisely what we confirm numerically in Fig. 16 (lower panels), in which we can appreciate how δ r reaches ∼ 40% at most -in contrast to δ r ∼ 100% that was spotted for Set II (cf. Fig. 15). In turn, Set VI (in Fig. 17) embodies again a rather light Higgs boson spectrum which resembles that of Set II, even though the present choice of masses is obtained upon artificially engineering a SUSY setup. Here, owing to the fact that M H ± is heavier than in Set II, unitarity conditions settle a more stringent upper bound on |λ 5 |, and hence on the maximum enhancement of the triple neutral Higgs couplings. As a result, radiative corrections in this case are still sizable, but retreat to the ∼ 40% level.
Of great relevance is the study of the behavior of σ(e + e − → h 0 A 0 ) as a function of the Higgs boson masses. We present these results conveniently in Figs. 18-20, √ s = 500 GeV √ s = 1 TeV α = β − π/2 α = β − π/3 α = β − π/6 α = 0 α = β − π/2 α = β − π/3 α = β − π/6 α = 0  where again we superimpose both the tree-level and the loop-corrected cross-sections for the process e + e − → h 0 A 0 . The dependence with M h 0 is indicated on the left panels and that with M A 0 on the right panels. The numerical results are shown for the Sets II,III and IV and have been obtained, once more, by setting tan β = 1, α = β, and making allowance for |λ 5 | to take on the highest permissible value in each set (see Table VII).
The center-of-mass energy is settled at the fiducial value √ s = 500 GeV.
The plots under consideration evince that the raise of the Higgs boson masses exerts a twofold effect on σ, one kinematical and the other dynamical. The former simply means that owing to the reduction of the available phase space the cross-section obviously falls down. We can check it in the figures, where we see that both the tree-level σ (0) and the loop-corrected σ (0+1) indeed decrease monotonically with growing M h 0 and M A 0 . As far as the quantum corrections themselves are concerned, the situation is a bit more subtle. Although heavier Higgs bosons imply larger propagator suppressions of the Higgs-mediated one-loop diagrams, as we have just described in the previous paragraph, the 3H self-couplings can partially offset this situation since they get invigorated in such circumstances. This counterbalance feature is manifest in the explicit expressions for the 3H couplings in Table II -the uttermost limitation to their growing being the unitarity constraints. In particular, 3H couplings involving neutral fields (which can be shown to carry the dominant effect on e + e − → h 0 A 0 ) turn out to grow either with M h 0 or M A 0 . Although both effects partially balance each other, we are ultimately left with the expected decoupling behavior in the limit of large Higgs boson masses, with the proviso that such decoupling takes place in a softer way than that of the Bornlevel cross section. In fact, in Figs 18-20 we can verify that the one-loop corrected cross-section remains most of the time on top of the tree-level result when we increase the Higgs boson masses, and moreover the slope of the former is clearly milder, particularly for M h 0 200 GeV. This translates into a growing trend of δ r up to values of M h 0 ≃ 275 GeV. Let us also mention in passing that the sudden spike in the δ r (M A 0 ) curve at M A 0 ≃ 340 GeV (which is also barely visible for σ) corresponds to the tt pseudo-threshold in the one-loop vertex amplitude (cf. the first row of diagrams in Fig. 5).
Analogous trends are encountered for sets III and IV of Higgs boson masses, the only reportable difference being that the decoupling behavior of σ (0+1) with M h 0 , M A 0 is more pronounced, signaling that the increase of these masses does not boost the 3H self-couplings as efficiently as in the case of Set II. The shaded regions in Fig. 19-20 correspond to mass ranges already excluded by experimental measurements of the δρ parameter (cf. section III).
C. e + e − → A 0 H 0 On general grounds, the basic phenomenological features are common to the h 0 A 0 and H 0 A 0 channels, with the only proviso that, in the latter case, the resulting production rates are slightly smaller simply because M H 0 > M h 0 . The final cross-sections are nevertheless comparable and reach up to a few dozen fb. In Table IX we collect, in a nutshell, the essential results from a dedicated numerical analysis of σ(e + e − → H 0 A 0 ) for all sets of Higgs boson masses (cf. Table IV) and different values of the tree-level A 0 H 0 Z 0 coupling, again at fixed tan β = 1 and using the largest permissible value of |λ 5 | in each set (see Table VII). The relative size of the quantum corrections δ r is also included in Table IX. As in the previous channel, we have used fiducial center-of-mass energies √ s = 500, 1000 GeV respectively. The maxi- mum production rates at √ s = 500 GeV are achieved for relatively light Higgs boson masses (Set I) and for α = β − π/2, in which case the tree-level A 0 H 0 Z 0 interaction is not suppressed, see Eq. (61). Within this setup, one is left with σ (0+1) ∼ 40 fb. Indeed the production rates lie at the O(10 fb) level in a broad range around α = β − π/2 for all the mass sets under analysis. If we nevertheless consider α = β − π/6, when the tree-level coupling gets suppressed by sin(β − α) = 1/2, the resulting cross sections are still at the level of a few fb -each femtobarn rending 500 events per 500 fb −1 of integrated luminosity.
Radiative corrections thus may leave a solid imprint into the final predicted σ (0+1) of the current process. In particular, for α = β−π/2 quantum effects rubber-stamp a characteristic 50−100% boost upon the tree-level crosssection, which could hardly be missed. Other choices of α amount to negative one-loop contributions in which A 0 H 0 Z 0 is dramatically suppressed. For instance, in the case of Set I for α = β − π/6, it renders δ r ∼ −65%. For center-of-mass energies from around √ s = 1TeV onwards, a powerful destructive interference is seen to be under way between the Born and the one-loop amplitudes. This translates into large (negative) quantum cor-rections which may well reach δ r ∼ −50 − 100%, thus literally stamping out the signal and pulling it down at the ∼ fb level or even below. To better illustrate this feature, in Figs. 21-23 we plot the behavior of σ(e + e − → A 0 H 0 ) (top panels) and δ r (bottom panels) as a function of √ s for the representative Sets II-IV of Higgs boson masses of Table IV and for different values of the tree-level coupling; viz. α = β − π/2 (left panels), which corresponds to the maximum of the tree-level A 0 H 0 Z 0 coupling, and sin α = 0 (right panels), corresponding to the fermiophobic limit for H 0 boson within type-I 2HDM (cf. Table  III). Again we fix tan β = 1 and we vary λ 5 within the range allowed by unitarity and vacuum stability restrictions. As it was already observed in the analysis of the h 0 A 0 channel, heavier Higgs-boson spectra furnish smaller production rates. We refrain from extending the discussion, which is entirely similar to that of e + e − → h 0 A 0 .
To our knowledge, no pre-existing study of the oneloop effects on the neutral Higgs boson channels (3) can be found in the literature. Let us, however, briefly comment on the prospects for the charged Higgs-pair channel, e + e − → H + H − . This was considered some time ago in references [41][42][43]. These old studies were devoted  to the computation of σ(e + e − → H + H − ) at one-loop in both the MSSM and the 2HDM. In the latter case, they also uncover very large quantum effects (of order 100% or even higher 13 ) which are attributed to enhanced 3H self-interactions involving charged Higgs bosons. Such quantum effects may lie well above the typical MSSM counterparts, which are found to render up to |δ r | ∼ 20%, and mostly negative. Let us, however, remark that the e + e − → H + H − channel could operate at sizable crosssections levels only for type-I models since the type-II models require a rather heavy mass spectrum. 13 We should nevertheless put a note of caution here in that the unitarity constraints that were employed in Ref. [41][42][43] are substantially less restrictive than those considered throughout our work. In this regard, an upgraded calculation of σ(e + e − → H + H − ) is mandatory before we can judiciously compare with the results that we have presented here for the neutral Higgs boson channels. However, due to the increase of the charged Higgs boson mass bounds in the last few years, the maximum cross-section rates should correspondingly fall down, thereby making this process less competitive.

VII. DISCUSSION AND CONCLUSIONS
In this work, we have undertaken a thorough investigation of the pairwise production of neutral Higgsbosons within the context of the general two-Higgsdoublet model (2HDM) at the ILC/CLIC linear colliders (linac) through the basic processes e + e − → h 0 A 0 /H 0 A 0 . To the best of our knowledge, this study is the first one existing in the literature which addresses these important processes at the one-loop level. To sum up the main objectives of our work, we can assert that we have concentrated on the following three tasks: i) we have, on the one hand, singled out the most favorable regimes for the Higgs-pair production processes once the radiative corrections are taken into account, specifically within the framework of the on-shell renormalization scheme with an appropriate definition of tan β and of the CPeven mixing angle α at the quantum level; ii) next, we have carefully quantified the importance of the leading quantum corrections, and not only so, we have pinpointed also the origin and relative weigh of the various sources; and, finally, iii) we have correlated the most powerful source of quantum effects with the strength of the Higgs boson self-interactions (essentially the triple Higgs boson self-couplings λ 3H ), to which the processes e + e − → h 0 A 0 /H 0 A 0 are highly sensitive through Higgsmediated one-loop diagrams. With these aims in mind,  Among the various subtle triggers that control the size of these important couplings, one of them stands out with overwhelming supremacy over the others, and this is the λ 5 coupling. This is the only one parameter in the list of seven free parameters -cf. Eq. (9) -conforming the structure of the general (CP-invariant) Higgs potential of the 2HDM that cannot be ultimately traded for a physical mass or a mixing angle. This parameter is, in addition, the only Higgs self-coupling that need not be renormalized at one loop in the current framework. Despite its special status of being apparently unfettered and without obvious physical connotations, one quickly discovers that it becomes severely restrained by subtle bounds connected with the quantum field theoretical consistency of the model, namely the stringent constraints imposed by tree-level unitarity and vacuum stability conditions. These restrictions, together with the limitations dictated by custodial symmetry, and of course also those emerging in the light of the direct collider bounds on the masses of generic Higgs boson models, not to mention the indirect phenomenological restraints (in particular, the upper bound on the charged Higgs mass in type-II models ensuing from radiative B-meson decays), do indeed eventually determine a fairly compact region of parameter space. In particular, λ 5 is compelled to be mostly negative and in the range |λ 5 | O (10). At the same time, tan β cannot be large when |λ 5 | is also large without violating the constraints. The bottom line is that the value that maximizes the processes under consideration is a rather equitable one: tan β ≃ 1.
Notwithstanding the sharp limits placed upon the physical region of the parameter space, our explicit calculation of radiative corrections has shown -and we feel it is truly remarkable -that the general 2HDM models are still able to unleash a good deal of their "repressed" power. If only we could catch a smoking gun of this stupendous potential through the accurate measurements that a linear collider should be able to render on fundamental processes like e + e − → h 0 A 0 /H 0 A 0 , we might be on the verge of experiencing an intense episode of Higgs boson physics beyond the SM. In the following, we describe the basic results and strategies that we could follow to this effect.
We have identified three main scenarios of potential phenomenological interest in the corresponding parameter space. In all of them we observe no remarkable dependence either on the details of the Higgs mass spec-trum nor on the pattern of Higgs-fermion Yukawa couplings (either for type-I or type-II 2HDM), not even on the particular channel h 0 A 0 /H 0 A 0 under study. In what follows we describe each one of these scenarios and spell out clearly their differences, in particular we provide the characteristic potential size of the relative quantum correction δ r -cf. Eq. (77) -typically associated to them (cf. Figs. 12-14): • Scenario 1. To start with, regions of tan β 1 and large |λ 5 | (specifically, λ 5 < 0, |λ 5 | ≃ 5 − 10) support the bulk of the radiative corrections. They may well reach δ r = ±(50 − 100)% at the largest attainable values of |λ 5 |. Furthermore, these effects critically depend on the actual size of |λ 5 | and on the center-of-mass energy √ s (cf. Figs. 15-17), to wit: δ r is positive at "low" center-of-mass energies around √ s = 500±100 GeV, while it is negative for √ s > 600 GeV and certainly also in the uppermost foreseeable energy segments √ s = 1 − 3 TeV, usually reserved for an upgraded ILC and specially for CLIC. Thus, interestingly enough, quantum corrections turn out to drastically invigorate the cross section of the basic processes e + e − → h 0 A 0 /H 0 A 0 near the standard booting energy of the ILC, while they tend to suppress the two-body Higgs boson signal for a center-of-mass energy some 20 − 30% beyond this initial regime -and of course for the TeV range and above -as a result of constructive/destructive interference effects, respectively, between the tree-level and the one-loop amplitudes. Incidentally, we clarify that this optimal scenario cannot be iterated for large and positive values of λ 5 because the combined set of constraints (mainly the vacuum stability bounds, in this case) preclude most of the region λ 5 > 0; • Scenario 2. On the other hand, we find a tail of subleading effects within a band of the (tan β, λ 5 ) subspace at moderate tan β = 1.5 − 5 and |λ 5 | 2 wherein the maximum radiative corrections turn out to stagnate and remain barely at the level of ∼ 10% or less. Should the parameters of the model inhabit this region, the practical possibilities to detect these effects would obviously become thinner than in the previous case. Still, 5 − 10% effects are not that alien to the high precision standards planned for a linac collider, and hence it should not deter us from the searching task. At this point, the detection strategies associated to the Higgs boson decay patterns could be essential (see below); • Scenario 3. Finally, there is the region of parameter space where the coupling λ 5 is very small or it just vanishes. Here the quantum corrections can still have a chance provided tan β < 1, and indeed we have found domains of this kind preserving all the basic constraints. While this region is not the most favored one, neither from the theoretical point of view nor from the experimental data, specially if tan β ≪ 1, we should keep it in mind as a possibility. Even for moderately small tan β, say for tan β 0.2, it can be the source of sizeable quantum effects of order 10 − 30% (and positive), which are larger than those in Scenario 2 and, in some cases, even comparable to those in Scenario 1. It is worth noticing that, in this region, there is an interesting collaborative interplay between the Higgs boson trilinear self-couplings and the top quark Yukawa coupling (which is the same in type I or type II 2HDM's). Only in this scenario a Yukawa coupling could play a significant quantitative role on equal footing with the trilinear Higgs boson self-couplings.
Obviously, the success of the whole Higgs boson search endeavor will also depend on the pattern of distinctive signatures available to these particles in the final state. For instance, for M h < 2M V 180 GeV, we should basically expect back-to-back pairs of highly energetic (roughly E ∼ √ s/4 100 GeV) and collimated b-quark and/or τ -lepton jets from A 0 , h 0 → bb/τ + τ − which, in contrast to the LHC, should not be overshadowed by the large QCD background which is inherent in hadronic machines. And of course similarly with the A 0 H 0 final state. This much clearer environment notwithstanding, we should be aware of the fact that other characteristic processes of linac physics, such as the fourjet cross-section, might be large enough to partially obscure the Higgs-pair production signal. However, the precise analysis of the corresponding signal distributions, which would be required in order to completely assess the real experimental possibilities, is beyond the scope of the present study.. And of course the same reasoning applies to the A 0 , H 0 final state. For M h > 2M V , on the other hand, signatures from h 0 , H 0 → V V and A 0 → bb/τ + τ − could also play a role (depending on the particular choice of α and β, cf. Table III). They could ultimately lead to signatures with two (up to four) charged leptons against a bb or τ + τ − pair in the final state (from W ± → ℓ ± + missing energy and, specially, from Z → ℓ + ℓ − ). Even if the branching ratios of gauge bosons into leptons are very small (B(V → ll) ∼ 0.03),the predicted number of Higgs-pair events is large enough so that, in the most favorable scenarios, one could collect O(10 2 ) events for these complementary signatures, which could be instrumental for a final tagging and identification of the Higgs bosons. Notice that we use the fact that the CP-odd A 0 "always" decays into bb or, to a lesser extent, to τ + τ − (even for tan β = 1) as it cannot couple to gauge bosons. So the patterns of signatures are rather characteristic and could not be missed in normal circumstances. But there is more than meets the eye in the whole strategy game, as we shall see later on.
In the most favorable regimes (Scenario 1), the experimental chances could be spectacular, especially in the case of relatively light Higgs boson mass spectrum and a large (negative) value of λ 5 , for which the predicted one-loop cross sections σ (0+1) (e + e − → h 0 A 0 /H 0 A 0 ) at √ s = 500 GeV may border the O(100 fb) level. This would translate into barely 5 × 10 4 events per 500 fb −1 of integrated luminosity at the linac. While the former situation describes perhaps the most optimistic possibility, cross sections of O(10 fb) at √ s = 500 GeV should be quite common in a wide patch of the 2HDM parameter space, and correspondingly of O(1 fb) at √ s = 1 TeV (see Table IV), thus delivering rates that range from a few hundred to a few thousand events respectively. With this level of statistics (plus a comparable number of events from the H 0 A 0 channel) it should be perfectly possible to insure a comfortable tagging of the Higgs bosons in the clean environment of a linac machine.
Most important, we have been able to demonstrate that the dominant quantum effects are, as we expected, primarily nurtured by the Higgs-mediated one-loop corrections to the A 0 h 0 Z 0 /A 0 H 0 Z 0 vertices, and can be ultimately traced back to the potentially enhanced 3H selfinteractions λ 3H (cf. Table II). This is indeed a trademark feature of the 2HDM, with no counterpart in the MSSM -in which SUSY invariance highly curtails the structure of the Higgs boson self-interactions and compels them to be purely gauge-like. Obviously, probing the structure of Higgs boson self-interactions -in the present case through the analysis of radiative corrections on direct production processes -is a most useful strategy to disentangle SUSY and non-SUSY Higgs physics scenarios.
To be specific, for the analysis of the e + e − → 2H cross sections we have faced the computation of the full set of O(α 3 ew ) quantum effects, among them the subset of O(α 2 ew λ 2 3H ) corrections, and supplemented them with the leading O(α 4 ew ) pieces arising from the square of the one-loop diagrams involving Higgs boson self-couplings. These (finer) corrections are of O(α 2 ew λ 4 3H ) (see, however, below). Their inclusion is actually a must at large center-of-mass energies (i.e. once √ s has left behind the startup situation, and of course also in the 1 TeV regime). The reason is that the tree-level O(α 2 ew ) and the oneloop O(α 3 ew ) amplitudes virtually cancel out each other at high energies, owing to the characteristic large and negative quantum effects in this regime. We can convince ourselves of this fact from the severe depletion exhibited by the one-loop curves -see e.g. Figs. 15-17 -as compared to the tree-level ones. Typically, the effect appears after we increase the center-of-mass energy by 100 GeV beyond the initial √ s = 500 GeV, and persists till 1 TeV and afterwards. At large energies, the inclusion of the O(α 2 ew λ 4 3H ) effects is a consistency requirement as it prevents the cross-section from becoming negative. These terms can actually be relevant also at the initial energy interval √ s = 500 − 600 GeV at which the linac will first operate; for example, for the scenarios with maximally enhanced 3H self-couplings, they may lead to an additional ∼ 20 % contribution to δ r , hence pushing the overall quantum correction even higher. In this respect, we note that our earlier estimate of these higher order terms, Eq. (82), was able to hit the right order of magnitude, although it actually underestimates the maximum size revealed by the exact numerical analysis -not too surprising from such a rough attempt at guessing their bulk size.
Let us clarify that the two-loop diagrams involving the Higgs boson-mediated A 0 h 0 Z 0 /A 0 H 0 Z 0 vertex corrections, lead of course also to amplitudes of order O(α 2 ew λ 4 3H ). These are actually finite, and originate from the interference of the tree-level amplitude with the twoloop one. An obvious concern is then the following: can we safely neglect them? Upon inspection of the Higgsmediated two-loop diagrams, and taking into account power counting and dynamical considerations, one can show that they are actually suppressed, and particularly so for those scenarios where σ (0+1) (e + e − → 2H) is optimized. But, most significantly, there is a deeper observation to be adduced at this point, which is well in the spirit of the effective field theory approach to QFT. The (finite) two-loop effects induced by the enhanced Higgs boson self-couplings can be conveniently reabsorbed in the form of one-loop corrections to the bare Higgs selfcouplings, which first appear at the one-loop level in the process under consideration. We can iterate this algorithm at any order in perturbation theory by further reabsorbing the higher order Higgs mediated loops into the triple and quartic Higgs self-couplings introduced at oneloop order. In this way we can define a collection of effective 3H and 4H self-couplings that encapsulate all these higher order effects.
The above remark is an interesting one, as it tells us that the application of the stringent tree-level unitarity relations that we have used can be directly reiterated for these effective couplings. Therefore, no matter how big are the quantum corrections induced by the 3H selfcouplings in subsequent orders of perturbation theory, their overall effect is constrained by the same set of formal unitarity relations. If expressed in terms of the original couplings, these relations will be equally stringent so that, on balance, the maximum enhancement capabilities of the 3H self-couplings are basically the same at any order of perturbation theory. This is the procedure adopted in the present work, and we believe it is perfectly consistent, in the sense that we have retained the contributions to the e + e − → 2H amplitudes at leading order in the 3H self-couplings while employing the unitarity constraints of Ref. [67]- [69], also derived at leading order in λ 3H .
How do our results compare with the expectations in the MSSM? As we have stressed in section II, a lot of work on Higgs boson production in e + e − colliders has been reported in the literature thus far, but mainly within the context of the MSSM. In particular, supersymmetric radiative corrections to σ(e + e − → h 0 A 0 ) have been considered in detail, although mostly in the context of LEP and including sometimes a more or less timid incursion into the TeV-class colliders, cf. Refs. [25,26,35,37,38,40]. Let us try to see how some of these results compare with ours. For the sake of definiteness, let us concentrate on the particular case of Ref. [37], in which σ(e + e − → A 0 h 0 ) is analyzed in the MSSM for a linear collider operating at LEP 2 energy and further extended up to √ s = 500 GeV. In the most favorable regimes, the loop-corrected cross sections may reach up to barely 30 fb at √ s = 500 GeV, which falls in the ballpark of the results that we have obtained here within the 2HDM. This is not surprising because in both cases the leading amplitude (i.e. the lowest order one) is purely gauge and hence the order of magnitude of the cross-section is independent of whether we compute it in the MSSM or in a generic 2HDM.
Quantum effects can be fairly efficient in the MSSM too, and they may also boost the tree-level predictions substantially (∼ 20%). However, their size is never as big as the maximal 2HDM contributions; and, no less important for our conceptual understanding, the fundamental origin of these quantum effects is completely different from that of the 2HDM. As already mentioned in the introduction, the leading quantum effects in the MSSM are due to the genuine Higgs-quark, Higgs-squark and quark-squark-chargino/neutralino Yukawa-like supersymmetric couplings, whose contributions carry enhancement factors of the fashion ∼ m b tan β/M W , ∼ m t (A t − µ/ tan β)/M 2 SUSY , m b (A b − µ tan β)/M 2 SUSY , none of them related to the structure of the Higgs potential and hence none of them associated to Higgs boson self-couplings.
In contradistinction to the MSSM, the 3H selfcouplings in the 2HDM embody the full potential for triggering quantum effects in physical processes, and these effects become thence tied to the very structure of the Higgs potential. Even if quantitatively similar in some scenarios, quantum effects in SUSY and non-SUSY two-Higgs doublet extensions of the SM are, therefore, prompted by intrinsically different mechanisms. In practice, the task of distinguishing whether the produced Higgs boson particles at a linac are supersymmetric or not can start by addressing the decaying signatures of the Higgs bosons into bb or τ + τ − pairs, as sketched above. However, in certain regions of parameter space the number of events can be by itself a direct signature of the sort of Higgs model that we have behind. For example, in the MSSM, the tree-level coupling C A 0 h 0 Z 0 ∼ cos(β − α) undergoes a severe suppression with growing values of tan β and/or M A 0 -as α depends on M A 0 in the MSSM. Thereby the predicted cross sections, either at the Born or at the one-loop level, are dramatically weakened as a function of M A 0 . For instance, σ (0+1) (e + e − → h 0 A 0 ) at √ s = 500 GeV and M A 0 = 300 GeV may end up as tiny as ∼ 0.01 fb (see e.g. Fig. 9 of Ref. [37]) while the 2HDM prediction is of O(10 fb), i.e. one thousand times bigger! But there is more to say here, as we have announced above. The identification strategy can be highly facilitated through the interplay of further 2H (and eventually 3H) channels, such as (1) and (2), which are also highly distinctive in the 2HDM -see [27] and [29]. In fact, for these channels the difference with the MSSM can be apparent already at the lowest order of perturbation theory.
Let us be a bit more quantitative here. In Table X, we compare the numerical predictions on the cross-sections of these various processes in truly equitable conditions, i.e. for different sets of common parameters for all the processes, and for three realistic values of the center-ofmass energy of the ILC/CLIC (the latter operating always at the highest edge of the energy band, although it could nominally reach the 3 TeV end).
To be more precise, the cross-sections that we are comparing on equal footing in Table X are the following: i) the one-loop corrected σ (0+1) (e + e − → A 0 h 0 ); ii) the leading-order O(α 3 ew ) cross-sections of the two triple Higgs production processes e + e − → h 0 H 0 A 0 /H 0 H + H − ; and iii) σ(e + e − → V * V * → h 0 h 0 + X) at leadingorder O(α 4 ew ). Notice that the bulk of the contribution from processes ii) and iii) comes from the O(α 2 ew λ 2 3H ) and O(α 3 ew λ 2 3H ) parts, respectively. In all cases we take tan β = 1, α = β and the maximum allowed value of |λ 5 | according to Table VII, that is, the most favorable scenario for h 0 A 0 pair production. We use Sets I-III for the Higgs boson masses. The resulting predictions for the different channels turn out to be highly illustrative of the complementary nature of such Higgs-boson production signatures. For a thorough discussion of the dynamical features and signatures of the multiparticle processes ii) and iii), we refer the reader to the original references [27,29]. Here we have recomputed the cross-sections ii) and iii) under the same kinematical conditions and in the very same region of the parameter space as the main processes i) under consideration, which is of course the indispensable first step that enables us to compare them all meaningfully. In this way, we are well prepared to infer the final strategy that can be designed to efficiently search for 2HDM Higgs bosons at the linear colliders.
The strategy that follows from the results encapsulated in Table X and the information that we have gathered in the figures of the last section should be by now crystalclear, and it comes down to the following two-step procedure: • Step 1. At the vicinity of the startup energy of the ILC (approximately in the range √ s = 500 ± 100 GeV), the exclusive neutral double Higgs boson channels e + e − → A 0 h 0 /A 0 H 0 -Eq. (3)prove to be the dominant ones as compared to the triple Higgs boson production processes e + e − → 3H and the inclusive double Higgs production ones e + e − → 2H + X, see equations (1) and (2). Therefore, in this lower energy band, the study of potential signatures of new physics, and in particular the identification of the nature of the produced Higgs bosons at the ILC, must be conducted through the careful analysis of the quantum corrections affecting the two-body channels e + e − → 2H. Here is precisely where the detailed results of the present work could be most useful, and indeed should be the first analysis to be implemented when the lin-Process σ( √ s = 0.5 TeV) fb σ(1.0 TeV) fb σ(1.   (1), (2) and (3). The results are obtained for the Sets I, II and III of Higgs boson masses (cf. Table IV) tan β = 1, α = β, and three different values of the center-of-mass energy. We observe a great complementarity between the different channels at different energies: the exclusive 2H channels (3) are dominant at the ILC startup energy √ s = 0.5 TeV whereas the others dominate at higher energies, specially the inclusive 2H +X one (1)-which is triggered mainly by weak gauge boson fusion.
ear colliders are set to work in the future; • Step 2. At higher energies, however, say for √ s > 600 GeV and, for that matter, in the entire TeV range (hence at the highest nominal energy regimes scheduled for the ILC, and specially for CLIC), the influence of the Z 0 -propagator stifles dramatically the cross-section of the exclusive neutral double Higgs boson channels e + e − → 2H. Therefore, in these high energy domains, they can no longer compete with the mechanisms providing multiparticle final states, such as (1) and (2). The latter process, e + e − → V * V * → h h + X, becomes indeed the most efficient one at the highest energies since it is not crippled by the s-channel propagator; quite on the contrary, its cross-section increases steadily with the energy (see Ref. [29] for details). Therefore, in the upgraded phase of the ILC, and certainly for the CLIC collider, one must concentrate all the search strategy power on looking for an anomalously large number of inclusive Higgs boson pairs of the type h 0 h 0 , H 0 H 0 , h 0 A 0 and A 0 A 0 , which should emerge mostly acollinear and dynamically thrusted along the beam direction, rather than appearing in a simple back-to-back geometry characteristic of the two-body process e + e − → 2H. Notice furthermore that, in contradistinction to the latter, the fused pairs consist of identical Higgs bosons. Unmistakeably, the predicted numbers in Table X could not, by any means, be accounted for if they were to be ascribed to a supersymmetric origin. Finally, the cross-correlation of these higher energy effects with the previous ones at lower energies, which we could track very well while upgrading the ILC collider from √ s = 0.5 TeV up to √ s = 1.5 TeV -and eventually till √ s = 3 TeV (through CLIC)should provide plenty of unambiguous evidence of non-supersymmetric Higgs physics 14 .
Let us finally introduce a remark that hints at the potentially far reaching implications of the 2HDM dynamics in different sectors of Particle Physics. The triple Higgs boson self-couplings could also have an indirect impact on the best high precision observables at our disposal, even well before the ILC/CLIC colliders are commissioned and put effectively to work. For example, they could have a bearing on the famous electroweak precision parameter ∆r -see Eq. (11). Since our main aim here has been to exploit the leading contributions from the Higgs boson self-couplings in the arena of the direct production processes (3), the ∆r part did not enter our quantum computation; it would only enter at two-loop level. However, we suspect [90] that the influence of the 3H self-couplings on ∆r might have a real interest per se, for it could induce a correction to δρ and ultimately trigger a shift in the W ± mass, and the result might well compete with the highly accurate calculations that have been performed on this observable within the alternative framework of the MSSM [14].
To summarize, if the genuine enhancement properties of the 2HDM effectively hold in the real world, the combined analysis of the exclusive double Higgs production (2H) processes at the startup energy of the ILC, √ s = 500 GeV, and subsequently of the triple (3H) and inclusive double Higgs production processes (2H + X) at the upgraded ILC/CLIC ( √ s 1 TeV), might reveal strong hints of Higgs boson physics beyond the SM. The bottom line of our study is that the physics of the linear colliders can be truly instrumental to unveil the nature of the Higgs bosons. These bosons may, or may not, have been discovered during the LHC era while the linear colliders take the floor, but the real issue at stake here is of highest priority and should strengthen the need for such machines. High precision Higgs boson physics in a linac can indeed provide a keen insight into the most sensitive building blocks of modern gauge theories of weak and electromagnetic interactions, viz. in the very core architecture of the Electroweak Symmetry Breaking.