Quantum effects on Higgs-strahlung events at Linear Colliders within the general 2HDM

The associated production of neutral Higgs bosons with the Z gauge boson is investigated in the context of the future linear colliders, such as the ILC and CLIC, within the general two-Higgs-doublet model (2HDM). We compute the corresponding production cross-sections at one-loop, in full consistency with the available theoretical and phenomenological constraints. We find that the wave-function renormalization corrections to the external Higgs fields are the dominant source of the quantum effects, which turn out to be large and negative, and located predominantly in the region around \tan\beta=1 and moderate values of the parameter \lambda_5 (being \lambda_5<0). This behavior can be ultimately traced back to the enhancement potential of the triple Higgs boson self-couplings, a trademark feature of the 2HDM with no counterpart in the Higgs sector of the Minimal Supersymmetric Standard Model. The predicted Higgs-strahlung rates comfortably reach a few tens of femtobarn, which means barely 10^3 - 10^4 events per 500 inverse femtobarn of integrated luminosity. Due to their great complementarity, we argue that the combined analysis of the Higgs-strahlung events and the previously computed one-loop Higgs-pair production processes could be instrumental to probe the structure of the Higgs sector at future linac facilities.


I. INTRODUCTION
After more than 40 years since the seminal ideas coined by a handful of theoretical pioneers [1], our present understanding of the Electroweak Symmetry Breaking (EWSB) phenomenon through the Higgs (Englert-Brout and Guralnik-Hagen-Kibble) mechanism is still rather incomplete and experimentally inconclusive. On the one hand, there is no compelling alternative to consistently embed the EWSB mechanism into the quantum field theoretical description of particle physics offered by the -in so many regards successful -Standard Model (SM) of the Strong and Electroweak interactions. On the other hand, we do not have a single phenomenological hint of the existence of elementary scalar fields, not to mention the fact that we do not understand how to make compatible the EWSB mechanism and its associated vacuum energy with fundamental problems of different scope, for example in the domain of cosmology. Still, the possibility to describe the inner theoretical structure of the SM through the EWSB mechanism is so successful within the restricted particle physics domain that it would not * Electronic address: dlopez@ecm.ub.es † Electronic address: sola@ecm.ub.es ‡ Electronic address: nicolas.bernal@cftp.ist.utl.pt be wise, not even advisable, to cease our pursue of the phenomenological implications of the Higgs mechanism paradigm till the limits of the current experimental possibilities.
It goes without saying that the quest for experimental evidences of the Higgs boson is a most preeminent milestone of the upcoming generation of collider facilities. Nonetheless, the future data might well reveal that the purportedly found Higgs boson actually belongs to a richer model structure, which might be grounded somewhere beyond the minimal conception of the SM, namely of a single, spinless, fundamental constituent of matter. If so, a few more fundamental spinless constituents could appear. A particularly well-motivated extension is the two-Higgs-doublet model structure encompassed by the Minimal Supersymmetric Standard Model (MSSM) [2]. Here the physical Higgs boson spectrum contains a couple of CP-even (h 0 , H 0 ), one CP-odd (A 0 ) and two charged Higgs bosons (H ± ), and the corresponding Higgs potential is highly constrained by the underlying Supersymmetry (SUSY). One particular consequence of the latter is that SUSY invariance restricts all Higgs boson selfinteractions to be gauge-like. While this represents an economy from the point of view of the number of couplings in the potential, it is strongly imbalanced by the exceeding number of parameters (mixings and masses) populating the other sectors of the MSSM. In addition, the pure gauge nature of the Higgs boson self-couplings make them highly inconspicuous from the practical point of view, in the sense that they are unable to trigger any outstanding phenomenological signature. The core of the enhancement capabilities of the MSSM Lagrangian resides, instead, in the multifarious pattern of Yukawa couplings between the Higgs bosons and the quarks, as well as between quarks, squarks and charginos/neutralinos. The rich interplay of opportunities that they give rise to has been extensively analyzed in the past within a plethora of processes, see e.g. [3][4][5][6][7], and also [8][9][10][11] for reviews on the subject.
The two-Higgs SU L (2)-doublet structure of the Higgs sector in the MSSM constitutes a genuine prediction of the SUSY dynamics. Nonetheless, a more general architecture can appear in the form of a non-SUSY framework through the so-called general (unconstrained) two-Higgs-doublet Model (2HDM) 1 . Although they share a common Higgs boson spectrum, the potentially most distinctive phenomenological features of both models are located in very different sectors. The most relevant observation here is that, in the absence of an underlying SUSY, the triple (3H) and quartic (4H) Higgs-boson selfinteractions are no longer restrained to be purely gauge. This can have a tremendous impact e.g. in the physics of the top quark in hadron colliders, as it was shown long ago in [12], and it can also trigger significant neutral flavor-changing interactions [13] that can perfectly compete with the SUSY ones [14,15].
Much attention has been devoted to Higgs boson production and decay in hadron colliders, extending from the ongoing Tevatron facility at Fermilab to the brand new LHC collider recently operating at CERN [8][9][10][11][12][13][14][15][16][17]. However, precision Higgs boson physics will greatly benefit from the interplay [18] of the forthcoming generation of linear colliders (linac), such as the ILC and CLIC projects [19]. Here a variety of processes can provide new clues to Higgs boson physics, e.g. the production of triple Higgs-boson final states, both in the MSSM [20] and in the general 2HDM [21]; the double Higgs-strahlung channels hhZ 0 [22]; and the inclusive Higgs-pair production via gauge-boson fusion [20,23]. In the same vein, also the γγ mode of a linac has been explored [24], in particular the loop-induced production of a single neutral Higgs boson [25] and of a Higgs boson pair [26]. In all the above mentioned cases, promising signatures were singled out and illustrate that, if effectively realized in Nature, such hints of a generic 2HDM structure could hardly be missed in the superbly clean environment of a linac. And, what is more, they could not be confused as having a SUSY 1 To be more precise, we refer to this model as "general" in the sense that we allow all possible operators leading to a renormalizable, gauge-invariant and CP-conserving Higgs potential. As usually done in practice [11], we impose an additional (softlybroken) Z 2 discrete symmetry Φ i → (−1) i Φ i , where i = 1, 2 denote each of the SU L (2) Higgs doublets, as a sufficient condition to banish the tree-level flavor changing neutral currents.
origin, because of the intrinsically different nature of the Higgs self-interaction sector. Besides, outstanding fingerprints of a generic 2HDM could also be stamped in the pattern of radiative corrections to Higgs production processes; for instance, quantum effects on the cross-sections for two-body Higgs boson final states: have been attentively investigated in the MSSM [27][28][29][30].
As for the general 2HDM, the efforts were first concentrated on the production of charged Higgs pairs [31]. This program has been recently brought to completion from a full-fledged study of the quantum effects on the production cross-sections in the neutral Higgs sector [32]. Alongside these processes we find the more traditional Higgs-strahlung events, in which a Higgs boson is produced in association with the Z 0 : (Notice that the A 0 Z 0 final state is forbidden by CPconservation.) Processes (2) are complementary to the (1) ones -cf. [33], and also [20,34] for a phenomenological analysis in the MSSM context. Let us also recall that the last available limit on the SM Higgs boson mass, placed by LEP searches, comes precisely from investigating the "Bjorken process" [35], i.e. e + e − → H Z 0 followed by Z 0 → ff , with the result: M H 114.4 GeV [36]. In this article, we aim at discussing the corresponding one-loop corrections to the generalized Bjorken processes (2) within the 2HDM. However, in contrast to the original SM case, where the radiative corrections are small, the quantum effects on the processes (2) can be large and are mainly driven by the 3H self-couplings. In fact, our main aim here is to identify the regions of the parameter space where this is so, and then quantify the impact of the potentially enhanced 3H self-couplings on the final cross-sections. In this way, joining this study with that of Ref. [32] on the pairwise Higgs boson production channels (1) in the 2HDM, a rather complete panorama of the genuine quantum effects associated to the basic neutral Higgs boson production channels becomes available.

