Sensitivity to Triple Higgs Couplings via Di-Higgs Production in the 2HDM at e+e- Colliders

An important task at future colliders is the investigation of the Higgs-boson sector. Here the measurement of the triple Higgs coupling(s) plays a special role. Based on previous analyses, within the framework of Two Higgs Doublet Models (2HDM) type~I and~II, we define and analyze several two-dimensional benchmark planes, that are over large parts in agreement with all theoretical and experimental constraints. For these planes we evaluate di-Higgs production cross sections at future high-energy $e^+e^-$ colliders, such as ILC or CLIC. We consider two different channels for the neutral di-Higgs pairs $h_i h_j=hh,hH,HH,AA$: $e^+e^- \to h_i h_j Z$ and $e^+e^- \to h_i h_j \nu \bar \nu$. In both channels the various triple Higgs-boson couplings contribute substantially. We find regions with a strong enhancement of the production channel of two SM-like light Higgs bosons and/or with very large production cross sections involving one light and one heavy or two heavy 2HDM Higgs bosons, offering interesting prospects for the ILC or CLIC. The mechanisms leading to these enhanced production cross sections are analyzed in detail. We propose the use of cross section distributions with the invariant mass of the two final Higgs bosons where the contributions from intermediate resonant and non-resonant BSM Higgs bosons play a crucial role. We outline which process at which center-of-mass energy would be best suited to probe the corresponding triple Higgs-boson couplings.