II. HIGGS-STRAHLUNG EVENTS AT ONE-LOOP: THEORETICAL SETUP
A good deal of attention has been devoted in the literature to multiple properties of the Higgs-strahlung process e + e − → HZ 0 in the SM, see e.g. [37] and references therein. Not surprisingly it was one of the gold-plated channels for Higgs boson search at LEP. Concerning the MSSM, the Higgs-strahlung channels (2) have also been discussed extensively -cf. Refs. [29,38] and part two of the review [9] -and are currently under investigation also in the MSSM with CP-violating phases [39] (the so-called complex MSSM [40]). The upshot of these analyses spotlights the following features: i) the dominant source of corrections at one-loop originate from the Higgs boson propagators, which can be reabsorbed into an effective (loop-corrected) mixing angle α eff (equivalently, a more generic mixing matrix in the complex MSSM case); ii) the corrections to the ZZ h vertex are in general small, although they can reach the level of 10 % for very low (or high) values of tan β, precisely in the regions where the Higgs Yukawa couplings to heavy quarks become enhanced; iii) the electromagnetic corrections to the initial state with virtual photonic corrections and initial-state radiation do not differ from the SM case, being in general large and positive (except near the production threshold). In the present article, our endeavor is to seek for the genuine phenomenological imprints associated to the generic (unconstrained) 2HDM dynamics, most particularly to the potentially enhanced 3H self-couplings.
Let us recall that the general 2HDM [11] is obtained by canonically extending the SM Higgs sector with a second SU L (2) doublet with weak hypercharge Y = +1, so that it contains 4 complex scalar fields. The free parameters λ i in the general, CP-conserving, 2HDM potential can be finally expressed in terms of the masses of the physical Higgs particles (M h 0 , M H 0 , M A 0 , M H ± ), tan β (the ratio of the two VEV's H 0 i giving masses to the up-and down-like quarks), the mixing angle α between the two CP-even states and, last but not least, the selfcoupling λ 5 , which cannot be absorbed in the previous quantities. Therefore we end up with a 7-free parameter set, to wit: . Furthermore, to ensure the absence of tree-level flavor changing neutral currents (FCNC), two main 2HDM scenarios arise: 1) type-I 2HDM, in which one Higgs doublet couples to all quarks, whereas the other doublet does not couple to them at all; 2) type-II 2HDM, where one doublet couples only to down-like quarks and the other doublet to up-like quarks. The MSSM Higgs sector is actually a type-II one, but of a very restricted sort (enforced by SUSY invariance) [11]. We refer the reader to Ref. [32] for a comprehensive account on the structure of the model and for notational details. In particular, Table II of that reference includes the full list of trilinear couplings within the general 2HDM that are relevant for the present calculation. Further constraints must be imposed to assess that the SM behavior is sufficiently well reproduced up to the energies explored so far. Most particularly, we take into account i) the approximate SU (2) custodial symmetry, which can be reshuffled into the condition |δρ 2HDM | ≤ 10 −3 [41]; and ii) the agreement with the low-energy radiative B-meson decay data (which demand M H ± 300 GeV for tan β ≥ 1 [42] in the case of type-II 2HDM).
Additional requirements ensue from the theoretical consistency of the model, to wit i) perturbativity [43]; ii) unitarity [44]; and iii) stability of the 2HDM vacuum [45]. Let us expand a bit more on the latter conditions, which turn out to play a crucial role in our analysis. As for the unitarity bounds, and following Ref. [44], the basic underlying strategy is to compute the S-matrix elements S ij for the possible 2 → 2 processes involving Higgs and Goldstone bosons in the 2HDM, and to subsequently restrain the corresponding eigenvalues U ik S kl U −1 lj = α i δ ij by the generic condition |α i | < 1/2 ∀i. The latter requirement translates into a set of upper bounds on the quartic scalar couplings λ i , i = 1 . . . 6 -and hence on combinations of Higgs masses, trigonometric couplings and, most significantly, the parameter λ 5 . As far as vacuum stability restrictions are concerned [45], they may be written as follows: More refined versions of these equations may be obtained from the Renormalization Group (RG) running of the quartic couplings λ i (µ 2 ) at high energies. Nonetheless, for our current purposes we need not assume any particular UV completion of the 2HDM. This would be an unnecessary additional assumption at this stage of the phenomenological analysis of our processes. Therefore, we are not tied to any specific UV-cutoff (which could be, in principle, as low as Λ ∼ 1 − 10 TeV). In this sense, we may just apply Eq. (3) with all λ i taken at the EW scale. This is after all the scale at which we fix all our input parameters and perform the phenomenological analysis (including the renormalization) of the processes under study. By the same token, the (cutoff dependent) triviality bounds are also circumvented. The latter kind of restrictions only impose tight upper limits on the Higgs boson masses when a very large UV-cutoff (e.g. the Planck scale or the GUT scale) is considered. Our choices of Higgs boson mass spectra (cf. Table I) lie, in any case, in a sufficiently low mass range so as to conform with the typical mass requirements allowed by triviality [46] -even for large cutoff scales. More details on the constraints setup are provided in Ref. [32].
The dynamics of the Higgs-strahlung processes under consideration is driven at leading order by the tree-level interaction Lagrangians: where s W ≡ sin θ W and c W ≡ cos θ W for the electroweak mixing angle θ W . Since  events is most likely insufficient to disclose their true nature. Similarly, in the limit α = β − π/2 the h 0 Z 0 Z 0 coupling coincides with the analogue coupling in the SM, HZ 0 Z 0 . It is, therefore, the pattern of radiative corrections the characteristic signature associated to each one of the possible models; most particularly, it should help to disentangle SUSY versus non-SUSY extended Higgs physics scenarios.
The leading-order O(α ew ) scattering amplitude follows from the s-channel Z 0 -boson exchange and renders Here p 1,2 and η 1,2 refer to the 4-momenta and helicities of the electron and positron, and ǫ (k 2 , σ 2 ) is the polarization 4-vector of the Z gauge boson with 4-momentum k 2 and helicity σ 2 . We have introduced also 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 notice that we do not include the finite Z 0 -width corrections, since they are completely negligible for the center-of-mass energies that we consider here. Finally, the total cross section σ(e + e − → h 0 Z 0 ) at the tree-level is obtained after squaring the matrix element (5), performing an averaged sum over the polarizations of the colliding e + e − beams and the outflowing Z 0 boson, and integrating over the scattering angle.
The calculation of σ(e + e − → h 0 Z 0 ) at one loop is certainly much more cumbersome. To start with, it is UV-divergent and it thus requires of a careful renormalization procedure in order to render finite results. We adopt here the conventional on-shell scheme in the Feynman gauge [47]. In the MSSM, these cross-section calculations can be mostly carried out in an automatized fashion through standard algebraic packages which allow a rapid and efficient analysis of the electroweak precision observables [48]. Several public codes are available, see e.g. [49][50][51][52]. However, in our case the calculation is non-supersymmetric and we must deal with the renormalization of the Higgs sector in the class of generic 2HDM models. A detailed description of the renormalization procedure for the 2HDM Higgs sector in the onshell scheme has been presented in [32] and we refer the reader to this reference for all the necessary details. We will make constant use of the framework described in this reference and sometimes we will refer to particular formulae of it.
With this renormalization setup in mind, the vari-  Table II of [32]. In such cases, we shall include these factors when assessing the order of magnitude of the diagram. We are now ready for sorting out the one-loop contributions in different categories: • Self-energy corrections to the Z 0 and γ − Z 0 mixing propagators, all of them of order O(α 2 ew ) with no enhancement factors at this order.
• Vertex corrections to the e + e − Z 0 interaction, which are also of order O(α 2 ew ). It should be noted that we do not include the virtual photonic O(α em α ew ) effects nor the real bremsstrahlung emission off the e + /e − legs. These pure QED corrections and the 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 focus. For this kind of processes involving electrically neutral Higgs bosons in the final state, the one-loop QED effects are confined to the initial e + e − vertex. In practice, the net outcome of the accompanying initial state radiation (ISR) is to lower the effective center-of-mass energy available for the annihilation process. All in all, these effects are unessential at this stage to test the presence of the new dynamical features triggered by the 2HDM in the Higgsstrahlung events under analysis.
• Vertex corrections to the e + e − H interaction. Since we have explicitly set m e = 0 throughout our calculation, the e + e − H tree-level Yukawa coupling is absent and the corresponding one-loop diagrams automatically render a UV finite contribution which is, in any case, very small.
• The loop-induced γZ 0 h 0 interaction. This one is order O(e α ew λ 3H ) and therefore includes an enhancement factor λ 3H . Due to the γ − Z 0 mixing at one-loop, the following counterterm is needed so as to render a UV-finite vertex: and similarly for the counterterm associated to the effective L H 0 Z 0 γ interaction: • The vertex correction for hZ 0 Z 0 encompasses different contributions of order O(α 2 ew , α ew e λ 3H ) (see the sample diagrams in Fig. 3). And, most important, there are, in addition, corrections of O(α ew λ 2 3H ) (in fact, the dominant ones) which are characterized by the highly conspicuous enhancement factor squared λ 2 3H . The latter originates from the associated vertex counterterm, most particularly from the Higgs field renormalization constant Z 1/2 As we shall discuss in more detail below, δ Z h 0 is sensitive to the scalarscalar self-energy and hence it involves products of two triple Higgs self-couplings (see Fig. 4). Of course, an equivalent discussion holds for the complementary channel, H 0 Z 0 . The full form of the associated one-loop vertex counterterms ensues from the usual splitting of bare fields and parameters into the renormalized ones and associated counterterms, to wit: For the h 0 Z 0 Z 0 and H 0 Z 0 Z 0 case one gets respec-tively: where we can spot in the structure of these formulae the presence of the above mentioned Higgs field renormalization counterterms δ Z h 0 and δ • The finite wave-function (WF) correction to the external Higgs fields h 0 or H 0 (Fig. 4). These are related to the fact that we have chosen the residue of the A 0 -propagator at its pole to be one, and therefore there is no more freedom to make the same choice for the other Higgs bosons. This entails a finite WF renormalization correction, see Ref. [32] for details. These finite renormalization effects are of utmost importance in this case, as they trigger leading contributions of order O(α ew λ 2 3H ). They actually drive the very bulk of the quantum effects (see below for a more detailed discussion).
• Finally, the box-type O(α 2 ew ) diagrams ( Figure 2), whose contribution is non-negligible, in particular at large center-of-mass energies -since they are not suppressed as 1/s. The dominance of the WF corrections is a very characteristic feature of the processes under analysis. It originates from a non-trivial cancelation between the Higgs boson field counterterms, which appear in two different pieces of the overall one-loop amplitude, though with opposite signs. Let us further elaborate on this important point. Without loss of generality, we concentrate on the h 0 Z 0 channel. The complete set of contributions to the scattering amplitude up to the one-loop level may be split in the following manner: wherein the finite WF corrections to the external Higgs field are introduced as At this point we are making explicit use of the on-shell renormalization conditions defined in [32], in particular ℜeΣ h 0 H 0 (M 2 h 0 ) = 0. Moreover, we only retain those contributions that are leading-order in the triple Higgs self-couplings. Let us notice that the counterterm piece δ M (1) is sensitive to O(λ 2 3H ) effects through the WF renormalization of the Higgs fields. More specifically, from the explicit expression of the h 0 -field counterterm as a linear combination of the two Higgs doublet counterterms (see Section IV of Ref. [32] for details) we may single out the following O(λ 2 3H ) contribution from Eq. (9), But, as warned, the finite WF renormalization factor of the h 0 field in (12) is sensitive to O(λ 2 3H ) terms too: where the dots stand for those terms with dependencies other than O(λ 2 3H ). Notice, therefore, that part of the O(λ 2 3H ) dependence cancels out between the counterterm diagram associated to the h 0 Z 0 Z 0 vertex -cf. Eq. (14) -and the finite h 0 WF-factor -cf. Eqs. (12) and (15). Specifically, it is the piece ℜe Σ ′ A 0 (M 2 A 0 ) that exactly cancels between the two terms. The remaining O(λ 2 3H ) contributions come from ℜe Σ ′ h 0 (M 2 h 0 ) and are generated by the Feynman diagrams displayed in Fig. 4. When sufficiently enhanced, these pieces account for the bulk of the one-loop quantum corrections. A typical diagram in the class of 2-point functions provides the following contribution 2 : wherein the scalar two-point function is defined as in Ref. [49]: µ denoting the 't Hooft mass unit. In the expression (16), M denotes the typical mass scale(s) which appear in these one-loop 2-point functions. Owing to the derivative with respect to the external momentum, the 1-point functions do not contribute to the WF renormalization, and hence the quartic Higgs boson self-couplings are not involved in this calculation. This means that the first line of diagrams in Fig. 4 does not actually contribute. We emphasize that the overall sign for this expression depends on the sign of such 2-point function. Notice that 2 Notice that we employ the convention that iΣ equals the Higgs boson self-energy diagram, such that Σ does not contain the global imaginary part emerging from the loop integral. This is why we have multiplied the loop integral in (16) by −i before taking the real part.
ℜ eB ′ 0 (p 2 , M 2 , M 2 ) > 0 for p 2 < 4M 2 , which is the case we wish to focus on for the h 0 Z 0 channel, and also for H 0 Z 0 (as long as the resonant decay H 0 → h 0 h 0 is forbidden by kinematics). Therefore, in most of the scenarios of interest, such finite WF corrections carry an overall minus sign. Being proportional to λ 2 3H , they become the leading quantum effects in the region where the trilinear couplings are enhanced, and as a result they generate a characteristic pattern of quantum effects in which a systematic suppression of the tree-level cross-section is predicted.
The presence of the large negative corrections induced by the Higgs boson self-energies brings forward a characteristic signature for the production cross-sections of the Higgs-strahlung processes (2). This feature is in marked contradistinction to the situation with the Higgs boson pair production mechanisms(1), where the corresponding corrections are just opposite in sign, i.e. large and positive, see [32]. We shall further comment on these interesting and correlated features in the next section. From the general structure of Eq. (16), one may anticipate the typical (maximum) size of the quantum effects on the Higgs-strahlung processes as follows: f being a dimensionless form factor (basically accounting for the behavior of the B ′ 0 function). The notation . . . stands for the various operations of averaging and integration of the squared matrix elements. Taking into account that unitarity limits let the trilinear couplings reach values as large as λ3H MW ≃ |λ5| e ≃ 30, and assuming M ∼ 200 GeV and f ∼ O(1), the above estimate typically predicts a strong depletion of the tree-level signal by δ r ≃ −90%. This prediction falls in the right ballpark of the exact numerical results that will be reported in the next section, which point to a maximum depletion of δ r ∼ −60%.
Before closing this section, let us recall that the counterterm amplitude δ M (1) e + e − →h 0 Z 0 derives from a number of renormalization conditions that determine the renormalized coupling constants and fields in a given renormalization framework. As we have said, in our calculation the renormalization is performed in the conventional on-shell scheme in the Feynman gauge, appropriately extended to include the 2HDM Higgs sector. For the latter, we need in particular a renormalization condition for the parameter tan β = v 2 /v 1 . We adopt the following [53]: This condition insures that the ratio v 2 /v 1 is always expressed in terms of the true vacua after the renormalization of the Higgs potential. The corresponding counterterm resulting from tan β → tan β + δ tan β can then be computed explicitly: This counterterm is involved in Eqs. (9)-(10), and it also determines the WF mixing term δ Z h 0 H 0 that appears in these equations, as follows: δ Z h 0 H 0 = sin 2α(δ tan β/ tan β). We refer the reader once more to the exhaustive presentation of Ref. [32] for the renormalization details.

III. NUMERICAL ANALYSIS
In this section, we present the numerical analysis of the one-loop computation of the Higgs-strahlung processes (2) within the 2HDM. We shall be concerned basically with the following two quantities: i) the predicted crosssection at the Born-level σ (0) and at one loop σ (0+1) ; and ii) the relative size of the one-loop quantum corrections: We have carried out our analysis with the help of the standard algebraic and numerical tools FeynArts, Form-Calc and LoopTools [49]. The different Higgs boson mass sets that shall be used hereafter are quoted in Table I. Let us highlight that, due to their mass splittings (and also to the fact that M h 0 ≥ 140 GeV), Sets C and D can only be realized in a general (non-SUSY) 2HDM framework. Furthermore, in view of the value of the charged Higgs mass, Sets A-C are only suitable for type-I 2HDM, whereas Set D is valid for either type-I and type-II 2HDM's.
Apart from reflecting a variety of possible situations in the 2HDM parameter space, some of these sets can be mimicked by the Higgs boson mass spectrum in supersymmetric theories. For instance, Sets A and B can be ascribed to characteristic benchmark scenarios of Higgs boson mass spectra within the MSSM; in particular, Set B lies in the class of the so-called maximal mixing scenarios [54], for which M h 0 takes the highest possible values within the MSSM. The numerical mass values for the MSSM-like sets have been obtained with the aid of the program FeynHiggs by taking the full set of EW corrections at one-loop [50].
We remark that the more massive the Higgs bosons are, the stronger the constraints that unitarity imposes over |λ 5 |. The maximum (negative) values are roughly attained for λ 5 ≃ −9 (Set A); λ 5 ≃ −10 (Sets B and C); and λ 5 ≃ −8 (Set D).  Before coming to grips with the analysis of the Higgsstrahlung events in the general 2HDM, it is interesting to briefly reconsider the analogue process in the simpler framework of the SM. In In Figs. 6-9 we illustrate the fundamental phenomenological features associated to the process e + e − → h 0 Z 0 in the scope of the general 2HDM. Fig. 6 summarizes the pattern of radiative corrections δ r projected onto the (tan β, λ 5 ) and (sin α , tan β) planes, the former at fixed α = β−π/2 and the latter with λ 5 = −2. We remark that for α = β − π/2 the h 0 Z 0 Z 0 tree-level coupling takes on the SM form and therefore is maximal. For definiteness, these plots have been generated for a Higgs boson mass spectrum as in Set B (cf. Table I), and at the fiducial ILC startup center-of-mass energy, √ s = 500 GeV. 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. Fig. 6 illustrates that the allowed region in the (tan β, λ 5 ) plane is severely restrained by the theoretical constraints stemming from the perturbative unitarity and vacuum stability. On the one hand, λ 5 > 0 values are strongly disfavored by the vacuum stability condition; on the other hand, the unitarity constraints tend to disfavor moderate and large values of tan β (specially tan β values significantly larger than 1) as well as of tan β ≪ 1. Altogether these constraints set an approximate lower bound of λ 5 ∼ −10 for tan β ∼ 1 and a rigid upper bound excluding almost all positive values of λ 5 for any tan β. The combined set of constraints builds up a characteristic physical domain, with a valley-shaped area centered at tan β ∼ 1 that sinks into the λ 5 < 0 region and becomes narrower with growing |λ 5 |.
The fact that the curves of constant δ r do not depend on tan β (left panel of the Fig. 6) is related to the choice α = β −π/2, which we have made in order to consider the situation where the tree-level cross-section for h 0 Z 0 production is maximal. For this choice of the mixing angle α, it turns out that no particular enhancement shows up for any value of tan β, neither from the 3H self-couplings nor from the Higgs-top quark Yukawa couplings. Indeed, under the condition α = β − π/2, the h 0 tt Yukawa coupling does not depend on tan β: and at the same time the 3H self-couplings which are relevant for this channel become independent of tan β; notice, for example, that The final outcome is that the potential dependence of the computed observables on tan β vanishes as long as we stick to these α = β − π/2 configurations -in which  Table I. the tree-level coupling h 0 Z 0 Z 0 is maximum and formally equivalent to that of the SM 3 . As a result, the only feasible mechanism able to significantly enhance the quantum effects in these scenarios is by increasing the value of the |λ 5 | parameter (towards more negative values, so as to be consistent with vacuum stability). Once we depart from the α = β − π/2 setting, we recover of course the expected dependence of the computed observables with tan β (cf. right panel of Fig. 6), but then the lowest order h 0 Z 0 production cross section becomes smaller. For tan β < 1, radiative corrections are boosted as a result of the enhanced 3H self-couplings, and partially also due to the Higgs-top quark Yukawa couplings. Either way, their overall effect is to suppress the tree-level signal, as such leading corrections are triggered primordially by the finite Higgs boson WF corrections. Another source of enhancement of δ r appears near the region where β ∼ α. However, this effect is not caused by a real increment of the one-loop cross-section σ (0+1) , but by a mere suppression of the cross-section at the Born-level, i.e. σ (0) in the denominator of Eq. (21), and in this sense it is an uninteresting situation.
In Fig. 7 we explore the evolution of the cross section as a function of the center-of-mass energy. We include, in each plot, the tree-level contribution σ (0) and the corresponding loop-corrected value, σ (0+1) , for different values of λ 5 . The right panel tracks the related behavior of the quantum correction δ r as a function of √ s, for the same set of λ 5 values. The plots are generated for Set B of Higgs boson masses, at fixed tan β = 1 and α = β − π/2. The leading-order cross section σ (0) curve exhibits the expected behavior with √ s, as it scales with the s-channel Z 0 -boson propagator, namely proportional to 1/(s − M 2 Z ). A similar pattern is also followed by the full loop corrected cross-section.
The range where the relative one-loop correction δ r is positive is very reduced and it is confined to a regime where the center-of-mass energy is around the startup value for the ILC ( √ s 500 GeV) and where λ 5 adopts the (small) positive values allowed by the constraints. It should not come as a surprise that this behavior is similar to the SM result found in Fig. 5; indeed, being λ 5 small there can not be significant 2HDM enhancements (not even from tan β, which is 1) with respect to the corresponding SM process. Let us recall that for Set B the maximum value of λ 5 allowed by the vacuum stability bounds is λ 5 ∼ 0.65. Large negative λ 5 values give rise to significant (negative) quantum effects, and translate into loop-corrected cross sections σ (0+1) depleted 40% to 70% with respect to the tree-level predictions in the entire range from √ s = 500 GeV to √ s = 3 TeV. Since the tree-level h 0 Z 0 Z 0 coupling is equivalent to the SM coupling HZ 0 Z 0 for α = β − π/2, the right panel of Fig. 7 also illustrates the departure of the 2HDM loopcorrected cross section with respect to the tree-level SM cross-section, namely (σ SM . The behavior of σ(e + e − → h 0 Z 0 ) as a function of the Higgs boson mass M h 0 is presented in Fig. 8. This figure is the 2HDM counterpart of Fig. 5 corresponding to the SM case. We superimpose the tree-level and the loop-corrected cross-sections for the Set B of Higgs boson masses, by setting tan β = 1, α = β −π/2 and using three different values of λ 5 . The center-of-mass energy is settled at the fiducial value √ s = 500 GeV. Obviously, the raise of the Higgs boson mass M h 0 implies a reduction of the available phase space, so that the cross-section falls Finally, in Fig. 9 we display the tree-level contribution and the total (one-loop corrected) cross section (left panel), together with the relative radiative correction δ r (right panel), as a function of the parameter λ 5 . These plots have been generated in the same benchmark conditions as in Fig. 8 and using a fiducial ILC start-up center-of-mass energy √ s = 500 GeV. Let us remember that the red shaded area on the right hand side of that figure (specifically for λ 5 2.4) stands for the region excluded by the vacuum stability bounds, whereas the light gray shaded area on the left (λ 5 −10.5) signals the domain excluded by the perturbative unitarity bounds. The cross section is seen to follow the expected behavior, which we have derived from the estimate of the leading effect -cf. Eq. (16) -in combination with the structure of the trilinear coupling (23), namely a negative quadratic dependence of the scattering amplitude, which transates into σ(e + e − → h 0 Z 0 ) ∼ (a − bλ 2 5 ) 2 , emerging ultimately from the finite WF corrections to the external Higgs boson leg (notice the inclusion of the leading quartic corrections). For moderate negative λ 5 , the relative size of the one-loop quantum corrections can reach up to ∼ −50%. Positive δ r , however, can only take place for λ 5 close to zero.
A rather comprehensive survey of the predicted crosssections for different values of the tree-level coupling and different Higgs mass setups is presented in Table II. We display the results for the one-loop corrected cross-sections (in fb), together with the relative radiative corrections (cf. Eq. (21)), for both neutral Higgsstrahlung channels h 0 Z 0 and H 0 Z 0 . Once more, we set √ s = 500 GeV; and we work at tan β = 1 and maximum allowed |λ 5 |, since we are mostly interested in spotlighting the imprints of the 3H self-couplings in the quantum effects associated to the Higgs-strahlung mechanism. Let us also recall in passing that α = β − π/2 maximizes the tree-level h 0 Z 0 Z 0 coupling, whilst H 0 Z 0 Z 0 is optimal in the complementary regime, α = β. From the table we may read out that radiative corrections in these regimes are certainly large, regardless of the details of the chosen mass spectrum, the tree-level coupling and the actual channel under consideration. In a nutshell: the characteristic signature of the enhanced 3H self-couplings in the pattern of quantum effects on the Higgs-strahlung processes is rather universal and manifests in the form of a substantial depletion of the tree-level signal (typically in the range of δ r ∼ −20/ − 60 %). Fortunately, even under such a dramatic suppression the final loop-corrected cross sections stay at the level of a few tens of fb -thereby amounting to a non-negligible yield of ∼ 10 3 − 10 4 events per 500 fb −1 .
A note of caution should be given at this point in regard to some particular configurations of the H 0 Z 0 channel. It turns out that certain choices of Higgs boson masses (e.g. for MSSM-like mass splittings) may lie close, or simply beyond, the kinematical threshold for the decay process H 0 → h 0 h 0 ; that is to say, M H 0 ≃ 2 M h 0 . If we then consider a regime wherein α = β and sizable |λ 5 |, then the H 0 WF-correction undergoes a remarkable boost, which results from the combination of two independent sources of enhancement, namely i) the actual strength of the H 0 h 0 h 0 coupling -which is maximally enlarged in this regime; and ii) the kinematical enhancement due to the vicinity of the H 0 → h 0 h 0   Table I. threshold, which is reflected as a sharp peak in δ . Out of such corner in the parameter space, the higher-order effects involved in the above resummed form of the WFcorrections become harmless and totally inconspicuous and can be safely discarded. If, alternatively, we were considering a mass setup such that M H 0 > 2M h 0 , then the decay H 0 → h 0 h 0 would be open and, in the scenario of maximum H 0 h 0 h 0 coupling, it would furnish a very large width for H 0 . In this particular setup, the process would effectively boil down to the double Higgs-strahlung channels e + e − → H 0 Z 0 → h 0 h 0 Z 0 previously explored in the literature [22].
The following comment is in order to clarify the role played by the remaining (potential) sources of enhancement. We have seen that the various constraints restrict very significantly the range of allowed values of tan β, and enforce it to stay around 1 when |λ 5 | is maximum in the region λ 5 < 0. As a result, the contributions from the top quark and bottom quark Yukawa couplings cannot be augmented in the domain where the trilinear couplings are maximal. Remarkably enough, let us also emphasize that they cannot be enhanced even in the region where |λ 5 | is small or zero. To see why, notice that although in h 0 Z 0 H 0 Z 0 α = β − π/2 α = β − π/3 α = β − π/6 α = β − π/3 α = β − π/6 α = β  such region the parameter tan β can be very large or small and still be compatible with the constraints (cf. Fig. 6), then all the terms in the trilinear couplings (cf. Table  II of Ref. [32]) which are not proportional to λ 5 become of order one (i.e. cannot be promoted to high values) in the regime where the tree-level production cross-section for h 0 or H 0 is optimal. Not only so, in the very same regime the Yukawa couplings of both the top and bottom quarks cannot be enhanced either. One can easily check all these features explicitly by observing that, in the regions where the corresponding tree-level processes are maximal, the tan β-enhancements are canceled in all the couplings. We have checked numerically that the residual corrections (positive and negative) are merely of a few percent for any value of tan β. At the end of the day, we conclude that the only sizeable and eventually measurable quantum effects on the processes under study are those stemming potentially from the enhancement of the λ 5 coupling, not from the Yukawa couplings.
Let us finally emphasize that the complementarity between the neutral final states h 0 Z 0 /h 0 A 0 and H 0 Z 0 /H 0 A 0 , i.e. processes (1) and (2), is a unique chance for analyzing potentially big correlations between quantum effects in the 2HDM production cross sections, as there is no similar opportunity in the charged sector. Indeed, the final state charged counterparts H ± W ∓ are suppressed owing to the fact that the ZH ± W ∓ vertices are forbidden at the tree-level in any 2HDM extension of the SM [56] and therefore they can only be studied in more general extensions of the Higgs sector [57], or through loop-induced H ± W ∓ vertices in the 2HDM [58] and in the MSSM [59]. In these loop-induced mechanisms, the cross section for the associated H ± W ∓ production in a linear collider is rather meager -generally below 1 fb at the startup ILC energy (for charged Higgs masses comparable to the ones we have considered for the h 0 Z 0 /h 0 A 0 production) and within the region allowed by the current constraints. We therefore deem quite difficult to use the H ± W ∓ channel to extract additional information. In our opinion, the main task should be focused on performing preci-sion tests of the h 0 Z 0 /h 0 A 0 final states, together with the H 0 Z 0 /H 0 A 0 ones (if the mass of the heavy CP-even Higgs is not too large).