Introduction
The discovery by the ATLAS and CMS collaborations of a new scalar particle with a mass of ∼ 125 GeV [1][2][3] is consistent with the existence of a Higgs boson in the Standard Model (SM). To date, all the measurements performed by the LHC of this Higgs boson are in agreement with the SM predictions. However, the uncertainties of the Higgs couplings measurements only reach a precision of roughly a ∼ 20%, hence there is room to beyond Standard-Model (BSM) Physics. There are plenty of BSM models that lead to scalar sectors with different features than the SM. In consequence, a main goal of the High Energy Physics community is to determine the nature of the Higgs mechanism and whether the discovered Higgs boson belongs to an extended BSM scalar sector. One of the main properties of the SM-like Higgs that remains yet undetermined is the value of its triple self-coupling, namely λ hhh , that it is only constrained to be inside the range −2.3 < λ hhh /λ SM < 10.3 at the 95% C.L. [4]. Many BSM models can induce important deviations in λ hhh with respect of the SM, therefore a more accurate measurement of the triple Higgs coupling will constitute a strong test of the SM. For recent reviews on the measurement of the Higgs couplings at future colliders see for instance [5,6]. In the case a BSM Higgs sector manifests itself, it will be a prime task to measure as well the BSM trilinear Higgs-boson couplings.
One of the most natural extensions of the SM Higgs sector is the "Two Higgs Doublet Model"(2HDM) (for reviews see, e.g., Refs. [7][8][9]) that consists in the addition of a second Higgs doublet to the SM one, with a ratio of the corresponding two vacuum expectation values given by tan β := v 2 /v 1 . This implies the existence of five physical Higgs bosons: two CP-even bosons h and H, usually with m h < m H , one CP-odd boson A and two charged Higgs bosons H ± . In Ref. [10] we performed an analysis of the possible size of triple Higgs couplings in the 2HDM being compatible with all the present theoretical and experimental constraints. For that analysis we assumed that the light CP-even Higgs-boson h is SM-like with a mass of m h ∼ 125 GeV. All other Higgs bosons were assumed to be heavier. To avoid flavor changing neutral currents at tree-level, a Z 2 symmetry is imposed [11], possibly softly broken by the parameter m 2 12 . Depending on how this symmetry is extended to the fermion sector, four types of the 2HDM can be realized: type I and II, lepton specific and flipped [8].
In Ref. [10] the 2HDM type I and II have been analyzed. We investigated the allowed ranges for all triple Higgs couplings involving at least one light, SM-like Higgs boson. The analysis was performed in several two-dimensional benchmark planes (i.e. all but two 2HDM parameters were fixed according to our definitions). We focused on the regions where all relevant theoretical and experimental constraints were fulfilled. For the SM-type triple Higgs coupling w.r.t. its SM value, κ λ = λ hhh /λ SM , we found allowed intervals of roughly [−0.5, 1.5] in the 2HDM type I and [0, 1] in type II. For the coupling of two light and one heavy CP-even Higgs bosons we found an approximate allowed interval of λ hhH ∈ [−1.4, 1.5] in the 2HDM type I and [−1.6, 1.8] in type II. Concerning the triple Higgs couplings involving two heavy 2HDM Higgs bosons, we found large allowed values for both 2HDM type I and II. For λ hHH , λ hAA and λ hH + H − we found maximum values of up to 15, 16 and 32, respectively. For further explorations several benchmark points were proposed as examples that exhibited large deviations from the SM-type triple Higgs coupling and/or large values of the triple Higgs couplings involving either two light and one heavy or one light and two heavy 2HDM Higgs bosons. Consequently, here our main motivation is to explore the potential of future e + e − colliders for the measurements of all these BSM triple Higgs couplings.
Future e + e − linear colliders, like the ILC [12] and CLIC [13] will play a key role for the measurement of the Higgs potential and to detect possible deviations from the SM with high precision [5,6,[14][15][16]. In particular, double Higgs boson production at future e + e − colliders appears to be the best way to explore possible deviations from the SM Higgs self-couplings (for recent reviews see, e.g., Refs. [16,17], and references therein). Within the 2HDM framework the largest effects of triple Higgs couplings on e + e − cross sections are also expected to be found in double Higgs production. These and other processes involving the 125 GeV Higgs boson at e + e − colliders have been explored in Refs. [18,19] in order to find signals from the 2HDM heavy Higgs bosons (for related work before the Higgs-boson discovery, see Refs. [20][21][22][23]).
In this paper we will analyze in detail the production of two neutral Higgs bosons in the 2HDM for the foreseen center-of-mass energies and luminosities of ILC and CLIC in the two main production channels: e + e − → Z * → Z + hh, hH, HH, AA, (1) e + e − → νν + hh, hH, HH, AA. (2) The first one is similar to the "Higgs-strahlung" channel of single Higgs production. The second one has an important contribution from the vector-boson fusion mediated subprocess, where the W W pairs (virtual, in general) are radiated from the initial e + e − together with the neutrinos: e + e − → W * W * νν. The processes in Eq. (2) also receive a contribution from the Z ( * ) mediated subprocesses, e + e − → Z ( * ) h i h j → ννh i h j , which are usually smaller than the contribution from W W fusion at the high energy colliders. Our main focus in this paper is addressed to study not just the cross section as a function of the 2HDM parameters, but also to explore the sensitivity to the various triple Higgs boson couplings via these two double Higgs production channels, at both ILC and CLIC. Our aim is to disentangle the role of these triple Higgs couplings in the mentioned double Higgs production processes and to investigate which process at which energy is best suited for an experimental determination. The calculation presented here of the total cross sections for these double Higgs production channels is performed in several two-dimensional benchmark planes (i.e. with all but two of the 2HDM parameters fixed), based on our results in Ref. [10]. The computations done here, as in Ref. [10], are performed with the help of the public codes MadGraph [24], FeynRules [25] and 2HDMC [26]. We demonstrate that it is possible to find regions with a strong enhancement of the hh production channels and/or with large production cross sections involving one light and one heavy or two heavy 2HDM Higgs bosons. First, we analyze in detail the mechanisms leading to these enhanced production cross sections. Second, we illustrate that sizable effects due to triple Higgs couplings can be seen in the cross section distributions as a function of the invariant mass of the di-Higgs final state (where we employed the code ROOT [27]). Finally, we discuss which process at which center-of-mass energy would be best suited to probe the corresponding triple Higgs-boson couplings.
Our paper is organized as follows. In Sect. 2 we briefly review the 2HDM, fix our notation and define the five benchmark planes used later for our investigation and summarize the constraints, which are (apart from updates) the same as in Ref. [10]. The di-Higgs production cross sections in the benchmark planes are presented in Sect. 3 and analyzed w.r.t. their dependence on the triple Higgs couplings. Finally, in Sect. 4 we analyze how the various triple Higgs couplings may be accessed experimentally through the invariant mass of the di-Higgs final state at the high-energy e + e − colliders. Our conclusions are given in Sect. 5.

The Model and the constraints
In this section we give a brief description of the 2HDM to fix our notation. We also review the theoretical and experimental constraints, which are (apart from updates) the same as in Ref. [10]. Finally we will define the benchmark planes for our analysis.

The 2HDM
We assume the CP conserving 2HDM (see Refs. [7][8][9] for reviews) whose potential can be written as: We denote the two SU (2) L doublets as Φ 1 and Φ 2 , where v 1 , v 2 are the real vevs acquired by the fields Φ 1 , Φ 2 , respectively, with tan β := v 2 /v 1 and they satisfy the relation v = (v 2 1 + v 2 2 ) where v 246 GeV is the SM vev. The eight degrees of freedom above, φ ± 1,2 , ρ 1,2 and η 1,2 , give rise to three Goldstone bosons, G ± and G 0 , and five massive physical scalar fields: two CP-even scalar fields, h and H, one CP-odd one, A, and one charged pair, H ± . Here the mixing angles α and β diagonalize the CP-even and -odd Higgs bosons, respectively.
To avoid the occurrence of tree-level flavor changing neutral currents (FCNC), a Z 2 symmetry is imposed, which is softly broken by the m 2 12 in the Lagrangian. The extension of the Z 2 symmetry to the Yukawa sector forbids tree-level FCNCs. This results in four variants of 2HDM, depending on the Z 2 parities of the fermions. We focus on type I and II, where the coupling to fermions are listed in Tab. 1 and Tab. 2. We will study the 2HDM in the physical basis, where the free parameters in Eq. (3) can be expressed in terms of the following set of parameters: c β−α , tan β , v , m h , m H , m A , m H ± , m 2 12 .
3 type I type II ξ u h s β−α + c β−α cot β s β−α + c β−α cot β ξ d,l From now on we use sometimes the short-hand notation s x = sin(x), c x = cos(x). In our analysis we will identify the lightest CP-even Higgs boson, h, with the one observed at m h ∼ 125 GeV. The couplings of the Higgs bosons to SM particles are modified w.r.t. the SM Higgscoupling predictions due to the mixing in the Higgs sector. In particular, the couplings of one neutral Higgs boson to fermions and to gauge bosons are given by: where m f , m W and m Z are the fermion mass, W mass and Z mass, respectively, and the factors in the couplings to fermions, ξ f h,H,A , and to gauge-bosons, ξ V h,H,A , are summarized in Tab. 2, for the 2HDM of type I and of type II 1 .
An important role in this paper is played by the couplings of the lightest CP-even Higgs boson with the other BSM bosons, concretely λ hhh , λ hhH , λ hHH and λ hAA . We define these λ hh i h j couplings such that the Feynman rules are given by: x L h j I v K 8 P m H g j N T K l w G f x 1 H X T g M W 9 5 a k u B C 5 C C W L c + b W E + D 7 6 n l a X W Q H i x f b I C B 4 r P U j G R C I a T d 6 H C L g 0 i e I 7 1 Q R F b B G o 8 p 3 3 G J M J P v t p X x S v l 2 p f Q 8 u 1 r M Y a Y g f 1 x Z o N 6 m 1 6 N B D D r F u D l 7 e 2 a H T 2 r Z L h w 3 c O 0 5 Q R H / u 2 c z b Q b T m r Y 7 g / 6 g K X q 5 8 d d N h 6 z r e N r + P Y k l K 3 M U h m W g d e A P C h N W o A x n m T W c l B o L Y H O Y Y W B b A T n q s G r + 6 5 I + S q S i J k X a z J v a C n K t F 3 l k N T m Y V O 9 y N f g / L i h N 8 j y s u C h K g 4 J Z i e W S M q N G 0 v p 5 0 J g r Z C Z b 2 A a Y 4 n Z L y l J Q w I x 9 R F t O T O Y F m J H f t 7 c O m 5 R 6 F 8 H Z m P z d U C 4 3 J 8 O + / 7 Q / f D 3 s H H X X g e 2 R B + Q h 6 R K f P C N H 5 B U 5 J m P C n I / O J + e L 8 9 X 9 7 H 5 3 f 7 g / V 1 L X W X 9 z n 2 y V + + s P E i / 5 + A = = < / l a t e x i t > = −i v n! λ hh i h j 1 Notice that we are using the same notation for the factors ξ f h,H,A as in our previous work of Ref. [10], and we are correcting here a typo in ξ u H and ξ d,l H of Table 2 of that reference.

4
where n is the number of identical particles in the vertex. The explicit expressions for the couplings λ hh i h j can be found in the Appendix of Ref. [10]. We adopt this convention in Eq. (7) so that the light Higgs triple coupling λ hhh has the same normalization as λ SM in the SM, i.e. −6ivλ SM with λ SM = m 2 h /2v 2 0.13. We furthermore define κ λ := λ hhh /λ SM . The couplings of the CP-even Higgs bosons strongly depend on c β−α . In particular, if c β−α = 0 one can recover all the interactions of the SM Higgs boson for the h state, what is known as the alignment limit. However, in the alignment limit in general one can still have BSM physics related to the extended Higgs sector, like hHH or ZHA interactions for example.

Experimental and theoretical constraints
In this section we will briefly summarize the various theoretical and experimental constraints considered in our scans.

• Constraints from electroweak precision data
Constraints from the electroweak precision observables (EWPO) can for "pure" Higgssector extensions of the SM, be expressed in terms of the oblique parameters S, T and U [28,29]. In the 2HDM the T parameter is most constraining and requires either m H ± ≈ m A or m H ± ≈ m H . In Ref. [10] we explored three scenarios: Here we will focus on scenario C. The 2HDM parameter space is explored with the code 2HDMC-1.8.0 [26].

• Theoretical constraints
Here the important constraints come from tree-level perturbartive unitarity and stability of the vacuum. They are ensured by an explicit test on the underlying Lagrangian parameters, see Ref. [10] for details. The parameter space allowed by these two constraints can be enlarged, in particular to higher BSM Higgs-boson mass values by the condition, • Constraints from direct searches at colliders The 95% confidence level exclusion limits of all important searches for BSM Higgs bosons are included in the public code HiggsBounds v.5.9 [30][31][32][33][34], including Run 2 data from the LHC. Given a set of theoretical predictions in a particular model, HiggsBounds determines which is the most sensitive channel and determines, based on this most sensitive channel, whether the point is allowed or not at the 95% CL. As input the code requires some specific predictions from the model, like branching ratios or Higgs couplings, that we computed with the help of 2HDMC [26].
• Constraints from the SM-like Higgs-boson properties Any model beyond the SM has to accommodate the SM-like Higgs boson, with mass and signal strengths as they were measured at the LHC. In our scans the compatibility of the CP-even scalar h with a mass of 125.09 GeV with the measurements of signal strengths at Tevatron and LHC is checked with the code HiggsSignals v.2.6 [35][36][37].
HiggsSignals provides a statistical χ 2 analysis of the SM-like Higgs-boson predictions of a certain model compared to the measurement of Higgs-boson signal rates and masses from Tevatron and LHC. Again, the predictions of the 2HDM have been obtained with 2HDMC [26]. Here, as in Ref. [10], we will require that for a parameter point of the 2HDM to be allowed, the corresponding χ 2 is within 2 σ (∆χ 2 = 6.18) from the SM fit: χ 2 SM = 84.73 with 107 observables. Many of the recent LHC Higgs rate measurements are now given in terms of "STXS observables". Contrary to our previous analysis in Ref. [10] the 2HDMC output can now allow the application of the (newly in HiggsSignals implemented) STXS observables. This results in substantially stronger limits on c β−α , particularly in the 2HDM type II. This will be visible in the fourth benchmark plane, see Sect. 3.1.
• Constraints from flavor physics Constraints from flavor physics have proven to be very significant in the 2HDM mainly because of the presence of the charged Higgs boson. Various flavor observables like rare B decays, B meson mixing parameters, BR(B → X s γ), LEP constraints on Z decay partial widths etc., which are sensitive to charged Higgs boson exchange, provide effective constraints on the available parameter space [38,39]. Here we take into account the decays B → X s γ and B s → µ + µ − , which are most constraining. This is done with the code SuperIso [40,41] where the model input is given by 2HDMC. We have modified the code as to include the Higgs-Penguin type corrections in B s → µ + µ − [42][43][44], which were not included in the original version of SuperIso. These corrections are indeed relevant for the present work since these Higgs-Penguin contributions are the ones containing the potential effects from triple Higgs couplings in B s → µ + µ − .
The effects of the various constraints on the 2HDM parameter space is discussed in detail in Ref. [10].

Benchmark planes
Based on the analysis in Ref. [10] we define five benchmark planes that exhibit an interesting phenomenology w.r.t. the di-Higgs production cross sections, Eqs. (1)  3 Cross section results

General considerations
In this section we present our analysis of the various di-Higgs production cross sections at e + e − colliders in the benchmark planes defined in Sect. 2.3. The analysis will be done for several collider options and official expected integrated luminosities, as summarized in Tab. 3. Concretely, we will take into account the ILC500 and the ILC1000 [12], as well as CLIC1500 and CLIC3000 [13]. The cross section predictions presented in this and the next section are calculated at tree-level precision with the help of the public code MadGraph5-2.7.2 [24] which generates and evaluates all contributing diagrams. All the information of the 2HDM regarding the model Lagrangian required by MadGraph was implemented with the Mathematica package FeynRules-2.3 [25]. In the computation, electrons and positrons were considered massless, so the diagrams where a Higgs boson couples directly to a fermionic line vanish. Our computation contains all the possible diagrams and does not rely at any point on the narrow width approximation. The width of all the Higgs bosons were calculated with the public code 2HDMC-1.8.0 [26]. We first start with a discussion on the cross sections as functions of the collider energy and next we continue with the results of the cross sections in the various benchmark planes introduced in the previous section for the different projected collider energies.  The complete set of diagrams that contribute to the two processes of our interest, e + e − → h i h j Z and e + e − → h i h j νν can be classified according to two different types of configurations: diagrams mediated by a virtual Z boson and diagrams mediated by virtual W W pairs. In e + e − → h i h j Z all contributing diagrams are of the first type, with the virtual Z attached to the initial e + e − pair, i.e., e + e − → Z * → h i h j Z. In e + e − → h i h j ν eνe , in Figure 1: Diagrams that include triple Higgs couplings in the di-Higgs production in the studied processes.  contrast, the two configurations contribute. There are diagrams mediated by a virtual Z attached to the initial e + e − , i.e. e + e − → Z * → h i h j ν eνe , and diagrams mediated by virtual W W pairs, i.e. e + e − → W * W * ν eνe → h i h j ν eνe . In this latter case, the subprocess involved is the so-called Vector Boson Fusion (VBF) given by W W → h i h j . The total number of diagrams in e + e − → h i h j Z are 7, 6, 7, and 7 for hh, hH, HH and AA, respectively. The total number of diagrams in e + e − → h i h j ν eνe and their separation in the two mentioned types (Z * mediated + VBF mediated) are 14 (7+7), 12 (6+6), 14 (7+7) and 12 (7+5) for hh, hH, HH and AA, respectively. While we do not show the diagrams here explicitly, they have all been included into our computations. It is important to note that only one diagram type among all these diagrams is the one carrying the triple Higgs couplings, concretely those depicted in Fig. 1. Thus, in order to get access to these couplings, a strategy has to be designed to disentangle the contributions from these particular diagrams among all the participating contributions in the total cross section. The discussion of the accessibility to the various triple Higgs-boson couplings will be presented in Sect. 4. In this subsection we focus first on the total cross sections.
The importance of the previously commented classification is well known in the literature, within the context of the SM, since their respective contributions to the total cross section behave very differently with the center-of-mass energy of the e + e − collisions, √ s. The Z * mediated configurations decrease with energy, whereas the W * W * mediated ones increase with energy. This is indeed the reason why for e + e − colliders with energies at and above the TeV scale the diagrams with VBF configuration may dominate the h i h j νν production rates. The VBF dominance is well known in the case of di-Higgs production e + e − → hhνν in both the SM (h = H SM ), with κ λ = 1, as well as in BSM models (where h represents the BSM Higgs boson at ∼ 125 GeV) with κ λ = 1 (see, for instance, Refs. [16,17] and references therein). Consequently, this is expected to also happen in the present case within the 2HDM. One relevant difference with respect to the SM case is that in the 2HDM the extra Higgs bosons (other than h) also participate in these processes: H and H ± , participate in the W * W * mediated subprocess, and H and A enter in the Z * mediated subprocess. One illustrative example of the different behavior of the cross sections with the collider energy, for the processes analyzed here, is presented in Fig. 2. Here we show e + e − → Zh i h j (left plot) and e + e − → ννh i h j (right plot) for a particular point of the 2HDM (type I) parameters given by m H = m A = m H ± = 500 GeV, tan β = 10, c β−α = 0.2, and m 2 12 = 24000 GeV 2 . Other examples will be presented in Sect. 4. (See Ref. [45] for an additional discussion.) The plots in Fig. 2 illustrate the commented different behavior with energy of the cross sections above about √ s ∼ 1 TeV. In particular, we see that hhνν grows with energy, whereas hhZ decreases with energy (compare the two light blue lines in Fig. 2). Indeed they have a similar behavior with √ s than the corresponding SM cross sections (red lines in Fig. 2). The 2HDM hh production rates are enhanced with respect to the SM ones by a factor of about 3, for this particular choice of 2HDM input parameters. This enhancement is due to the additional contributions from the extended 2HDM Higgs sector, including the effects from the new diagrams with intermediate heavy Higgs bosons (H and A in hhZ and H, A and H ± in hhνν). The size of the enhancement depends strongly on the choice of the 2HDM parameters, most prominently on c β−α away from the alignment limit. As will be discussed below, this enhancement mainly originates from the diagrams where the intermediate heavy Higgs bosons, H and A, can be produced on-shell (i.e. in the case of hhZ for , but values of κ λ = 1 can also play a role. Consequently, the largest contributions to the enhancement w.r.t. the relevant invariant mass distribution occur in the resonant region, i.e. for m hh ∼ m H or m hZ ∼ m A , respectively. In contrast, the contributions from an intermediate H ± in hhνν production, which appear only in t-channel like diagrams, do not exhibit a resonant behavior in this process. The case of hH shows a different pattern with √ s for higher center-of-mass energies above ∼ 1200 GeV in this particular point. Both cross sections for hHZ and hHνν decrease with energy (green lines in these plots). The cross section for hHZ is above the hhZ one by a factor of about 10, and the decrease with energy is well understood in terms of the process being mediated by a virtual Z * propagating in the s-channel. The cross section for hHνν is below the hhνν one, and it shows a qualitatively different behavior with √ s at high center-of-mass energies. The hHνν cross section decreases with √ s and follows a similar pattern as in the hHZ channel, in clear contrast with the hhνν cross section that, as discussed above, increases with √ s. This suggests that at large energies, hHνν is dominated by the Z * mediated diagrams whereas hhνν is dominated by the W * W * mediated diagrams. Indeed, we have checked this dominance explicitly by comparing the cross sections h i h j νν with electron neutrinos (which include both types of mediated diagrams), versus those with muon neutrinos (which include only Z * mediated diagrams). We have found that at high energy, above around 1 TeV, σ(hhνν) σ(hhν eνe ), and σ(hHνν) 3σ(hHν µνµ ), confirming our reasoning above. Besides, in order to understand better the VBF mediated contributions we have explored the cross section of the corresponding subprocess of the two cases, hh and hH. For the case of hhνν, we have explicitly checked that the cross section of the subprocess W W → hh tends to a constant value at high energies (as in the SM case) due to a strong cancellation of the contributions from the t, u and contact subprocess diagrams, which are separately growing with energy. The contribution from the charged Higgs is very suppressed in this case of large m H ± . For the case of hHνν, however, there is no contact diagram in W W → hH and the strong cancellations occur instead between the W and the charged Higgs boson contributions to the t and u channels separately. This leads to a remnant contribution that decreases with energy, as they do the subprocess diagrams with the h and the H propagators in the s-channel. The charged Higgs contribution is important at high energies in this case.
The production cross sections of HH and AA in both the Z and the neutrino channels are numerically nearly equal due to our choice m H = m A , that implies λ hHH ∼ λ hAA and λ HHH ∼ λ HAA . There are some diagrams present only in the AA production, involving the interaction vertex hAZ, whereas some others appear only in HH production with ZZH and W W H vertices (the latter only present in the neutrino channel). However, all the diagrams mentioned above carry a c 2 β−α factor and consequently, they are suppressed with respect to the rest of the diagrams, that are equivalent for HH and AA productions if m H = m A . Furthermore, the cases of HHνν and AAνν production are clearly dominated by the VBF subprocess and their increase with energy are the most pronounced ones (purple and yellow lines in Fig. 2, which appear indeed superimposed). However, the rates are lower than for the hhνν and hHνν cases due to the large values of m H,A for this particular point (set in this plot to m H,A = 500 GeV. To gain access to them will require the highest energy option, i.e. 3 TeV as will be discussed in Sect. 4. Due to the strong similarity of HH and AA production, below we will often discuss only ZHH and HHνν.
Looking at all the production modes, the important contribution of VBF in the h i h j νν channels play an important role in the present study. The high energy collider options, particularly for 3 TeV, may offer an efficient window to improve the sensitivity to the triple Higgs couplings, which are mediated via the diagram on the right in Fig. 1, belonging to the VBF diagrams subset. This will be analyzed in detail in Sect. 4. In the following, we present our results for the total cross section. We will first discuss the three 2HDM type I benchmark planes, followed by a discussion of the two 2HDM type II planes.

2HDM type I
We present the results obtained in the three benchmark planes for the various production cross sections as given in Eqs. (1) and (2) and various energies and integrated luminosities according to Tab. 3. Some planes are omitted when the results do not show any relevant variation. In each plane we indicate with (the interior of) a solid black line the part of the parameter space that is allowed taking all theoretical and experimental constraints as given in Sect. 2.2. Small differences w.r.t. Ref. [10] are due to updated experimental constraints.
The predictions for the cross sections depend on the various triple Higgs couplings. Therefore in Fig. 3 we reproduce the results from Ref. [10] for κ λ = λ hhh /λ SM (left), λ hhH (2nd column), λ hHH λ hAA (3rd column) and we add new predictions for λ HHH λ HAA (right) in the 2HDM type I benchmarks planes 1 (upper row), 2 (middle row), 3 (lower row). The triple Higgs couplings involving h are defined according to Eq. (7). The corresponding definitions are set for the other Higgs triple couplings not involving h, which we do not include here explicitely for shortness. Shown in black solid lines are the exterior boundaries of the allowed regions of the 2HDM parameter space (see above). The red line in the left column indicates κ λ = 1.

hh production
We start with the production of two light (SM-like) Higgs bosons. In Fig. 4 we show the production cross sections at √ s = 500 GeV for e + e − → Zhh (left) and e + e − → hhνν (right) in the benchmark planes 1-3 (top to bottom). The SM cross sections for the two processes are 0.158 fb and 0.034 fb, corresponding to 630 and 135 events, respectively, for the luminosity given in Tab. 3. At this low energy of √ s = 500 GeV, the e + e − → hhνν channel is dominated (as in the SM case) by the Z * mediated configurations, i.e. by e + e − → Z * → hhνν, and the VBF configurations being subdominant do not play a relevant role. Therefore, the potential accessibility to the two involved triple Higgs couplings λ hhh and λ hhH will mostly come from the diagrams to the left in Fig. 1.
Three main effects can change the 2HDM prediction w.r.t. the SM predictions. The first one are deviations in λ hhh from its SM value, i.e. κ λ = 1. This will modify the contributions of type e + e − → Z * → Zh * → Zhh where the intermediate h * is off-shell. The second one are additional diagrams involving in particular an intermediate H, like e + e − → Z * → ZH → Zhh, which can be produced on-shell in contrast to the previous case. Correspondingly, large effects are expected for large values of λ hhH , especially when the H can be produced on-shell. The last contribution comes from the sub-processes e + e − → Z * → Ah → Zhh (→ ννhh), that can also give sizable contributions to the cross section if the virtual A can propagate on-shell, despite the fact that they are ∝ g ZhA ∝ c β−α . However, these A mediated diagrams do not carry any sensitivity to the triple Higgs couplings.
In the alignment limit, reached for c β−α = 0, the 2HDM cross sections reproduce exactly the SM cross sections. Here one has κ λ = 1 and λ hhH = 0, as can also be seen in the first two columns of Fig. 3. Also, as mentioned above, g ZhA ∝ c β−α and therefore the A * mediated diagrams also vanish in the alignment limit. However, as can be seen in the color code of Fig. 4, SM-like cross sections are also reached away from c β−α = 0 in all three planes. In these regions λ hhH ∼ 0 (see the 2nd column of Fig. 3) due to cancellations, as has been discussed in Ref. [10]. Furthermore, λ hhh does not strongly deviate from the SM value, as can     be seen in the left column of Fig. 3, again due to cancellations in the various contributions to λ hhh , see the discussion in Ref. [10]. On the other hand, large enhancements of the cross section predictions are found for large values of |λ hhH |. Furthermore, in particular for a near resonant di-Higgs production, m H ∼ 2m h ≈ 250 GeV, the cross sections are found to be strongly enhanced. Correspondingly, the largest cross sections inside the allowed regions are found in the benchmark plane 3 for m H ∼ 250 GeV. Here the enhancement goes up to ∼ 6.5 (13) times the SM cross section for hhZ (hhνν). Here it is important to note that in these regions of the parameter space it is not sufficient to calculate the cross section in the narrow width approximation, (NWA) i.e. via HZ or Hνν production with the subsequent decay H → hh. As will be discussed in detail in Sect. 4, the reason for the failure of the NWA to provide an accurate prediction of the total cross section for the di-Higgs bosons production in e + e − collisions can be understood as follows. Outside the resonant region in the invariant mass of the final Higgs pairs the remaining non resonant diagrams are very relevant and account for a sizeable contribution to the total cross section that cannot be neglected.
Also cross sections substantially smaller than in the SM can be found. This can happen in particular for hhZ in the benchmark plane 1 in the "tip" of the allowed region where ∼ 0.3 times the SM cross section is found. In the planes 2 and 3 the largest suppression are found for c β−α ∼ 0.15, where ∼ 0.7 and ∼ 0.5 times the SM prediction is found. These regions correspond to the smallest allowed values of λ hhh , see the left column of Fig. 3.
We now turn to the results for √ s = 1000 GeV, which are shown in Fig. 5 with the same color coding and line styles as in Fig. 4. At this center of mass energy, both channels have comparable cross sections in the SM. For hhZ and hhνν production one finds SM cross sections of 0.120 fb and 0.097 fb, corresponding to 962 and 779 events, for the luminosity given in Tab. 3. One difference respect to the previous case is that at this higher energy of √ s = 1000 GeV, the VBF subprocess enters in hhνν more relevantly than for √ s = 500 GeV.
Then the sensitivity to λ hhh and λ hhH enters via the two types of diagrams (left and right) in As for √ s = 500 GeV the largest enhancements within the allowed regions are found for m H ∼ 250 GeV and larger values of c β−α . In the benchmark plane 3 we find enhancements of ∼ 6 and ∼ 13 w.r.t. the SM predictions for hhZ and hhνν, respectively. Within benchmark plane 2 enhancements of ∼ 2 are found for c β−α ∼ 0.25 (0.15) for the two production cross sections, respectively. In benchmark plane 1 with m H = 1000 GeV the impact of diagrams involving the H are negligible, and the cross sections are controlled by κ λ . The different signs of interference with the diagrams not involving the hhh vertex lead to a different numerical behavior of the two cross sections. In the "tip" of the allowed regions, where κ λ reaches the smallest possible values (in the allowed region, see the upper left plot of Fig. 3), we find a suppression of ∼ 0.6 w.r.t. the SM cross section for hhZ, while for hhνν we find an enhancement of ∼ 6.   In Fig. 6 we show the results for √ s = 1500 GeV. At this center of mass energy, the hhνν channel, due to its t-like-channel nature, is larger than the s-like-channel dominated hhZ production. The cross section for hhνν is clearly dominated by VBF configurations, as in the SM case. For hhZ and hhνν production one finds SM cross sections of 0.077 fb and 0.239 fb, corresponding to 193 and 597 events, for the luminosity given in Tab. 3. However, despite the substantially larger cross section for hhνν because of the smaller anticipated integrated luminosity at √ s = 1500 GeV (at CLIC) w.r.t. the high luminosity expected at √ s = 1000 GeV (at ILC) the absolute number of hhνν events remains smaller. The overall pattern of the cross sections relative to the SM cross sections, as given by the color code in Fig. 6 is qualitatively similar to the one for √ s = 500, 1000 GeV, but the deviations from the SM tend to be larger in the 1500 GeV case than in the previous ones of 500 GeV and 1000 GeV. In the c β−α -tan β plane shown in the first row within the allowed parameter space we find enhancements up to 1.3 and 3 w.r.t. the SM predictions for the hhZ and hhνν channel, respectively. These largest values are realized, as before, in the "tip" of the allowed region, where κ λ deviates most from unity. While the interference with the SM-type contributions is similar to the previous center-of-mass energies now also the additional contributions from an intermediate H play a role, leading, e.g., to an enhancement of hhZ, contrary to the case of √ s = 1000 GeV. In the c β−α -m 2 12 plane shown in the middle row relative enhancements of ∼ 5 can be found. As will be discussed in the next section, the contribution of an on-shell H with H → hh plays an important role here, and thus does the size of λ hhH . We find that the impact of the latter coupling extends up to m H < ∼ 1000 GeV, as can be seen in the c β−α -m (m = m H = m A = m H ± ) plane in the third row of Fig. 6. The maximum effect, as for smaller center-of-mass energies is found for the resonant di-Higgs production at m H ∼ 2m h ≈ 250 GeV. Here we find maximum enhancements of ∼ 6 and ∼ 16 for the hhZ and hhνν channel, respectively. The corresponding values of λ hhH do not vary strongly in this part of the parameter space and are found to be O(2), see the 2nd column of Fig. 3. For m H > ∼ 1000 GeV, where λ hhH becomes less relevant, despite reaching its largest possible values (see the lower plot in the 2nd column of Fig. 3) we even find a decrease of the hhZ channel up to ∼ 0.8, whereas the hhνν channel remains enhanced in the allowed parameter space.
We finish the analysis of hh production with √ s = 3000 GeV, the highest energy potentially reachable at CLIC, in Fig. 7. At this center-of-mass energy hhνν clearly dominates over hhZ, with SM production cross sections of 0.033 fb and 0.819 fb, corresponding to 164 and 4097 events, for the luminosity given in Tab. 3. The hhνν rates are now fully dominated by the VBF configurations where the new window to test λ hhh and λ hhH is clearly open. As before, the overall pattern of the cross sections relative to the SM cross sections, as given by the color code in Fig. 7 is qualitatively similar to the one for smaller energies, but now with deviations from the SM tending to be larger for 3000 GeV than for the 500, 1000 and 1500 GeV cases. In the c β−α -tan β plane we find enhancements of up to ∼ 2 and ∼ 2.5 for hhZ and hhνν production, respectively. Slightly larger enhancements are found in the c β−α -m 2 12 plane. As for √ s = 1500 GeV, close to the resonant region, the on-shell decay H → hh contributes substantially to the total cross section and this will be important to reach sensitivity to λ hhH , as will be shown in Sect  we are close to the alignment limit, where λ hhH is very tiny. In the lowest row of Fig. 7 we find enhancements of up to ∼ 4.5 and ∼ 10, respectively.

hH production
In this subsection we analyze di-Higgs production of one light and one heavy CP-even Higgs boson. As expected, we find that in the alignment limit the production cross sections become exactly (analytically) zero. Regarding the values of the two triple Higgs couplings involved, λ hhH and λ hHH , the first one vanishes in the alignment limit, as was discussed above, but the second one does not vanish in this limit. The vanishing of the contribution from the H * mediated diagram to this hH cross section, in the alignment limit, occurs in this case because of the vanishing of the H coupling to the gauge bosons in this limit. Consequently, a key point here is to investigate the accessibility to these λ hhH and λ hHH away from the alignment limit.
In largest cross sections that we find in the allowed parameter region is about ∼ 8 fb and ∼ 2 fb for hHZ and hHνν production, respectively. As before, λ hhH and λ hHH do not play an important role here.
In Fig. 10 we present the results for √ s = 1500 GeV. This center-of-mass energy is sufficiently high such that all benchmark planes can in principle exhibit relevant production cross sections. This large cross sections for both hHZ and hHνν processes at 1500 GeV can also be seen in Fig. 2. However, as discussed above, these large cross sections found occur mainly due to the intermediate on-shell A production that decays to Zh. Consequently, we find that for the parameters in our benchmark plane 1 with m H = m A = 1000 GeV no large cross section can be observed within the allowed parameter space. Only for very large values of c β−α and tan β (outside the allowed region) cross sections reaching the fb region are reached. This is due to the extreme values that λ hHH λ hAA can take.
In the c β−α -m 2 12 plane as shown in the middle row of Fig. 10 we have m A +m H = 1200 GeV and thus relevant cross sections are indeed found, reaching up to ∼ 0.5 fb for hHZ and ∼ 0.1 fb for hHνν production. The results are found to be effectively independent of m 2 12 , but depend only on c β−α . This reflects the parametric dependence of the production cross  section on the couplings g ZHA g ZhA ∝ s β−α c β−α . The results in the benchmark plane 3, presented in the lower row of Fig. 10 show the same pattern as for the previously analyzed center-of-mass energies. As expected, relevant production cross sections are now found up to m A = m H ∼ 750 GeV. The largest cross sections are found around m H = 300 GeV, reaching ∼ 4 fb and ∼ 0.8 fb for hHZ and hHνν production, respectively. The small but non-zero cross section found above m A = m H ∼ 750 GeV for |c β−α | 0.2 (i.e. outside the allowed region) originate in regions of the parameter space, where λ hHH is very large, as can be seen in the lower left plot of Fig. 3, which can give a relevant contribution to the total cross section. As it can be seen in Fig. 3, this coupling can reach values up to ∼ 15 in this part of the parameter space.
We finish the hH analysis with Fig. 11, where we show the results for √ s = 3000 GeV.
This high center-of-mass energy allows for on-shell intermediate A production in all three benchmark planes and consequently leads to relevant cross sections in all plots. In the c β−α -tan β planes shown in the first row we find, as for √ s = 1500 GeV the largest cross sections for large c β−α and tan β (or very small tan β in the case of hHνν), but outside the allowed regions. Within these regions the largest cross sections found are ∼ 0.5 fb and ∼ 0.2 fb for hHZ and hHνν production, respectively. As before, this is correlated with the extreme values that λ hHH λ hAA take in this part of the parameter space.
The cross sections shown in the c β−α -m 2 12 planes in the middle row of Fig. 11 exhibit the same pattern as for √ s = 1500 GeV, i.e. dominated by the contributions originating in e + e − → HA and thus no relevant dependence on m 2 12 . The largest cross sections are found to reach ∼ 0.9 fb and ∼ 0.2 fb for hHZ and hHνν, respectively.
In the third benchmark plane we find relevant cross sections up to m H = m A ∼ 1500 GeV, but for smaller masses, < ∼ 700 GeV, now smaller cross sections are observed w.r.t. √ s = 1500 GeV. In this part of the parameter space again the contributions originating in e + e − → HA dominate, which are suppressed with 1/s. The largest cross section values are found again for m H = m A ∼ 300 GeV and reach ∼ 1.25 fb and ∼ 0.25 fb. However, the overall pattern of the predicted cross sections, in particular for hHνν, differs from the ones for √ s = 1500 GeV, which can be understood as follows. While the hHZ channel is again dominated by the e + e − → HA → H(hZ) resonant subprocess, in the hHνν channel it can be seen that there is a similar (or even larger) prediction than in the hHZ channel for m H = m A > ∼ 900 GeV (albeit outside the allowed area). The main difference between the hHZ and the hHνν channels at high energies of O(3 TeV) is the different influence of the triple Higgs couplings. The hHZ channel is dominated, as discussed above, by the intermediate HA state and the subsequent decay A → hZ, which are not sensitive to triple Higgs couplings. In contrast, the hHνν not only receives relevant contributions from the diagrams shown in the left of Fig. 1 with the subsequent decay Z → νν, but as well from VBF subprocesses such as W W → hH, where the triple Higgs couplings do enter via the s-channel diagrams with either virtual h or virtual H propagating. This type of topology, mediated by W W fusion, grows with energy and thus contributes relevantly to the e + e − → hHνν cross section. Consequently, at high energies one expects to gain access to the triple couplings involved, λ hhH and λ hHH . Since for larger values of m H one finds larger values for λ hHH than for λ hhH , as can be seen in the lower row of Fig. 3, one expects a higher sensitivity to λ hHH in this part of the parameter space. The final reach to these couplings is the combined result  of the large coupling effect vs. the suppression of the corresponding virtual Higgs (either h or H) propagating in the s-channel of this VBF subprocess. For the parameter region explored in Fig. 11 λ hHH can reach large values up to ∼ 10 within the allowed area, as it can be seen in the lower 3rd column plot of Fig. 3. The virtual H boson, propagating close to its mass shell, can thus contribute strongly to the hHνν production cross section. Conversely, diagrams with a virtual h are suppressed by relatively smaller triple Higgs couplings of O(1), as can be seen in the lower 2nd column plot of Fig. 3, as well by the propagation of a very off-shell h boson. Overall, this yields a higher sensitivity to λ hHH than to λ hhH in the allowed parameter region, as will be discussed in more detail in Sect. 4.

HH (AA) production
In this subsection we analyze di-Higgs production of two heavy Higgs bosons. This can be either HHZ/HHνν or AAZ/AAνν. Since for our choice of m A = m H the di-CP-odd production is always very similar to the di-CP-even production we only describe the latter. In these channels, the alignment limit does not necessarily imply a zero cross section, and in general no resonant enhancement of cross sections can take place in our benchmark planes, since the intermediate Higgs bosons are always off-shell.
The results for √ s = 500 GeV are shown in Fig. 12. For this low center-of-mass energy only the low mass region of benchmark plane 3 yields kinematically allowed production cross sections. For HHZ production we find very small cross sections going up to ∼ 0.05 fb in the allowed region for m H ∼ 160 GeV. The HHνν channel does not exhibit any non-negligible cross sections. In the relevant parameter space λ hHH and λ HHH are of similar size, but do not exceed the h boson mediated one. Indeed we find that the contribution involving an intermediate A boson can be somewhat larger that the h mediated one.
In Fig. 13 the results for √ s = 1000 GeV are shown. As for √ s = 500 GeV, only the low mass region of benchmark plane 3 yields kinematically allowed production cross sections. As before, diagrams involving λ hHH give a larger contribution than the ones involving λ HHH . Higher cross sections are found for smaller m H : HHZ production reaches ∼ 0.

2HDM type II
We now turn to the two benchmark planes defined for the 2HDM type II, see Sect. 2.3 (planes 4 and 5). Other planes were found to not exhibit any relevant variation. Furthermore, most of the plots also in these two planes showing the results for the various production cross sections as given in Eqs. (1) and (2) and various energies and integrated luminosities according to Tab. 3, are also omitted, since the results do not shows any relevant variation either. Within the allowed parameter space, taking into account all theoretical and experimental constraints as given in Sect. 2.2, we find numerical results for the cross sections of e + e − → hhZ and e + e − → hhνν, at all the studied energies given in Tab. 3, are in general very close to the respective SM predictions, with a deviation of at most ∼ 5%. On the other hand, the hH production channels yield cross sections below 0.01 fb with no relevant variation in the allowed parameter space, for all center-of-mass energies. The same applies to the heavy di-Higgs production (where again the results for HH are very similar to the results for AA production)     for √ s = 500, 1000, 1500 GeV. Only for e + e − → HHZ and e + e − → HHνν at the highest energy, √ s = 3000 GeV, some relevant production cross sections exceeding ∼ 0.1 fb within the allowed parameter space are obtained. Thus, we choose this √ s = 3000 GeV center-of-mass energy in our discussion below. As we have said previously, in the two channels for HH production, HHZ and HHνν, the two involved triple couplings are λ hHH and λ HHH , and they enter via the diagrams with an intermediate off-shell h and H, respectively. Thus, the size of these two couplings will be the most relevant issues for the HH production rates. We show the results for these two triple Higgs couplings (with λ hHH λ hAA and λ HHH λ HAA ) in Fig. 16, for the two chosen planes, 4 and 5. As in the previous subsection in each plane we indicate with a solid black line the part of the parameter space (inside the black line) that is allowed taking all theoretical     in both benchmark planes 4 and 5. The regions where the cross section becomes larger correspond to the regions where λ hHH reaches its largest "allowed" values around ∼ 6 − 8. This indicates that the total cross section receives a relevant contribution from the diagrams containing λ hHH , in particular from the VBF-like contributions, which dominate at this large center-of-mass energy. On the other hand, over large parts of the allowed parameter space the cross section does not exceed 0.05 fb, and no access to triple Higgs couplings can be expected. Because of these small cross sections for the 2HDM type II, even in the best case of √ s = 3000 GeV, we will study only one "best-case" scenario in the following section.

Sensitivity to triple Higgs couplings 4.1 Benchmark points
After having explored the total cross sections at the future e + e − colliders for the several channels of two Higgs bosons production in the 2HDM of type I and II, we now turn to the potential sensitivity to the involved triple Higgs couplings. As was discussed above, each channel involves different triple Higgs couplings. The processes with two light Higgs bosons hh in the final state involve λ hhh (with κ λ = λ hhh /λ SM ) and λ hhH . The processes with hH involve λ hhH and λ hHH . Those with HH involve λ hHH and λ HHH , and those with AA involve λ hAA and λ HAA . Since the predictions for HH are very similar to those of AA we will focus here just on cases of hh, hH and HH. As we anticipated before, the total cross section of these processes is not sufficient to infer the effects of the triple Higgs couplings. In this work we propose to access to these couplings through another observable: the differential cross section with respect to the invariant mass of the final Higgs pair h i h j , which we will study in the following sections.  In order to study the sensitivity to the various λ h i h j h k we have chosen specific benchmark points (BPs) within the 2HDM, which are in agreement with all present data. These  experimental constraints such that sizeable values of the triple Higgs couplings are realized. The five selected BPs and their input parameters are specified in Tab. 4. The first four, BP1, BP2, BP3 and BP4 are for type I and the fifth one, BP5, is for type II. We have also included in this table the corresponding output values of the relevant triple Higgs couplings and the relevant heavy Higgs boson widths. It should be noted that in all these points the largest triple Higgs coupling is λ hHH , so that one might expect the largest sensitivity to this coupling. However, this is not the case, as will be discussed in this section, since the other input parameters, especially the mass of the heavy Higgs bosons and the value of c β−α , also play a very relevant role and will enter significantly in the conclusion to the sensitivity to the triple Higgs couplings.
To complete the description of these five BPs we show the predictions of the total cross section for the various channels as a function of the e + e − collider energy in Fig. 18, for BP1 through BP4, and in Fig. 19 for BP5. As a general remark, one can see that a similar pattern of the cross section behavior with energy is found here for these five selected BPs as for the point shown in Fig. 2, see the discussion in Sect. 3.1. Generically, at high energies the rates for h i h j Z decrease with energy and for h i h j νν increase with energy. With the exception of hHνν that first decreases with energy (as seen in our figures) up to around 3 TeV, and then reaches a nearly flat (but slightly increasing) behavior at very high energies above 3 TeV (not shown in our plots). We have checked this flat behavior with energy also at the WBF subprocess level, i.e., the cross section σ(W W → hH) tends to a constant value at very high √ s. This flat behavior also happens in the σ(W W → hh) case, similarly to the SM case, but it is reached at lower energies in the hh channel than in the hH one.
We find the following hierarchies in the size of the cross sections: (i) the hhZ channel provides the largest σ(h i h j Z) in the low energy region from threshold up to ∼ 600 GeV for BP1, up to ∼ 1000 GeV for BP2, up to ∼ 1200 GeV for BP3, up to ∼ 2000 GeV for BP4, and in the full energy range for BP5; (ii) hHZ provides the largest σ(h i h j Z) at energies above those mentioned values; (iii) hhνν provides the largest σ(h i h j νν) for all the studied BPs and at all energies. However, as we will discuss below, just looking for the largest cross sections is not sufficient to reach the best sensitivity to the triple Higgs couplings.
Our strategy proposed here to explore the accessibility to the triple Higgs-boson couplings is outlined as follows. It is obvious that for this analysis the relevant quantity is not that of the total cross section for production of h i h j pairs, but another quantity where the effect from the triple Higgs coupling is enhanced. There are several of those quantities, but we choose here what we believe is the most promising one: the cross section distribution with respect to the invariant mass of the produced Higgs boson pair. The main idea with this distribution is to look for specific windows of this invariant mass where the 2HDM rates will be the most sensitive ones to the triple Higgs coupling involved. Obviously, the extreme case with maximum sensitivity will happen when the invariant mass of the produced h i h j pair, m h i h j is close to the intermediate Higgs boson mass, m h k , since in that case the diagram carrying the λ h k h i h j coupling dominates over the others due to the resonant behavior of this diagram (see Fig. 1). Here it should be noted that we are not assuming the on-shell production of this intermediate Higgs boson h k . As we have checked explicitly, it is not a good approximation to estimate the total rates for h i h j X production by using the naive NWA, i.e., by calculating the h k X production rates times the BR(h k → h i h j ). The reason for this bad approximation is the fact that outside the window of the m h i h j invariant mass that is close to m h k the rest of diagrams (not carrying λ h k h i h j ) contribute very significantly, as will be shown in the following. In the next subsections we will discuss separately the results for the cases of hh and hH, HH.

Sensitivity to triple Higgs couplings in hh production
In this section we analyze the sensitivity to triple Higgs couplings in the hh production channels. We present the results of the cross section distributions as a function of the invariant mass of the hh pair, m hh , for all the collider energies considered in this work. First we comment on the distributions for the points selected in 2HDM-type I and latter we present the results for the point selected in 2HDM-type II.
The set of plots in Fig. 20, Fig. 21, Fig. 22, Fig. 23, are the results in the 2HDM-type I case, for √ s = 500 GeV, 1000 GeV, 1500 GeV, and 3000 GeV, respectively. In the vertical axis on the right of these plots, the corresponding hh event rates using the luminosities of Tab. 3 are also shown. The plots in the left of all these figures correspond to the hhZ channel, and the plots on the right correspond to the hhνν channel. It should be noted that for each energy we have only included the predictions for the BPs where the corresponding heavy Higgs boson mass m H is clearly reachable. Concretely, BP1 with m H = 300 GeV is reachable at the four energies; BP2 with m H = 500 GeV and BP3 with m H = 600 GeV are reachable at 1000 GeV, 1500 GeV and 3000 GeV; and BP4 with m H = 1000 GeV is only reachable at 1500 GeV and 3000 GeV. For completeness, we also include the corresponding predictions of the 2HDM total cross section in all these cases, both in the interior boxes of these figures and in Tab. 5 and Tab. 6 for the hhZ and hhνν cases, respectively. The predictions for the SM case are included as well, for comparison. For the forthcoming discussion, it is also important to analyze separately the contributions from the various diagrams in these m hh distributions. Thus, in all these figures we display separately: the complete dσ/dm hh 2HDM rates (in red), the contributions from the diagrams containing the λ hhh (in light blue), those from the diagrams containing the λ hhH (in dark blue), the sum of these two latter containing the two triple Higgs couplings (in purple), the sum of all the rest of diagrams not containing the triple couplings (in yellow), and the complete dσ/dm hh SM rates (in black). Notice that comparing the purple lines with the dark blue and light blue lines one could also determine the effect of the interference between the h and H mediated diagrams. One can also see from this comparison that the mentioned interference effect is negligible in the mass invariant m hh region close to the resonant peaks where the H mediated diagrams clearly dominate. We start with the discussion regarding the sensitivity to λ hhh in all these distributions. First, we find that the maximum sensitivity to this triple coupling appears in the low m hh region, slightly above the 250 GeV threshold of hh production. This feature also happens in the SM case, and is well known in the literature, see, e.g., Ref. [16]. The contribution from λ hhh modifies the profile of the distribution in that region close to threshold, and the distortion is larger for larger values of κ λ . In fact the contribution from the diagram containing this λ hhh produces an important interference effect which is constructive in the hhZ channel, whereas it is destructive in the hhνν case. This is seen in our plots by comparing the red lines (total) with the yellow lines (without triple Higgs couplings) and with the light blue lines (with just λ hhh ) in that region close to the hh threshold. The red lines are above the yellow lines in hhZ, but they are below the yellow lines in hhνν. The size of this interference effect from λ hhh is larger for larger triple coupling values (or equivalently larger κ λ values) and is clearly correlated with the contribution from the diagram with an intermediate virtual h (light blue line). This interference effect from the triple Higgs coupling also happens in the prediction of the cross section distribution for the case of SM di-Higgs production, as it is well known in the literature (see, for instance, [16,17] and references therein). Indeed, it has been studied as a very efficient strategy at CLIC to determine the value of κ λ within the context of the SM, with a high precision of around ±0.1 [16]. In the present 2HDM case, the figures show that the effect from λ hhh is larger than in the SM in the hhνν channel, at the colliders with the largest energies, i.e. at √ s = 3 TeV. The biggest contribution from λ hhh and the largest size of the interference effect can be seen in Fig. 23 in comparison to Fig. 20, Fig. 21 and Fig. 22. Comparing the various BPs we see that the maximum effect is produced in the points with largest λ hhh (hence, largest κ λ ), as expected. In particular, we have checked that the point BP3 with κ λ = 1 provides similar rates as the SM in this close to threshold m hh region (compare red and black lines). Overall, we find that the effect from λ hhh within the 2HDM is similar to the corresponding one in the SM case, and therefore the sensitivity to this triple Higgs coupling is expected to be comparable to that concluded in the literature from the studies of the sensitivity to κ λ in di-Higgs SM production at e + e − colliders. The only exception can occur when m H is not much larger than 2 m h and the contribution from an intermediate H boson gives an important contribution in the relevant part of the m hh distribution (as discussed below). We next comment on the sensitivity to λ hhH that appears, as previously mentioned, via the diagrams with an intermediate H. The first clear signal from this intermediate heavy Higgs boson H contribution in all these figures is seen as a resonant peak emerging in the invariant mass region with m hh close to m H . The dark blue line accounts for this H mediated contribution which displays the resonant behavior and provides the dominant contribution to the total result (red line) in the narrow region around the resonance . The second observation, as already mentioned, is that using the NWA to compute the total cross section by σ(hhZ) = σ(HZ) × BR(H → hh) and σ(hhνν) = σ(Hνν) × BR(H → hh) does not provide an accurate prediction. We have checked this failure of the NWA prediction, using the corresponding cross section for single H production and the BR(H → hh) for each BP. The main reason for this failure is that the remaining diagrams, other than the one mediated by the intermediate heavy CP-even Higgs boson, contribute very significantly to the total cross section. These remaining diagrams together explain the non-resonant part of these distributions, i.e. they explain the red lines outside the resonant region (usually called 'the continuum'). In particular, we have checked explicitly that one of the most relevant contributions outside the resonant peak is provided by the diagram where the CP-odd A boson is the intermediate boson. We have checked that the 'big step shape' that is best seen to the right of the resonant peak in some of the distributions is due to the A mediated contribution, which indeed dominates the total production cross section. This can be seen, for instance, in the distributions in the hhZ channel for BP1 at 1000 GeV, 1500 GeV and 3000 GeV, that display most clearly these 'big steps'.
The beginning and ending of these 'steps' due to the A boson mediated diagrams depend on the momentum and energy of the two Higgs bosons in the final state coming from this type of diagrams, which can be written in short as e + e − → hA ( * ) → hhZ. In that case, we can differentiate one of the final light Higgs boson, labeled as h 1 , that comes from the process e + e − → h 1 A and the other Higgs boson, labeled as h 2 , that comes from the decay A → h 2 Z. On one side, the momenta of h 1 and A are completely determined by √ s, m A and m h and they are given by On the other side, the momentum of h 2 depends on the energy and momentum of its mother particle A, namely p A and E A , the mass of the final states m h and m Z and on the angle between the h 2 and the A, θ h 2 A : where p CM is the final momentum of h 2 and Z in the center-of-mass frame of the A. If m A p CM /m h p A < 1, the maximum and the minimum of p h 2 is reached for θ h 2 A = 0 for the positive and negative sign of the Eq. (10) respectively. Otherwise, if m A p CM /m h p A > 1 the only physical solution is the positive sign of Eq. (10) and the maximum is reached for θ h 2 A = 0, whereas the minimum is reached for θ h 2 A = π. Therefore, these extrema in p h 2 define the range in the invariant mass of the two final light Higgs bosons where the effects from the A boson manifest, and are given by where E h i is the energy of h i , with i = 1, 2, and the angle θ h 1 h 2 must satisfy θ h 1 h 2 + θ h 2 A = π. For instance, for BP1 at 1000 GeV we find the beginning and ending of the step at m min,max hh = 456, 907 GeV, respectively, in agreement with what is seen in Fig. 21.
On the other hand, this relevant A mediated diagram for the hhZ channel provides the largest contribution when the intermediate A is produced on-shell and then it decays to hZ. Indeed, we have found the corresponding emergent A resonant peak which can be seen clearly in the distributions with the alternative invariant mass variable, m hZ (but not in the distributions with m hh ). Here we restrict ourselves to one example for the m hZ invariant mass distribution displayed in Fig. 24 for BP1 at √ s = 1000 GeV. However, we do not investigate further the appearance of these peaks here because they are not sensitive to the triple Higgs couplings which are the main objective in our analysis. These A resonant peaks obviously could provide interesting information on the BSM physics induced by the A bosons within the 2HDM, other than the triple Higgs couplings.
Overall, the results of the various contributions to the differential cross section given by the color lines in our plots in Fig. 20, Fig. 21, Fig. 22 and Fig. 23 demonstrate explicitly that the dominant contribution in the region of m hh close to the m H resonance stems clearly from the diagram carrying the λ hhH (dark blue line), and therefore the height of these resonant peaks in comparison with the basis lying on the bottom of the resonance (given by the yellow line) provides the most efficient access to this particular λ hhH . We furthermore remark that these figures display some sensitivity to the sign of λ hhH , and not only to its absolute value. There is an asymmetry in the peak of the H boson since the non-resonant diagrams are not negligible and the relative sign of the resonant diagram changes when m hh = m H (caused by a sign flip in the H-boson propagator). This effect can be seen even more clearly in Fig. 25 (where we show the "family of BP3" for various values of c β−α , which will be introduced later at the end of this subsection): in the points where λ hhH > 0 (corresponding to the points with c β−α = 0.1, 0.12, 0.14) the cross section at the left of the H peak is visibly larger than at the right of the peak. On the other hand, in the cases with λ hhH < 0 (corresponding to the points with c β−α = 0.16, 0.18, 0.2) the cross section at the left of the peak is larger than at the right. The opposite effect due to the sign of λ hhH happens in the "H peak" of the hhZ channel, as can be seen for example in the left panels of Figs. 20, 21, 22 and 23.
To provide a more quantitative estimate of the sensitivity to the leading triple coupling λ hhH in this case of hh production, and in absence of a more realistic study including real backgrounds from detector effects etc., we propose here a more theoretical approach to define 'the significance of the signal'. 'The signal' here refers to the specific events under the resonant peaks which are the only ones really carrying the sensitivity to λ hhH . Concretely, we set this interval as the one in between the crossings of the yellow and dark blue lines.            The events considered in this study are those after the Higgs decays into bb pairs, therefore producing typically a final state signature with either four b-jets and a Z boson (presumably easily detectable at the clean e + e − collider environments) in the case of hhZ; or with four b-jets and missing transverse energy, E T , (from the undetected neutrinos) in the case of hhνν. The number of events under the resonant peaks, N R 4bZ for hhZ, and N R 4b & E T for hhνν, are then extracted from these plots for the four BPs and the various considered energies. We also extract from these plots the event rates from the non resonant contributions (the so-called continuum contributions), in the corresponding m hh resonant region, named N C 4bZ and N C 4b & E T , for hhZ and hhνν respectively, which are defined by the yellow lines in these plots. Finally, we also extract the predictions of the corresponding SM event rates, N SM 4bZ and N SM 4b & E T , respectively, for comparison. All these event rates are obtained for the corresponding luminosities in Tab. 3 and include the BR factors of the two Higgs decays, i.e. (0.58) 2 .
More realistically, this m hh invariant mass should rather be considered as the invariant mass of the 4-b-jets, m 4b out of which each bb pair reconstructs each of the two Higgs bosons. This could be done, presumably, by the proper cuts in the m 2b invariant mass window close to m h . But as stated above, we simplify our analysis and choose to work with the theoretical variable m hh . However, since the more realistic experimental detection of these signals will require the tagging of the b-jets and the tagging of the missing transverse energy in the hhνν case, or the Z boson in the hhZ case, we have re-evaluated all the above rates, taking into account the detection acceptance, A, and the b-jet tagging efficiency, b . Concretely, we evaluate the acceptance as follows: where N without cuts and N with cuts are the total event rates predictions without and with cuts applied, respectively. For this study we again use MadGraph and apply the following cuts to the b-jets (meaning cuts applied to the b-quarks since we are not dealing with jets), missing transverse energy (for the neutrino channel) and Z transverse momentum (for the Z channel), which are similar to those given in [14,46]: where, p b T , η b are the transverse momentum and pseudo rapidity of each of the four b-jets, p Z T the transverse momentum of the Z, and ∆R bb ≡ (∆η bb ) 2 + (∆φ bb ) 2 , with ∆η bb and ∆φ bb , are the separations in pseudorapidity and azimuthal angle of the two b-jets, respectively, which are identified coming from one light Higgs decay. With these cuts we get the values of the acceptances for the various BPs and for the four studied energies that are displayed in Tab. 5 for the hhZ channel and in Tab. 6 for the hhνν channel. As one can see, the acceptance varies with energy, as expected, and we find the best values of around 0.70 for the lowest energies. At the highest energy of 3000 GeV the acceptance for hhZ gets reduced to around 0.1 and for hhνν to around 0.5. These acceptances are close to the SM ones (slightly better, indeed) which we have also included in these tables, for comparison. For the numerical estimates of the more realistic event rates that are involved in our study of the resonant peaks for the various BPs,N , we then apply the reduction factors, given by A and b , to the previous event rates, N , as follows: where we set the value of the b-jet tagging efficiency in our numerical evaluations to b = 0.8 (see, e.g., Refs. [14,47]). To evaluate the size of the λ hhH effect we compute the following ratio R: which is our theoretical estimator of the 'sensitivity' to λ hhH . We collect all the results for the correspondingN R ,N C and R in Tab. 5 and Tab. 6 for hhZ and hhνν, respectively. As one can see in these tables, the number of signal events are quite significant in most cases. For the hhZ channel,N R 4bZ decreases with the collider energy and leads the largest R value for the BP1 at 1000 GeV. For BP2 and BP3 the rates are clearly lower and the largest R values found are again at 1000 GeV. For BP4 the rates are too low to be detected and we find no sensitivity to λ hhH in this hhZ channel. For the hhνν channel the signal rates are clearly larger (except for BP1 at 500 GeV) and, in contrast to hhZ, increase with the energy collider. We find high sensitivity to λ hhH in all the studied cases, except for BP4 at 1500 GeV, where the rates are too low. The highest sensitivities, corresponding to the highest R values above 100, are found for BP1 for the three center-of-mass energies 1000, 1500 and 3000 GeV. BP3 at 3000 GeV also shows a large R value of 100. BP2 and BP4 reach their maximum R values, 48 and 34 respectively, also at 3000 GeV.
All in all, we conclude that the process e + e − → hhνν seems very promising to give access to the triple λ hhH coupling at all the studied energies, since very high values of our estimator R are found in all the studied 2HDM points (except BP4 at 1500 GeV). The highest sensitivities, indicated by the highest values of R in this hhνν channel, are found at 3000 GeV. The channel hhZ can also give access to this triple λ hhH coupling, but at the  Table 7: R 4b & E T for BP3 for different c β−α at √ s = 3000 GeV. Also shown are λ hhH , the total widths Γ H , total cross sections, σ, and event numbers,N , as defined in the text. "-" indicates values that cannot be evaluated due to a too small number of events.  Table 8: GeV. Also shown are λ hhH , the total widths Γ H , total cross sections, σ, and event numbers,N , as defined in the text. "-" indicates values that cannot be evaluated due to a too small number of events.
lower energy colliders and for the 2HDM points with a relatively light H boson. The highest sensitivities in this case are reached for BP1 at 500 GeV and 1000 GeV.
After having explored the sensitivity to λ hhH in all the four BPs of the 2HDM type I, we will now analyze the relevance of the choice of c β−α for the found sensitivities. We have then repeated our previous study by means of the R variable, but now varying c β−α for each BP, whereas the other input parameters are kept at their original value. In this part of the analysis we focus now on the channel and energy with the highest sensitivity, namely hhνν at 3000 GeV, and present two examples of BP related points: the case of 'the family of BP3' in Tab. 7, and the case of 'the family of BP4' in Tab. 8. We have also included in these tables the corresponding predictions of the triple coupling λ hhH for those 'family points'. The corresponding predictions for the cross section distributions with the invariant mass m hh are displayed in Fig. 25 and Fig. 26 for 'the family of BP3' and the family of BP4', respectively. In these figures, we see clearly how both the resonant peak and 'the continuum' evolve with c β−α in each of the two studied cases. To conclude on the significances of these new resonant signals and their sensitivity to λ hhH , we again compare the rates under the resonant peak (red lines) with the rates under the continuum (yellow lines), in the m hh interval defined by the two crossings of the yellow and dark blue lines. Then we compute the corresponding R 4b & E T value using the same acceptances and b-tagging efficiencies as before. The corresponding results are shown in Tab. 7 and Tab. 8 for 'the family of BP3' and 'the family of BP4', respectively. One can see that the sensitivity, given by R 4b & E T , directly correlates with the value of |λ hhH |. Since the benchmark points were defined to reach large values of |λ hhH | (within each family), we find that the highest sensitivity found is precisely for the parent BP3 and BP4. Similar results and conclusions are found for the other 'family points' (not shown for brevity).

Sensitivity to triple Higgs couplings in hH and HH production
In this section we study the sensitivity to the triple Higgs couplings in hH and HH production. We have selected in these two cases the highest collider energy of 3 TeV, where the effects of the triple couplings have been found to be the largest ones. We start with the hH case. In this channel the two couplings involved are λ hhH and λ hHH , and their effects appear via their contributions in the diagrams that are mediated by either an intermediate h or H, respectively. However, these intermediate Higgs bosons are always off-shell. Therefore, they do not produce resonant peaks in the relevant invariant mass distribution, which in this case refers to the invariant mass of the final hH pair, m hH . The results of the distributions dσ/dm hH and the corresponding event rates for both channels hHZ and hHνν are shown in Fig. 27 in the left and right column, respectively. The color code in Fig. 27 is as follows: red lines are the complete dσ/dm hH taking into account all diagrams, dark blue lines are the contributions from the diagrams mediated byh, green lines are those from the diagrams mediated byH, purple lines are the contributions from the sum of the two latter mediated by h and H, and yellow lines are the contributions from the rest of diagrams other than the h and H mediated ones. The first conclusion from these plots is that the hHZ channel does not provide sensitivity to any of the two involved triple Higgs couplings. This channel is dominated by the diagram e + e − → Z * → H A → H hZ, as we commented in Sect. 3.2.2. In fact, the cross section distributions shows a "step" shape similar to the one analyzed in the case of hh production and the beginning and ending of this contribution can be reproduced if m H is introduced accordingly in Eq. (11). The second conclusion from Fig. 27 is that the hHνν channel can provide sensitivity to the triple couplings. The relevant effect from these two couplings, λ hhH (dark blue lines) and of λ hHH (green lines) appear in the close to threshold region, i.e. with m hH slightly above m h + m H . The two contributions together provide the largest contributions in this region, leading to the purple lines in these plots, and interfere negatively with the remaining contributions (yellow lines), producing a decrease in the total rates (red lines lying below the yellow lines). Comparing the contributions of the two triple couplings separately, we see that in all the BPs the largest contribution is from the green line, indicating that the largest sensitivity found in this channel is from to the triple coupling λ hHH . This feature is in clear correlation with the fact that for these chosen points this particular λ hHH coupling gets large values, varying from 2 in BP1 and BP2 to 6 in BP3 and BP4. The relative interference of the h and H contribution depends on the relative sign of λ hhH and λ hHH . A positive (negative) interference is found for λ hhH · λ hHH positive (negative).
Next we comment on the sensitivity to the triple couplings in HH production. In this     Fig. 28, in the left and right column, respectively. The color code in these plots is as follows: red lines are the complete dσ/dm HH taking into account all diagrams, green lines are the contributions from the diagrams mediated by h, pink lines are those from the diagrams mediated by H, purple lines are the contributions from the sum of the two mediated by h and H, and yellow lines are the contributions from the rest of diagrams other than the h and H mediated ones. The most relevant effects from these two couplings, λ hHH (green lines) and of λ HHH (pink lines) appear as before in the region close the threshold of HH production, which in this case means m HH slightly above 2m H . The first important conclusion from these plots is that both channels, HHZ and HHνν, notice the effects of these triple Higgs couplings. Comparing the two involved couplings, we see that the largest effect is clearly from λ hHH which is the largest among all the studied triple Higgs couplings in all the four BPs. Its values range from 2 to 6, as can be seen in Tab. 4. Accordingly, the green lines in Fig. 28 provide a much larger contribution than the pink lines in the two channels, HHZ and HHνν for all the BPs. Comparing the event rates at 3 TeV (see the right vertical axis), both the total rates and the ones in the mentioned close to threshold region, one can see that they are sizeable (except for BP4) in HHνν. On the other hand, they remain relatively small for HHZ production. Consequently, the sensitivity to λ hHH at this large energy appears mainly in the channel with neutrinos whereas there is no sensitivity to λ HHH . Comparing the red, the yellow and the purple lines in these HHνν plots, we see that the purple ones dominate over the yellow lines, i.e. the contributions from the rest of diagrams other than the ones mediated by h or H, except for the largest values of m HH . Since the purple lines are slightly above the red lines in HHνν, we conclude that the interference from the triple Higgs coupling is destructive in this channel, but with a small effect due to the large numerical difference between the sets of contributing diagrams. In contrast, in the HHZ channel, the purple line is below the red line, therefore, the interference is constructive. As in the previous hH case, we do not discuss further the sensitivity to the triple Higgs couplings.
To estimate the effects of the triple Higgs couplings on the cross section distributions for the hH and HH channels more quantitatively, as we have done for λ hhH in the hh case, is interesting but involved. This would require a more refined analysis in the close to threshold m hH and m HH regions respectively, including the study of the more realistic final states considering the Higgs bosons decays and the corresponding backgrounds, which is beyond the scope of this work. Thus, we do not go further in these hH, and HH channels.

Sensitivity to triple Higgs couplings in the 2HDM type II
In this section we focus in the 2HDM type II and present the results of the relevant cross section distributions for h i h j production in the corresponding invariant mass m h i h j for the benchmark point BP5 (see input values in Tab. 4). As before we focus on the highest envisaged e + e − collider energy of 3 TeV. The predictions for all the channels are summarized In our analysis of BP5 in the 2HDM type II we start with the h i h j Z channels. We see some (potential) sensitivity to λ hhh in the close to hh threshold region for hhZ production, see the light blue line that clearly dominates in this region. However, the event rates are too low and therefore unable to provide sensitivity to this coupling. The case of hHZ is clearly not sensitive to the triple Higgs couplings at this high energy. In the case of HHZ production there is an appreciable effect from the green line, involving λ hHH that gives the dominant contribution in the close to threshold region. However the predicted rates are yet too small and do not provide sensitivity enough to this coupling. We now turn to the h i h j νν channels shown in the right column of Fig. 29, where, as expected, overall larger production cross sections are obtained. We find relevant effects in the hhνν channel from λ hhh , as can be seen in the light blue line in the region close to the m hh = 2m h threshold. These contributions interfere negatively with the diagrams not involving triple Higgs couplings (yellow lines), yielding a sizable number of events. However, the prediction for the complete differential cross section (red line) is very close to the SM one (black line), reflecting the fact that κ λ 1 in BP5. Consequently, we do not expect deviations in BP5 from the sensitivity expected to κ λ within the SM. In the case of hHνν, although there is an apparent effect from λ hHH , and to a lesser extent from λ hhH , they do not seem statistically significant, because the rates are also very low in this case. Finally, the channel HHνν seems more promising at this high energy. The green line, involving λ hHH clearly dominates over the rest and provide practically the complete cross section (red line). Here one should note that for BP5 λ hHH is the largest one among all triple Higgs couplings with a value of ∼ 7. The event rates are yet significant, indicating that the effect from this large λ hHH could be measured in this channel.

Conclusions
An important task at future colliders is the investigation of the Higgs-boson sector. Here the measurement of the triple Higgs coupling(s) plays a special role. Based on previous analyses [10], within the framework of Two Higgs Doublet Models (2HDM) type I and II we have analyzed several benchmark planes that are over large parts in agreement with all theoretical and experimental constraints. These comprise electroweak precision data, tree-level unitarity, stability of the vacuum, bounds from direct searches for BSM Higgs bosons, the rate measurements of the 125 GeV Higgs boson at the LHC, as well as flavor constraints. While three benchmark planes are based on the 2HDM type I, for type II we defined two benchmark planes. The relevant parameters are summarized in Sect. 2.3. As one important characteristic, we have set all heavy Higgs-boson masses equal, m H = m A = m H ± to simplify the analysis. For these planes we investigated the di-Higgs production cross sections w.r.t. triple Higgs couplings at future high-energy e + e − colliders, such as ILC or CLIC. We considered two different channels e + e − → h i h j Z and e + e − → h i h j νν with h i h j = hh, hH, HH, AA. The first one is similar to the "Higgs-strahlung" channel of single Higgs production. The second one has an important contribution from the vector-boson fusion mediated subprocess, W + W − → h i h j . It also receives a contribution from the Z mediated subprocess, e + e − → Zh i h j , with Z → νν, which is usually smaller than the contribution from W W fusion at the high energy colliders. The triple Higgs couplings relevant for our analysis are λ hhh (with κ λ = λ hhh /λ SM ), λ hhH , λ hHH and λ HHH . Due to the fact that HH production is very similar to AA production we investigated only the former. Conclusions for λ hAA (λ HAA ) are very similar as for λ hHH (λ hAA ).
In the first part of our study we have evaluated the h i h j Z and h i h j νν total production cross sections for √ s = 500, 1000, 1500 and 3000 GeV, i.e. the (relevant) energy stages foreseen for the ILC and CLIC. We have discussed the various production cross sections as a function of the 2HDM parameters in the selected benchmark planes.
As expected from the known SM cross sections, our study of the 2HDM cross sections as a function of the center-of-mass energy shows that, in general, the h i h j Z cross sections decrease with increasing center-of-mass energy, while h i h j νν cross sections increase, despite different types of diagrams can contribute for i = j. We furthermore related the pattern of enhanced or suppressed cross sections to the sizes of the relevant triple Higgs couplings. As anticipated, the 2HDM hh production cross section approaches the SM value in the alignment limit, c β−α → 0. On the other hand, large enhancements of hh production w.r.t. the SM cross sections are found at some 2HDM parameter configurations. Our analysis indicates that this enhancement is due to the additional contributions from the extended 2HDM Higgs sector, including the effects from the new diagrams with intermediate heavy Higgs bosons. The size of the enhancement depends importantly on the choice of the 2HDM parameters, most prominently on c β−α away from the alignment limit. These intermediate heavy Higgs bosons are H and A for the hhZ channel, and H, H + and A for the hhνν channel. We have found that the enhancement mainly originates from the diagrams where the intermediate heavy Higgs bosons, H and A, can be produced on-shell (i.e. in the case of hhZ for √ s > ∼ m H + m Z or > ∼ m A + m h , respectively), but values of κ λ = 1 can also play a role. Consequently, the largest contributions to the enhancement w.r.t. the relevant invariant mass distribution occur in the resonant region, i.e. for m hh ∼ m H or m hZ ∼ m A , respectively. The effect of the H ± boson, on the contrary, does not produce such a resonant behavior in the e + e − → hhνν process, and it does not involve any triple Higgs coupling. The effect of the on-shell A boson is resonant in the relevant mass invariant distribution in both hhZ and hhνν channels but it is not sensitive to the triple Higgs couplings either. Out of all these intermediate heavy Higgs boson contributions, the most relevant one in both hhZ and hhνν channels is the one mediated by an on-shell H, where the relevant mass invariant distribution shows the H resonance at m hh ∼ m H and displays the unique sensitivity to the triple Higgs coupling λ hhH .
Concerning hH production, these channels do also receive relevant contributions from diagrams with an intermediate A boson: e + e − → Z * → H A ( * ) → H hZ and e + e − → Z * → H A ( * ) → H hZ ( * ) → H hνν, and we find sizeable cross sections at all the studied collider settings. But these A mediated contributions are not sensitive either to the triple Higgs 52 couplings. Again we find sensitivity only via the H mediated contributions that in this case give access to probe λ hHH . Regarding HH production, the cross sections at the planned colliders are found to be small, < ∼ 1 fb, in the allowed parameters space, and therefore the sensitivity to the involved couplings, λ hHH and λ HHH , is correspondingly smaller.
In the second part of our analysis we investigated in detail the mechanisms leading to the enhancements of the production cross sections identified in the first part, and explored how to reach the best sensitivity to the involved triple Higgs couplings via the study of the cross sections distributions with respect to the invariant mass of the corresponding final Higgs bosons pair. We have concentrated on five benchmark points, which are contained in our previously defined benchmark planes (with the exception of BP2, which is inspired by benchmark plane 2, but has a slightly smaller value of m H = m A = m H ± and a correspondingly slightly adjusted value of m 2 12 ). Four are defined in the 2HDM type I (BP1, BP2, BP3 and BP4) while one is defined in type II (BP5). For all production channels we found, depending on the center-of-mass energy, that the contributions from an intermediate H Higgs boson play a crucial role, offering interesting access to (BSM) triple Higgs couplings involving this heavy Higgs boson. Specially, when this H is produced on-shell and appears as a resonant peak in the corresponding invariant mass distributions of the final Higgs bosons pair. Furthermore, we outlined which process at which center-of-mass energy would be best suited to probe the corresponding triple Higgs-boson couplings. In general the 2HDM type I allows for larger BSM triple Higgs couplings and thus better prospects for their measurements than the 2HDM type II.
Focusing on the 2HDM type I, hh production approximately offers access to λ hhh similar to the corresponding SM determination. On the other hand, in the region of m hh ∼ m H , a large (resonant) enhancement of the cross section is found, giving clear access to λ hhH . In order to quantify the sensitivity to these triple Higgs couplings we chose to analyze the event rates for the final states containing four b-jets, i.e. we assumed the final h bosons decays into bb pairs. We then evaluated the number of 'signal' events (corresponding to the resonant diagrams involving λ hhH ) and compared with the 'continuum' events (corresponding to the other diagrams, which do not resonate), convoluted the relevant acceptances and efficiencies. Only benchmark points with low m H have a significant dependence on λ hhH in the hhZ production channel. On the other hand, the hhνν channel offers more significant dependence on λ hhH , for all considered values of m H , in particular at the higher center-of-mass energies. Therefore, the hhνν channel shows significant sensitivity to λ hhH , at CLIC. These findings were further substantiated by varying c β−α for two of the benchmark points (within the region allowed by all constraints). The significance for λ hhH in the hhνν channel scaled directly with |λ hhH |.
Regarding hH and HH cross sections, they were evaluated for the highest foreseen energy at CLIC, √ s = 3000 GeV, as it is expected that this center-of-mass energy will yield the best sensitivities to the BSM triple Higgs couplings. Some access to λ hhH and λ hHH may be possible in the region close to threshold of m hH > ∼ m h + m H in the hHνν production channel. Here λ hHH was found to be more relevant than λ hhH . For HH production also the HHZ channel shows a relevant dependence on λ hHH , but the number of expected events for the anticipated CLIC luminosity remains too small for a measurement. HHνν, on the other hand, may give additional access to λ hHH , in particular in the region close to threshold of m HH > ∼ 2 m H .
Overall we find that higher energy e + e − colliders can offer access to λ hhh in the 2HDM similar to the SM case. Depending on the BSM Higgs-boson masses that lead to resonant BSM Higgs contributions in the mass invariant distributions of di-Higgs production also a significant dependence on BSM triple Higgs couplings is found. The best sensitivity found among all the triple Higgs couplings it to λ hhH via hhνν production and reaches the highest values at the highest energy colliders. Considering all triple couplings, ILC and CLIC may both shed light on the Higgs potential of models with extended Higgs-boson sectors.