IV. DISCUSSION AND CONCLUSIONS
In this article, we have concentrated on the analysis of the production of neutral Higgs bosons in association with the Z 0 gauge boson at the future linac facilities within the framework of the general (nonsupersymmetric) Two-Higgs-Doublet-Model (2HDM). Our basic endeavor has been twofold: i) on the one hand to study the impact of radiative corrections to the final predicted rates; and ii) on the other hand to correlate such quantum effects with the enhancement potential of the 3H self-interactions, which are a genuine dynamical feature of the 2HDM -in the sense that it is unmatched to its supersymmetric counterpart. The upshot of our analysis singles out sizable (although negative) quantum effects which are correlated to the enhancement properties of the 3H self-interactions in the general 2HDM. Such large effects are identified fundamentally in the region of parameter space with tan β ≃ 1 and where the coupling |λ 5 | is at its maximum possible value compatible with the various constraints. The quantum effects reach typically δσ/σ ∼ −20%/ − 60%, for |λ 5 | ∼ 8/ ∼ 10 (λ 5 < 0), this being the crucial parameter that tunes the actual size of the 3H self-couplings. Let us stress that the most stringent limits on λ 5 are placed by the conditions of perturbative unitarity and vacuum stability 4 .
Although the vertex corrections to the h 0 Z 0 Z 0 /H 0 Z 0 Z 0 couplings are sensitive to such 3H self-interactions through Higgs-mediated one-loop diagrams, the most significant quantum effects are induced by the wave-function renormalization corrections to the h 0 /H 0 Higgs boson external lines. These (negative) effects triggered by the trilinear couplings are of order α ew λ 2 3H , and are therefore very responsive to changes in the value of the parameter λ 5 in the Higgs potential of the 2HDM. In turn, the gauge-boson and fermion one-loop contributions remain subleading (including the effects from the Yukawa couplings) as they show no remarkable departures from their SM analogues in the relevant regions of parameter space -which yield δσ/σ at the level of a few percent. Most important is also the fact that the combined analysis of such Higgs-strahlung events with the previously considered Higgs-pair production processes (h 0 A 0 , H 0 A 0 ) -see [32] for details -could be instrumental as a highly sensitive probe of the underlying architecture of the Higgs sector. At the fiducial center-of-mass energy value of √ s = 500 GeV, and when the genuine enhancement mechanisms of the 2HDM are active, the hA 0 events are remarkably strengthened at the one-loop level while the hZ 0 channels are simultaneously depleted.
We can assert that the described pattern of leading quantum effects emerges as a kind of "universal" feature of the 2HDM as far as the predictions for the Higgs-strahlung processes are concerned, meaning that this pattern is virtually independent of the details of the Higgs mass spectrum and of the actual production channel (h 0 Z 0 , H 0 Z 0 ). The same is true for the pairwise production channels [32]. Therefore they all can be physically significant provided the produced Higgs boson is not too heavy for the actual center-of-mass energy. Focusing on the Higgs-strahlung processes, the rise of large (and negative) radiative corrections to their crosssections within the 2HDM can be regarded as a characteristic imprint of a general two-Higgs-doublet model structure. Furthermore, the presence of a Z 0 -boson in the final state (and hence of its clear-cut leptonic signature Z 0 → l + l − ) should enable a rather comfortable tagging of these Higgs events following a method entirely similar to the original Bjorken mechanism [35] in the SM. Although the corresponding cross-sections at the higher operating energies planned for the future linear colliders are smaller (as compared to LEP), they are nevertheless sufficiently sizeable (typically in the range 10 − 40 fb) for a relatively comfortable practical measurement at the higher luminosities scheduled for these machines. For example, the associated decay of the accompanying neutral CP-even Higgs-boson would manifest basically in the form of either i) collimated and highly-energetic bquark or τ -lepton jets, for M h 0 2 M V ∼ 180 GeV; could provide an essential quantum handle enabling us to perform a proper identification of the kind of Higgs bosons produced in a linear collider.
In summary, we have analyzed the classical Higgsboson strahlung processes at linear colliders in the light of the general 2HDM, and we have elucidated very significant quantum imprints which spotlight, once more, the stupendous phenomenological possibilities that triple Higgs boson self-interactions could encapsulate in nonsupersymmetric extended Higgs sectors. To be sure, regardless of whether the LHC is finally capable to discover the Higgs boson(s), a paramount effort will still be mandatory in order to completely settle its experimental basis and to uncover the nature of the spinless constituents behind the electroweak symmetry-breaking mechanism. As we have shown, experiments at the fu-ture linear collider facilities could play a crucial role in this momentous endeavor.