Exploring sizable triple Higgs couplings in the 2HDM

An important task at future colliders is the measurement of the triple Higgs coupling. Depending on its size relative to the Standard Model (SM) value, certain collider options result in a higher experimental accuracy. Within the framework of Two Higgs Doublet Models (2HDM) type I and II we investigate the allowed ranges for all triple Higgs couplings involving at least one light, SM-like Higgs boson. We take into account theoretical constraints (unitarity, stability), experimental constraints from direct Higgs-boson searches, measurements of the SM-like Higgs-boson properties, flavor observables and electroweak precision data. We find that the SM-type triple Higgs coupling w.r.t. its SM value, $\lambda_{hhh}/\lambda_{\rm SM}$, can range between $\sim -0.5$ and $\sim 1.5$. Depending on which value is realized, the HL-LHC can compete with, or is clearly inferior to the ILC. We find the coupling $\lambda_{hhH}$ between $\sim -1.5$ and $\sim 1.5$. Triple Higgs couplings involving two heavy Higgs bosons, $\lambda_{hHH}$, $\lambda_{hAA}$ and $\lambda_{hH^+H^-}$ can reach values up to ${\cal O}(10)$, roughly independent of the 2HDM type. This can lead to potentially strongly enhanced production of two Higgs-bosons at the HL-LHC or high-energy $e^+e^-$ colliders.


Introduction
One of the key points of our study when exploring the parameter space of the 2HDM type I and II under the given constraints is the following: the primary focus of our explorations was to find allowed parameters that lead to either large non-SM triple Higgs boson couplings, or to large deviations from unity in the ratio of the light triple Higgs-boson coupling w.r.t. its SM value, λ hhh /λ SM . In particular, we have explored scenarios with relatively heavy masses m H , m A and m H ± near 1 TeV, but not enforcing the so-called alignment limit, cos(β − α) → 0 (see, e.g., [51]). Furthermore, we have investigated the dependences of the allowed triple Higgs couplings on the soft Z 2 -breaking parameter m 2 12 . As we will see, m 2 12 plays a very important role in our search of sizable triple Higgs couplings. Finding a way to obtain large values for m 2 12 , still being allowed by experimental data and by the theoretical constraints, turned out to be crucial in the course of this work. This also constitutes one of the main differences between our present study and previous studies on constraints in the 2HDM, from LHC physics [52][53][54], EWPO [48,55,56], flavor physics [57] and global fits [51,[58][59][60]. The relevance of m 2 12 in the context of large triple Higgs couplings in the 2HDM type II was also studied in [61,62] (with the then available data). Furthermore, in this paper, we also explore special choices for m 2 12 in relation with other 2HDM parameters, like m H , tan β and cos(β − α). In particular, we explore the implications of the setting m 2 12 = m 2 H cos 2 α/ tan β (as considered previously, e.g., in [63]).
Our paper is organized as follows. In Sect. 2 we discuss details of the 2HDM and fix our notation. The experimental expectations for the measurement of λ hhh are briefly reviewed in Sect. 3. We discuss in Sect. 4 the theoretical and experimental constraints applied to our sampling of the 2HDMs. The numerical results are presented in Sect. 5. Here we show the maximum deviations of λ hhh from the SM that are still allowed taking into account all constraints. We also present the values that can be reached for the other triple Higgs couplings involving at least one h. Our conclusions are given in Sect. 6.

The Two Higgs Doublet Model
We assume the CP conserving 2HDM. The scalar potential of this model can be written as [8]: where Φ 1 and Φ 2 denote the two SU (2) L doublets. To avoid the occurrence of tree-level flavor changing neutral currents (FCNC), a Z 2 symmetry is imposed on the scalar potential of the model under which the scalar fields transform as: This Z 2 , however, is softly broken by the m 2 12 term 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. Tab. 1 lists the couplings for each type of fermion allowed by the Z 2 parity in four different types of 2HDM. Table 1: Allowed fermion couplings in the four types of 2HDM.
However, one can use the two minimization conditions of the potential at the vacuum to substitute the bilinears m 2 11 and m 2 22 for v and tan β: Furthermore, the couplings λ i in Eq. (1) can be replaced by the physical scalar masses and mixing angles: v 2 λ 2 = 1 sin 2 β m 2 h cos 2 α + m 2 H sin 2 α −m 2 cos 2 β , where m h ≤ m H denote the masses of the CP-even Higgs-bosons, m A , m H ± denote the masses of the physical CP-odd and charged Higgs bosons respectively and, for later convenience, we have defined a new mass squared parameterm 2 , derived from m 2 12 , given by: sin β cos β .
We will study the 2HDM in the physical basis, where the free parameters of the model, which we use as input, are chosen as: 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 ∼ 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. It is convenient to express the couplings of the neutral scalar mass eigenstates h i normalized to the corresponding SM couplings. We therefore introduce the coupling coefficients c h i V V such that the couplings to the massive vector bosons are given by: where g is the SU (2) L gauge coupling, c w the cosine of weak mixing angle, c w = m W /m Z , s w = 1 − c 2 w , and m W and m Z the masses of the W boson and the Z boson, respectively. For the CP-even boson couplings we have that c hV V = s β−α and c HV V = c β−α whereas the CP-odd is c AV V = 0.
In the Yukawa sector, the discrete Z 2 symmetry leads to the following Lagrangian: where the coefficients ξ f h i are defined in Tab. 2 for type I and II. The parameters ξ f h,H can be interpreted as the ratio of the Higgs coupling with the fermions w.r.t. the SM coupling.
The potential of the 2HDM produces new interactions in the scalar sector. In this paper we will study in detail the couplings of the lightest CP-even Higgs boson with the other BSM type I type II ξ u h s β−α + c β−α cot β s β−α + c β−α cot β ξ d,l h s β−α + c β−α cot β s β−α − c β−α tan β ξ u H c β−α − s β−α tan β c β−α − s β−α tan β ξ d,l bosons, concretely λ hhh , λ hhH , λ hHH , λ hAA and λ hH + H − . We define these λ hh i h j couplings such that the Feynman rules are given by: 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 > where n is the number of identical particles in the vertex. The explicit expressions for the couplings λ hh i h j are shown in Appendix A. We adopt this notation so the light Higgs trilinear has the same definition as in the SM, i.e. −6ivλ SM with λ SM = m 2 h /2v 2 0.13. It should be noted that all 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. This limit is very interesting because, as we will discuss in Sect. 4.3, the Higgs measurements in colliders seem to overall agree with the SM values. However, in the alignment limit in general one can still have BSM physics related to the Higgs sector, like hH + H − or ZHA interactions for example. On the other hand, the parameter m 2 12 may have a relevant impact on the triple Higgs boson couplings. In the alignment limit, it does not affect the couplings λ hhh and λ hhH , but there are potentially relevant effects on the other couplings, λ hHH , λ hAA and λ hH + H − . Outside the alignment limit (i.e. for |s β−α | < ∼ 1) the effect of m 2 12 can also enter in a relevant way into λ hhh and λ hhH .

Experimental expectations for λ hhh
A determination of λ hhh (at different degrees of precision) will be able at future collider experiments. Various production cross sections show different dependences on λ hhh , making those channels complementary to each other. Most evaluations of the anticipated experimental precision in λ hhh focus on the SM value. However, as we will analyze below, substantially different values of λ hhh are possible in the 2HDM (and other BSM models). The potential for the measurement of λ hhh at a future collider experiment thus strongly depends on the value of κ λ := λ hhh /λ SM that is realized in nature.
In Fig. 1 we show the the various double Higgs production cross sections in the SM in pp collisions with √ s = 14 TeV at (next-to) leading order ((N)LO) QCD, see Ref. [64] for details. The largest cross section is given by gg → hh, 2 which will be most relevant for the measurement of λ hhh at the HL-LHC. One can see that the production cross section has a minimum around κ λ ∼ 2. Consequently, if such a value was realized, it is expected (see below) that the future experimental precision would be worse than for, e.g., κ λ = 1, where a determination at the level of ∼ 50% is anticipated [65]. Largest production cross sections, on the other hand, are for negative κ λ . Consequently, a BSM model with very small or even negative values of κ λ is expected to result in a better determination of λ hhh . A similar behavior is observed for the second largest production channel, the WBF channel pp → hhjj (where j denotes a jet), with a minimum around κ λ ∼ 1.5. Different dependences are observed for the other, less relevant channels.
Similarly, in Fig. 2 we show the dependence on δκ λ := κ λ − 1 for the Higgs-strahlung process, e + e − → Zhh (left) and the weak-boson fusion (WBF) channel, e + e − → ννhh (right) for various center-of-mass energies, √ s, at the ILC and CLIC [66]. Also indicated as horizontal colored bands are the anticipated experimental accuracies at the ILC500 (left) and ILC 1TeV, CLIC 1.4TeV and CLIC 3TeV (right). As for the HL-LHC, also at e + e − colliders the different production channels exhibit a different dependence on λ hhh . For the Higgs-strahlung process smaller (larger) cross sections are obtained for smaller (larger) κ λ . Higher values of √ s yield a weaker dependence on λ hhh , as well as a smaller absolute cross  Figure 2: Higgs-strahlung (left) and WBF production (right) of a pair of SM Higgs bosons as a function of λ hhh at the ILC and CLIC [66]. It should be noted that the experimental precision on the total cross section indicated by the horizontal bands is valid only for the SM case.
section (as typical for s-channel processes). Consequently, a determination of λ hhh based (only) on the Higgs-strahlung channel is expected to be best at lower √ s (e.g. at the ILC500) and for larger values of κ λ . The WBF channel exhibits a minimum at δκ λ ∼ 0.5. As for the Higgs-strahlung channel the dependence becomes weaker for larger values of √ s, whereas the absolute values of the cross section increase with √ s (as typical for t-channel processes).
Consequently, a case-by-case study is necessary to take into account the different, opposing effects.
The results of such a case-by-case study are shown in Fig. 3 [67]. Depicted are the relative (left) and absolute (right) accuracies of a determination of λ hhh ("λ meas /λ true ") as a function of κ λ ("λ true /λ SM ") in the range of -0.5 . . . 2. Compared are the anticipated HL-LHC precision (based on a scaling of the results for κ λ = 1, the ILC500 precision (i.e. using only the Higgs-strahlung channel) and the ILC500+ 1 TeV accuracy (i.e. also using the WBF channel results). Here it should be kept in mind that the HL-LHC analysis assumes that the other Higgs-boson couplings take their SM value, whereas for the ILC analysis it has been shown that the inclusion of the variation of the other Higgs-boson couplings does not lead to a degradation of the anticipated precision. It is worth mentioning that in all these analyses, other possible channels that might contribute to double Higgs production in extensions of the SM, like for instance the 2HDM, are not considered.
The achievable precisions follow the cross section dependences discussed above. At the HL-LHC the most (im)precise determination is expected for smaller (larger) values of κ λ . A ∼ 35(70)% relative precision is anticipated for κ λ = −0.5(2.0), as can be seen in the left plot of Fig. 3. Using only ILC500 results better (worse) experimental determinations are expected for larger (smaller) values of κ λ , ranging from ∼ 65% at κ λ = −0.5 to ∼ 15% at κ λ = 2.0. The large relative uncertainties close to κ λ = 0 are caused exactly by the smallness of the triple Higgs coupling. As can be seen in the right plot of Fig. 3, the absolute determination of κ λ continuously improves with smaller κ λ . The combination with the WBF measurements at √ s = 1 TeV yields a substantially better determination for all values of κ λ , but no monotonous behavior is found, owing to the different opposing effects, as discussed above. Future precisions between ∼ 5% and ∼ 30% are expected, depending on the value of κ λ realized in nature. Again the largest relative uncertainties of up to 30% are found close to κ λ = 0, whereas the absolute determination exhibits a nearly constant very precise determination of κ λ in the interval [−0.5, 1.0]. These results clearly show that the physics potential of a future collider experiment strongly depends on the actual value of λ hhh realized in a BSM model. This motivates the analysis presented in the following sections showing which values of λ hhh (and other triple Higgs couplings) can be realized in 2HDMs, taking into account all existing experimental and theoretical constraints.

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

Constraints from electroweak precision data
Constraints from the electroweak precision observables (EWPO) can, in a simple approximation, be expressed in terms of the oblique parameters S, T and U [46,47]. This approximation holds if the BSM effects enter mainly via corrections to gauge boson self-energies, as it is the case for extended Higgs sectors. Under these assumptions, the corrections are independent of the Yukawa sector of the 2HDM, and therefore the same for all types.
In 2HDMs there is a strong correlation between T and U , and it is known that T is by far more constraining than U [49]. Hence, U can safely be dropped in the present analysis. Specifically, our criterion to accept a point in the 2HDM parameter space, as being in agreement with the EWPO data, is as follows. For a given choice of input parameters in Eq. (13) to be allowed by the experimental observation, we require that the prediction of the S and the T parameter are in agreement with their experimental values S = 0.02 ± 0.10 and T = 0.07 ± 0.12 [45]. In this section we will study and compare the requirement of agreement at the 1 σ and 2 σ level. In our posterior numerical analysis in Sect. 5 we will require agreement at 2 σ.
In the 2HDM, as mentioned above, the most constraining oblique parameter is T , thus, we will focus in the following of this section on the constraints from the T parameter. In the forthcoming analysis in Sect. 5 we have checked that once the allowed regions by T are set, these are also allowed by S and U , i.e. effectively it is sufficient to require agreement of T with its experimental value. One peculiarity of the T parameter in the 2HDM is that it depends on the relative mass squared differences of the scalar Higgs bosons. This can be seen in the explicit expression for the T parameter in the CP conserving 2HDM that is given by [48]: where F (x, y) = x+y 2 − xy x+y log x y , and it satisfies that F (x, x) = 0. Therefore, the contributions to T become small when either the mass of H or A is sufficiently close to the mass of the charged Higgs boson H ± [55,56]. This motivates us to define three different simplified scenarios to explore the parameter space that is allowed by the EWPO in the 2HDM: scenario A, where m A = m H ± ; scenario B, where m H = m H ± and scenario C where the masses of all the BSM Higgs bosons are equal, m H = m A = m H ± . One can see from Eq. (17) that in scenario A the main contributions to the T parameter vanish for any value of c β−α , whereas in scenario B a contribution proportional to c 2 β−α F m 2 A , m 2 H ± still survives, that will remain small close to the alignment limit.
Our study of the impact of these scenarios in the prediction for the T parameter for different values of c β−α and m H ± is summarized in Fig. 4. The 2HDM parameter space is explored with the 2HDMC code [50]. For a given set of input parameters and a given Yukawa type of the 2HDM, the code computes as output the mass spectrum, decay widths and branching ratios of all the Higgs bosons. It furthermore calculates the S, T and U parameters and contribution to the anomalous magnetic moment of the muon, (g − 2) µ . In Fig. 4 it can be seen that for scenario A any mass splitting between m H and m A = m H ± is allowed inside the 2σ region even far from the alignment limit. However, this is not the case for scenario B, where it can be seen that the prediction for T is only inside the 2σ region close to the alignment limit. If one goes to higher values of c β−α (plots on the right) there are some values of m A that are disallowed, for example when m H ± = 650 GeV and c β−α = 0.25 the allowed region is m A − m H ± < 350 GeV (upper right plot). This effect becomes stronger for larger values of m H ± . For instance, for m H ± = 1000 GeV and c β−α = 0.25 (lower left plot) the allowed region shrinks to −380 GeV < m A − m H ± < 200 GeV. In general, scenario A and C (as a subset of scenario A) is broadly allowed by T for any value of c β−α and mass splitting among the Higgs bosons, whereas scenario B can lead to a large deviation if c β−α and m H ± increases (which is taken into account in Sect. 5 as discussed above).

Theoretical constraints
Like all models with extended scalar sectors, the 2HDM also faces important constraints coming from tree-level perturbartive unitarity and stability of the vacuum. We briefly describe these constraints below (for a discussion of higher-order effects and other considerations regarding the alignment limit, see, e.g., [68,69]).
• Tree-level perturbative unitarity Perturbative unitarity is achieved by demanding that the eigenvalues of the lowest partial wave scattering matrices of the 2 → 2 processes in the scalar sector of the 2HDM, at the tree level, remain below 16π. This leads to the following constraints [11,12]: It should be noted that the above requirement of tree level perturbative unitarity, limiting the maximum size of the given combinations of λ i 's, also ensures indirectly that the potential remains perturbative up to very high scales. Hence, in the present paper we do not incorporate additional constraints from other alternative criteria to require perturbativity that are based on limiting the size of the separate λ i 's which could be a priori more restrictive than the one applied here.

• Stability
First, we require the boundedness from below criterion. This criterion demands that the potential does not go to minus infinity when the field values approach infinity. This is fulfilled if the following conditions are satisfied [10,11,13]: Besides those inequalities, we will also demand that the minimum of the theory is a global minimum of the potential that can be achieved if [13] According to equations (7) to (11) the size of the triple couplings λ i are closely related to the size of the masses of the Higgs bosons and m 2 12 . In general, the size of the triple Higgs couplings involving one h and two heavy Higgs bosons grow with the corresponding heavy Higgs mass and, therefore, they can be large for large heavy masses, near the TeV scale. Consequently, unitarity sets limits on the maximum allowed size of these large heavy masses, as can be seen in Fig. 5. One finds that only in scenarios where the heavy masses are large but nearly degenerate that these unitary bounds can be relaxed. The parameter m 2 12 also plays an important role in that concern. The plots on the right in Fig. 5 show that by setting the value of this parameter to m 2 12 = m 2 H cos 2 α/ tan β, a diagonal corridor opens up allowing  Figure 5: Allowed areas in two selected Higgs masses of the 2HDM parameter space, delimited by the theoretical constraints from unitarity (green areas), stability (red areas), and both together (dotted areas), obtained from equations (18) to (28). The alignment limit is assumed and tan β is fixed to tan β = 1.5 for scenario A (top plots) and scenario B (bottom plots). m 2 12 is set to 0 (left plots), 100000 GeV 2 (middle plots) and m 2 12 = m 2 H cos 2 α/ tan β (right plots).
for larger values of these heavy masses above 1500 GeV and with a considerable splitting. On the other hand, m 2 12 enters with a negative sign in some of the stability conditions (Eqs. (24) - (27)) and Eq. (28) imposes m 2 12 ≥ 0. Therefore if m 2 12 is large the Higgs boson masses should be also large to compensate those negative contributions. In fact, setting m 2 12 to large values reduces considerably the allowed region by stability and shrinks it to the upper right corner in these two dimensional mass plots. This reduces as well the intersection area with the unitarity allowed region (dotted areas), as can be seen in the two plots in the middle with m 2 12 = 100000 GeV 2 . Here λ 1 plays an important role, as it contains a negative contribution ∝ m 2 12 that grows with tan β, see Eq. (7). This can drive λ 1 to negative values and yield disagreement with the stability condition in Eq. (24). One way to minimize this effect on λ 1 is to fix m 2 12 such that the two last terms in Eq. (7) cancel each other. This condition leads to the above commented equation allowing for the diagonal corridor in the right plots of Fig. 5 where the intersection region (dotted area) is clearly expanded. Therefore, to enlarge the allowed region by unitarity and stability in our forthcoming analysis we will consider this as a special interesting case where to explore the maximum allowed size of the triple Higgs couplings. This condition on m 2 12 has been considered previously [63] and can also be translated into a condition onm 2 , using Eq. (12), Regarding the comparison of the allowed regions for the two considered scenarios A and B, we show in Fig. 5 some specific examples, for tan β = 1.5, where one can clearly see the impact of m 2 12 = 0 and compare it with imposing Eq. (29). In the case when m 2 12 = 0 (left) all masses are allowed by stability but they are restricted by unitarity, and the final allowed dotted region is, in scenario A, for masses m H ± = m A 1000 GeV and m H 650 GeV and, in scenario B, for masses roughly below 750 GeV. When m 2 12 increases (center) the allowed region by unitarity is similar to the previous situation, but due to the large value for m 2 12 , now to get stability, the masses should be larger than approximately 500 GeV in both scenarios A and B. The situation is completely different in the right plots where Eq. (29) is adopted. In these cases masses can get very large values as well as m 2 12 and also splitting between the two free masses is allowed. This splitting stretches in both scenarios A and B as the masses grow and the final allowed region by stability and unitarity is confined to a diagonal corridor which is narrower in scenario B than in scenario A. It should be noted that in cases where Eq. (29) is satisfied, in order to cope with the theoretical constraints scenario A demands that m H ± = m A ≥ m H and scenario B that m A ≥ m H ± = m H . The allowed region by both theoretical constraints in the left and center columns would dramatically shrink for a larger value of tan β because of the size of λ 1 , but the right plots would remain similar. In some sense, Eq. (29) gives an upper limit for m 2 12 for large masses and large tan β.

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.3.2 [14][15][16][17], 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 the 2HDMC code (see Sect. 4.1). In Fig. 6 plotted in blue are shown the allowed regions of the 2HDM in the (c β−α , tan β) plane for the case where all the masses of the heavy Higgs bosons are set to 650 GeV, i.e. in the simplest scenario C. In the upper (lower) row we show the results for the 2HDM type I (II) with m 2 12 = 0, 100000 GeV and set via Eq. (29) in the left, middle and right column, respectively. The particular exclusion channel that sets a bound limiting this blue region is specified with a Latin letter and corresponds to one of the following channels: In broad terms, the 2HDM type I seems to be less constrained than type II by the searches of heavy Higgs bosons. Both types have a lower bound on tan β ∼ 1.4 from channel (d), whereas type II also has an upper bound given by channel (h). In type I negative values of c β−α are constrained by channels (a), (b) and (c) while for positive values the more relevant channel is (e). On the other hand, in type II for a negative c β−α channels (e) and (g) become the most restrictive ones, and in the positive c β−α region channel (i) is the most sensitive one.
It is also worth to notice that for m 2 12 = 100000 GeV 2 (center plots) there are more stringent bounds than in the other cases coming from channel (f) in type I and from channel (g) in type II. This is an example of how m 2 12 can be relevant in some situations when the contributions from the scalar sector are important. Clearly, the experimental bounds on BSM Higgs searches strongly depend on the masses of such particles, so the allowed contours and the exclusion channels shown in Fig. 6 will change for a different value of the masses. In general, for smaller values of the input masses the parameter space would be more constrained.

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 [1][2][3]. 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.2.3 [27,28]. 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 the 2HDMC code. The complete list of implemented experimental data can be found in Ref. [29]. Here and in our posterior analysis 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 = 43.6. In Fig. 6 we present the results of HiggsSignals for scenario C with m H = m A = m H ± = 650 GeV as a function of tan β and c β−α , which are the most relevant parameters to determine the couplings of the h boson to the SM particles. In yellow are shown the allowed regions from HiggsSignals. (In blue are shown the allowed regions from HiggsBounds, as discussed in the previous subsection). In this figure we show the contours from HiggsSignals corresponding to a 1σ (dashed lines) and 2σ (solid lines) distance from the SM fit and the contours that have the same fit as the SM (dotted grey lines). In consequence, the regions inside these dotted grey lines have a better agreement with the experimental results that the SM. It can be seen that the parameter space is strongly constrained for c β−α to be close to the alignment limit, such that h behaves sufficiently SM-like. In particular, the 2σ allowed region for the Yukawa type II (bottom) is substantially smaller compared to type I (top). In particular for type II, we find that negative values of c β−α are very disfavored. The maximum deviation from the alignment limit takes place for tan β ∼ 1, where values between c β−α = 0.13 and c β−α = −0.03 can be found inside the 2σ region from the SM. However, as tan β increases the model is forced to be very close to the alignment limit to agree with the experimental data. This is caused by an enhancement of the coupling of h to b-quark (see Tab. 2). It should be noted that in the type II fits a new allowed branch appears in the upper right part of the plot which corresponds to ξ d h = −1, known as the wrong sign Yukawa region. For type I the constraints are weaker, specially for tan β > 3, where we can accommodate inside the 2σ region values for c β−α up to ±0.3. Fig. 6 also captures the role of m 2 12 in the fits. In type I m 2 12 barely changes the fits for tan β 3 region. However, the increment of m 2 12 narrows the 1σ, 2σ contours around the alignment limit, notably for m 2 12 = 100000 GeV 2 (upper center) where the fit forces c β−α ∼ 0 when tan β is large. In the case of type II the fits seems to be roughly independent of m 2 12 , except again for m 2 12 = 100000 GeV 2 where the model is completely outside the 2σ region for tan β > 25.
In addition, it can be seen that an extensive region exists for both types that gives a better fit to the experimental data than the SM i.e. χ 2 < χ 2 SM , even though for type I m 2 12 is required to be different from zero. Such regions are expected due to the additional freedom in the 2HDM to accommodate the LHC measurements. For the sake of completeness, we would like to comment that the impact of m H , m A and m H ± could be important for the fit when they are low, because only then they could give sizable contributions to the light Higgs measurements, specially for the H boson.
Other recent studies from LHC data analysis [52][53][54], also set similar constraints on the (c β−α , tan β) plane, since these are the most relevant 2HDM parameters (entering the Higgsboson couplings) at the LHC. One of the main differences to our study is that, as emphasized in the introduction, we have a strong focus on the role played by the m 2 12 parameter, which turns out to be relevant in our search of sizable triple Higgs couplings.

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 [57,58]. Here we will take into account the decays B → X s γ and B s → µ + µ − , which we find to be the most constraining ones and whose experimental values are (we use the average from [45]): We will set our bounds in the 2σ region from the central value according to the experimental value.
In order to compute the theoretical predictions in the 2HDM we have used the public code SuperIso [30,31] with the model input given by 2HDMC. Moreover, we have included in SuperIso the contributions to the Wilson coefficient C P from the Higgs-penguin diagrams, that are missing in the public version and that can be relevant for the BR(B s → µ + µ − ) prediction [32][33][34].
In Fig. 7 we present the allowed regions from the flavor constraints in the (m H ± , tan β) plane in the alignment limit for scenario C (all masses of BSM bosons degenerated) for Yukawa types I (upper row) and II (lower row) for m 2 12 = 0, 100000 GeV and set via Eq. (29) in the left, middle and right column, respectively. We show the regions allowed by B → X s γ (pink areas) and by B s → µ + µ − (teal areas). Dotted areas are the intersections of these two allowed regions. The 2HDM contribution to the process B → X s γ depends on the couplings of the b and s quarks with the other u-type quarks through a charged Higgs boson. As the Yukawa couplings of the charged Higgs bosons in the 2HDM scale like the ones of the CP-odd Higgs boson, this coupling is given by a combination of ξ u,d A (see Tab. 2) and the quark masses. In the case of model type I those couplings are enhanced for large values of cot β, in consequence the region of low tan β is forbidden in the top plots of Fig. 7 and softly fades as the mass of the charged Higgs increases. On the contrary, in type II it is found a well known tan β independent constraint of m H ± > 500 GeV. The BSM contributions to B → X s γ are induced from the Yukawa coupling and therefore neither c β−α or m 2 12 affects the bounds, as it can be seen in the figure. Focusing on B s → µ + µ − one finds a similar constraint for low tan β on both model types due to analogous arguments discussed before for B → X s γ. Nevertheless, in model type II there is a disallowed region for large tan β and low  masses. This is due to the contributions from the Higgs-penguin diagrams (mediated by H and h) to the process B s → µ + µ − which are sensitive to m 2 12 , via λ HH + H − and λ hH + H − from the loops involving charged Higgs bosons, and that are enhanced at large tan β (see also [33]). The largest effect from m 2 12 on B s → µ + µ − is from λ HH + H − since the H-penguin diagram goes as tan 3 β, and this leads to relevant constraints in the large tan β and low m H + region. If, however, m 2 12 is fixed to Eq. (29) and if the alignment limit is taken, then the coupling λ HH + H − vanishes and in consequence the Higgs penguins contributions are not large enough to give a bound in that region.

Numerical results
In this section we analyze numerically which intervals (or extreme values) of the various triple Higgs boson couplings are still allowed, taking into account all experimental and theoretical constraints as discussed in Sect. 4. In the case of λ hhh this will give a guideline to which collider option may be needed to perform a precise experimental determination. For the triple Higgs couplings involving heavy Higgs bosons this will indicate in which processes large effects, e.g. possibly enhanced production cross sections, can be expected due to large triple Higgs couplings.
We perform our evaluation in both type I and type II models (and leave the other types for future investigations). We start our exploration with the "simplest" scenario C, but later also explore scenario A and B. In the headlines of our plots we indicate which type and which scenario are chosen. The other parameters are chosen such as to maximize either the deviations of λ hhh from it SM value (where the plots below show κ λ := λ hhh /λ SM ), or to maximize (positive or negative) the size of the triple Higgs couplings involving the heavy Higgs bosons (where the plots below show the triple Higgs couplings as defined in Eq. (16)).

Scenario C
We start with scenario C, i.e. m H ± = m H = m A 3 , and m h = 125 GeV. In Fig. 8 we show the (c β−α , tan β) plane in the 2HDM type I, where m 2 12 is fixed by Eq. (29) to maximize the regions allowed by unitarity and stability of the potential, see Sect In the (c β−α ,tan β) plane this intersection defining the total allowed area starts at tan β ∼ 2 up to the highest investigated values, where we stopped at tan β = 50. c β−α = 0, is allowed for all tan β values, with a roughly triangular shape, extending up to c β−α ∼ 0.2.
The results for κ λ = λ hhh /λ SM are presented in the lower plot of Fig. 8(A), with the total allowed area discussed above being now marked by the bounding black solid line. The red solid line indicates κ λ ≡ 1. This is either the alignment limit for c β−α = 0, or the "wrong sign limit" in the upper right corner. For the latter, see the discussion in Sect. 4.4. The color code shows the values reached by κ λ . In the area allowed by all experimental and theoretical constraints, values of κ λ < ∼ 1 are realized, going down to κ λ ∼ −0.4 in the "tip" to the right of the allowed area. The corresponding implications will be discussed in Sect. 5.4.
We now turn to the triple Higgs couplings involving at least one heavy Higgs boson. In Fig. 8(B) we show the results for λ hhH , λ hHH , λ hAA and λ hH + H − in the upper left, upper right,  lower left and lower right plot, respectively. As before, the area allowed by all experimental and theoretical constraints is indicated by a black solid line, and the color code shows the values reached by the triple Higgs couplings. In all four cases we find positive couplings with the minimum values reached for c β−α = 0. The larger values are found in the right edge of the allowed area, with largest values (as in the case of λ hhh ) in the "tip" to the right of the allowed area. λ hhH is found to be larger around tan β ∼ 8 and c β−α ∼ 0.1. The maximum values found for the rest of the triple Higgs couplings in this case are λ hHH ∼ 12, λ hAA ∼ 12 and λ hH + H − ∼ 24. It should be noted that here and in the following λ hH + H − always reaches the maximum values of all the considered triple Higgs boson couplings. The corresponding phenomenological implications will be discussed in Sect. 5.4.
We continue the exploration of scenario C, type I in the (c β−α , m 2 12 ) plane for m H ± = m H = m A = 650 GeV and tan β = 7.5, as shown in Fig. 9. The sequence and the color coding of the plots is the same as in Fig. 8. The overall allowed area is restricted, particularly by the requirement of unitarity and stability, to be within a curved band around m 2 12 = 55000 GeV 2 , ranging from c β−α ∼ 0 to c β−α ∼ 0.28. Here the purple solid line in the middle left plot  indicates that Eq. (29) is satisfied. The lower plot in Fig. 9(A) presents the results for κ λ , which show a weak dependence on m 2 12 . Values of κ λ ∼ 1 are found around c β−α = 0 (as required by the alignment limit), but also around c β−α ∼ 0.26. The lowest value of κ λ ∼ 0.5 is realized for c β−α = 0.2, whereas the highest value of κ λ ∼ 1.2 are found for c β−α ∼ 0.28. Contrary to the (c β−α , tan β) plane shown in Fig. 8, we now also encounter values of κ λ larger than 1. However, these are realized for the largest departure of the alignment limit, and thus will be under scrutiny by the next round of Higgs-boson rate measurements at the LHC.
The results for the triple Higgs couplings involving heavy Higgs bosons are shown in Fig. 9(B), analogous to Fig. 8 5,14]. It should be noted that due to the contribution from m 2 12 here these couplings can also be slightly negative. As before, the maximum of λ hhH is found for c β−α ∼ 0.1 whereas for the other couplings, which can be of O(10), the largest values are realized for the largest departure of the alignment limit, and thus will be under scrutiny by the next round of Higgs-boson rate measurements at the LHC.
We finish our analysis of the scenario C, type I in Fig. 10  The values that can be reached by κ λ , as shown in the lower plot of Fig. 10(A), range from κ λ ∼ 0.07 for c β−α ∼ 0.1 and large m H ± close to 1200 GeV to about κ λ ∼ 1.2 for the largest allowed c β−α values and m H ± ∼ 300 GeV. The ranges reached by the triple Higgs couplings involving at least one heavy Higgs boson, as shown in Fig. 10  the edge for larger c β−α and m H ± > ∼ 800 GeV. We finish our analysis of scenario C with the (c β−α , m 2 12 ) plane in the 2HDM type II for m H = m A = m H ± = 1000 GeV and tan β = 0.9, as presented in Fig. 11. The sequence of the plots and the color coding are as in Fig. 9. The total allowed area is found, roughly between c β−α ∼ −0.05 and c β−α < ∼ 0.1, as well as m 2 12 > ∼ 2 · 10 5 GeV 2 and m 2 12 < ∼ 6 · 10 5 GeV 2 . The values that can be reached by κ λ , as shown in the lower plot of Fig. 11(A), range from κ λ ∼ −1.0 for c β−α ∼ 0.13 and low m 2 12 to κ λ = 1 for the alignment limit. The ranges reached by the triple Higgs couplings involving at least one heavy Higgs boson, as shown in Fig. 11 values, and are nearly independent on c β−α . In contrast, λ hhH shows dependence on both variables where its maximum is found around m 2 12 ∼ 400000 GeV 2 and c β−α ∼ 0.08 and its minimum is found around m 2 12 ∼ 550000 GeV 2 and c β−α ∼ −0.03. As for the 2HDM type I, the phenomenological interpretation of these intervals will be given in Sect. 5.4.

Scenario A
We continue our numerical investigation by relaxing the conditions for the heavy Higgs-boson masses and evaluate the triple Higgs-boson couplings in scenario A, as defined in Sect. 4.1, m A = m H ± = m H 4 , and m h = 125 GeV. In Fig. 12 we present the (m H ± = m A , m H ) plane with m 2 12 = (m 2 H cos 2 α)/(tan β), to maximize the parameter space allowed by unitarity and stability of the Higgs potential, and c β−α = 0.2 and tan β = 10. The upper two rows show the various constraints, with the same color coding as in Fig. 8. One can see that the LHC searches and measurements, as well as the flavor observables allow for the whole plane. Unitarity and stability roughly select a square bounded from above by m H ± ∼ m H ∼ 1000 GeV. The results for λ hhh are not explicitly shown, as they vary only very weakly in the chosen scenario. The values reached are in the interval κ λ ∼ [0.98, 1.02]. The lower two rows in Fig. 12 show the results for the triple Higgs couplings involving at least one heavy Higgs boson. The upper left plot (of the two lower rows) shows λ hhH , which is independent of m H ± . Lowest (highest) values are reached for high (low) values of m H , following the analytic result in Eq. (33). They range from 0.02 to -1.5.
The upper right plot depicts the results for λ hHH , again independent of m H ± . Here lowest (highest) values are reached for low (high) values of m H , following the analytic result in Eq. (35). For λ hHH they range from 0.2 to 16. The lower row shows the results for λ hAA (left) and λ hH + H − (right), which exhibit a similar behavior, see Eq. (37) and Eq. (39). The values are nearly independent of m H , where lowest (highest) values are found for low (high) m A = m H ± . They range from 0 to 16 for λ hAA , and from 0 to 32 for λ hH + H − . As in Sect. 5.1 we leave the phenomenological discussion to Sect. 5.4.
Analogous results in the 2HDM type II are presented in Fig. 13, with the color codings as in Fig. 12. As before m 2 12 is fixed by m 2 12 = (m 2 H cos 2 α)/(tan β). In order to maximize the results for the triple Higgs couplings we have chosen c β−α = 0.025 and tan β = 6.5. The overall allowed region, as depicted in the upper two rows, can be found on the strip roughly around the diagonal m H ± = m A ∼ m H . The results for λ hhh again vary only weakly in this region, and are found in the interval κ λ ∼ [0. 8,1]. The third row shows the results for λ hhH (left) and λ hHH (right), which follow similar patterns and are independent of m H ± , see Eq. (33) and Eq. (35). Lowest (highest) values are found at low (high) m H , ranging from 0 to 1.25 for λ hhH and from 0.15 to 3 for λ hHH . The fourth row presents the results for λ hAA (left) and λ hH + H − (right), which again follow similar patterns and are nearly independent of m H , see Eq. We finish our analysis in the scenario A with the 2HDM type II presented in Fig. 14, with the color codings as in Fig. 12. In comparison with the previous analysis we have chosen a relatively low value of tan β = 0.9, and fixed m 2 12 = 100000 GeV 2 , while for c β−α a relatively large value (for the 2HDM type II) of c β−α = 0.05 was chosen. The overall allowed region, as depicted in the upper two rows can be found roughly around 800 GeV < m H ± = m A < 1100 GeV and m H between 500 GeV and 1000 GeV. As before,   the results for λ hhh vary only weakly in this region, and it takes values for κ λ ∼ 0.9 in the whole plane. The third row shows the results for λ hhH (left) and λ hHH (right), and as before both are independent of m H ± , see Eq. (33) and Eq. (35). λ hhH exhibits a small variation between -0.14 to 0.23. λ hHH , on the other hand, can reach very large values for large m H , and it is found to be in the range of 0.3 and 15. The fourth row presents the results for λ hAA (left) and λ hH + H − (right), which again follow similar patterns and are independent of m H , see Eq. (37) and Eq. (39). The lowest (highest) values are found at the lowest (highest) allowed values for m H ± = m A ∼ 800 (1100) GeV. They range from 8 to 16 for λ hAA and from 16 to 32 for λ hH + H − . The phenomenological implications are discussed in Sect. 5.4.

Scenario B
We finish our numerical investigation with the third scenario suggested by the electroweak precision observables, scenario B, as defined in Sect. 4.1, m A = m H ± = m H 5 , and m h = 125 GeV.
In Fig. 15 we present the (m H ± = m H , m A ) plane with the other parameters chosen as in the corresponding scenario A, with m 2 12 = (m 2 H cos 2 α)/(tan β) to maximize the parameter space allowed by unitarity and stability of the Higgs potential, and for c β−α = 0.2 and tan β = 10. The upper two rows show the various constraints, with the same color coding as in Fig. 8. Besides, in the upper left plot, where we indicate the regions allowed by the LHC measurements, we also indicate the bound arising from the EWPO, see Fig. 4. While the LHC measurements of the SM-like Higgs boson as well as the direct searches for BSM Higgs bosons do not yield restrictions in the parameter space, the EWPO favor a broad region roughly around the diagonal m H ± = m H ∼ m A . As in the corresponding scenario A, unitarity and stability roughly select a square bounded from above by m A ∼ m H ∼ 1000 GeV. The results for λ hhh are again not explicitly shown, as they vary only very weakly in the chosen scenario. The values reached are in the interval κ λ ∼ [0.98, 1.03]. The lower two rows in Fig. 15 show the results for the triple Higgs couplings involving at least one heavy Higgs boson. The upper left plot (of the two lower rows) shows λ hhH , which is independent of m A . Lower (higher) values are reached for high (low) values of m H ± , following the analytic result in Eq. (33). They range from -1.5 to 1. The upper right plot depicts the results for λ hHH , again independent of m A . Here, the lowest (highest) values are reached for the lowest (highest) allowed values of m H ± , following the analytic result in Eq. (35). For λ hHH they range from 0.2 to 15. The lower row shows the results for λ hAA (left) and λ hH + H − (right), where the latter exhibits a similar behavior as λ hHH , see Eq. (37) and Eq. (39). The values of λ hAA (λ hH + H − ) are nearly independent of m H ± (m A ), where the lowest (highest) values are found for the lowest (highest) allowed values of m A (m H ± = m H ). They range from 0.2 to 16 for λ hAA , and from 0.5 to 30 for λ hH + H − . As in Sect. 5.1, we leave the phenomenological discussion to Sect. 5.4.
The final scenario analyzed is scenario B analogous to the last example in scenario A as presented in Fig. 15, with the color codings as in Fig. 12. As in scenario A we have chosen a relatively low value of tan β = 0.9, and fixed m 2 12 = 100000 GeV 2 , while for c β−α a relatively large value (for the 2HDM type II) of c β−α = 0.05 was chosen. The overall allowed region, as depicted in the upper two rows can be found roughly around m H ± = m H ∼ 900 GeV and m A between 500 GeV and 1000 GeV, analogous to the corresponding scenario A. As before, the results for λ hhh vary only weakly in this region, and it leads to values of κ λ ∼ 0.9 in the whole plane. The lower two rows in Fig. 16 show the results for the triple Higgs couplings involving at least one heavy Higgs boson. The upper left plot (of the two lower rows) shows λ hhH , which is independent of m A and varies only weakly in the parameter plane. Lower (higher) values are reached for high (low) values of m H ± , following the analytic result in Eq. (33). They range from -0.14 to 0. The upper right plot depicts the results for λ hHH , again independent of m A . Here, the lowest (highest) values are reached for the lowest (highest) allowed values of m H ± , following the analytic result in Eq. (35). For λ hHH they range from 8 to 15. The lower row shows the results for λ hAA (left) and λ hH + H − (right), where the latter exhibits a similar behavior as λ hHH , see Eq. (37) and Eq. (39). The values of λ hAA (λ hH + H − ) are independent of m H ± (m A ), where the lowest (highest) values are found for the lowest (highest) allowed values of m A (m H ± = m H ). They range from 0.14 to 16 for λ hAA , and from 15 to 30 for λ hH + H − .
Finally, to close the numerical results section, we present in Tab. 3 some examples of interesting configurations that maximize the size of the triple Higgs couplings. In the Yukawa type I, all examples have tan β > 1 while for type II all the points are around tan β ∼ 1. This is mainly due to the constraints from the LHC data, because it is easier to accommodate a SM-like Higgs in those regions (see Sect. 4.3). In addition, particularly in type I, flavor observables disallow low values of tan β. In type I we recover a larger allowed parameter region by choosing m 2 12 according to Eq. (29), especially for the larger values of tan β, which are easier in conflict with the theoretical constraints, as we discussed in Sect. 4.2. Furthermore, in type II the tight constraint from B → X s γ that sets m H ± > ∼ 500 GeV should be also satisfied. However, as we have discussed in the previous subsections, the main constraint that prevent from obtaining large triple Higgs couplings are the theoretical constraints.
For both types, I and II, points with larger triple Higgs couplings are also the ones with the heavier Higgs masses around 1 TeV (where we have not explored values above ∼ 1.6 TeV). The only exception to this are λ hhh and λ hhH . In the case of λ hhh in the alignment limit the SM value is reproduced. On the other hand, λ hhH is proportional to c β−α , see Appendix A. Consequently, their extrema are both found outside the alignment limit. As a consequence those points are stronger tested and possibly "easier" excluded in the future by more precise measurements. Overall, and particularly for type II, κ λ is close to unity. While this does not correspond to large enhancements of di-Higgs production, the deviations are still large enough to be tested at future colliders, see Sect. 3. On the other hand, it is possible to find large allowed values of couplings involving more than one heavy Higgs boson near the alignment limit. Those triple Higgs couplings always have positive and large values for all scenarios A, B and C. In fact, they can be larger in scenarios A and B due to the allowed splitting of two of the masses, as we have seen in Sects. 5.2 and 5.3. In those cases, the larger mass is the one of the Higgs boson that appears in the vertex. Overall, we find values of O(10), with the maximum value corresponding to λ hH + H − ∼ 30 in both type I and II.
In Tab. 3 we also include points with smaller Higgs masses that also yield interesting sizes of the triple Higgs couplings. Due to the relatively smaller masses, these points are better kinematically accessible. These last kind of points are presumably the easiest to probe at future colliders. For these more moderate, but potentially more accessible masses, we find

Possible implications for future collider measurements
We now turn to the phenomenological implications of the allowed ranges found for the various triple Higgs couplings, as discussed in the previous subsections. As an overall result we find that the allowed intervals for the various triple Higgs couplings depend only weakly on the chosen EWPO scenario A, B or C. However, the 2HDM type I exhibits a substantially stronger variation in λ hhh than type II. This is mostly owed to the larger allowed deviation from the alignment limit, see Sects. 4.3, 4.4. In this section we will concentrate on the anticipated impact of the triple Higgs couplings on the various di-Higgs production cross sections, where we leave a full phenomenological analysis for future work [70].
For λ hhh we roughly find allowed intervals of [−0.5, 1.5] in the 2HDM type I and [0, 1] in type II. While the production of two SM-like Higgs bosons, both at pp and at e + e − colliders depends already at the tree-level on λ hhh and λ hhH , the dependence on λ hhh is expected to be substantially stronger due to the propagator suppression with the inverse of m 2 H of λ hhH . Consequently, over the possible parameter range of λ hhh the HL-LHC is not expected to yield a precision on κ λ better than 35%, and a deviation from λ hhh = 0 can not be established better than ∼ 2 σ. Comparing the HL-LHC to the ILC500, the HL-LHC performs better (worse) than the ILC500 for κ λ < ∼ ( > ∼ ) 0.5, where both intervals are still allowed in both types of the 2HDM. In other words, the HL-LHC results in comparison with the ILC500 may look a bit better than anticipated for κ λ = 1. However, in this comparison it must be kept in mind that the HL-LHC analysis is based on the variation of the Higgs triple coupling only, whereas for the ILC500 (at κ λ = 1) it has been shown that the analysis holds also for a variation of all Higgs-boson couplings within their anticipated experimental accuracies. Furthermore, deviations below κ λ ∼ 0.5 are realized for larger deviations from the alignment limit and may thus be tested in the next round of Higgs rate measurements at the LHC. Combining the ILC500 measurements with the final stage of the ILC1000, the Linear Collider shows a substantially better result than the HL-LHC for all the allowed λ hhh parameter space. Only around a vanishing trilinear Higgs coupling similar precisions are anticipated (but the above mentioned caveat of the differences in the HL-LHC and ILC analyses still holds).
The phenomenological implications of the allowed ranges for λ hhH are twofold. This coupling can enhance or suppress the contribution of the off-shell heavy Higgs in the hh production, which, however, are generally suppressed as mentioned above. On the other hand, a very large enhancement of this coupling would yield a relatively large cross section for hH production. However, we find that large values of λ hhH are not allowed taking all existing experimental and theoretical constraints into account.
The triple Higgs couplings involving two heavy Higgs bosons, λ hHH , λ hAA and λ hH + H − can have a very strong impact on the heavy di-Higgs production and possibly facilitate the discovery of such heavier Higgs bosons (see, e.g., [54,69]). Roughly independent of the EWPO scenario and the 2HDM type, we find values of up to 15, 16 and 32, respectively. Here it must be kept in mind that the larger values of a triple Higgs coupling of h to two heavy Higgs bosons are realized for larger values of the respective heavy Higgs-boson mass. Consequently, the effects of the large coupling and the heavy mass always go in opposite directions. A detailed study will be left to future work [70].

Conclusions
An important task at future colliders is the measurement of the triple Higgs coupling λ hhh . Depending on its size relative to the SM value, certain collider options result in a higher experimental accuracy. Similarly, large values of triple Higgs couplings involving heavy Higgs bosons can lead to enhanced production cross sections of BSM Higgs bosons.
Within the framework of Two Higgs Doublet Models (2HDM) type I and II we investigate the allowed ranges for all triple Higgs couplings involving at least one light, SM-like Higgs boson. We take into account all relevant theoretical and experimental constraints. From the theory side these comprise unitarity, and stability conditions. From the experimental side we require agreement with the direct BSM Higgs-boson searches, as well as with measurements of the SM-like Higgs-boson rate as measured at the LHC. We furthermore require agreement with flavor observables and electroweak precision data (where the T parameter plays the most important role). In this context we investigate more extensively the dependence of several of these constraints on the soft Z 2 -breaking parameter, m 2 12 . Here we find that large values of this parameter can affect notably the allowed parameter space, especially in the region of large tan β.
For theoretical constraints m 2 12 plays a key role: lower (higher) values are favored by the tree-level stability (unitarity) constraint, and the size of the intersection region is thus controlled by m 2 12 . Thus, to enlarge the allowed region by both unitarity and stability we have used Eq. (29) on several occasions.
Regarding the experimental constraints, BSM Higgs boson searches and measurements of the 125 GeV Higgs boson at the LHC can also be sensitive to the effects of m 2 12 in the scalar sector such like the h → γγ decay (via the hH + H − vertex) or the production of a heavy BSM boson that decays to two 125 GeV Higgs bosons, specially in the range of low masses. On the other hand, the triple Higgs couplings λ hH + H − and λ HH + H − also enter in the 2HDM prediction for B s → µ + µ − via the h and H Higgs penguins contributions with charged Higgs bosons in the loops, and they can be relevant (see also [33]). The largest effect from m 2 12 in B s → µ + µ − is due to λ HH + H − in the region of large tan β and low m H ± and, therefore, this region is correspondingly constrained by the B s → µ + µ − data.
Based on a parameter scan we investigated several mass and parameter planes. We demanded agreement with the above given constraints and evaluated the maximum and minimum values of the various triple Higgs couplings. For the SM-type triple Higgs coupling w.r.t. its SM value, κ λ = λ hhh /λ SM , we roughly find allowed intervals of [−0.5, 1.5] in the 2HDM type I and [0, 1] in type II. The production of two SM-like Higgs bosons, both at pp and at e + e − colliders depends already at the tree-level strongly on λ hhh . Consequently, over the possible parameter range of λ hhh the HL-LHC is not expected to yield a precision on κ λ better than 35%, and a deviation from λ hhh = 0 can not be established better than ∼ 2 σ. Comparing the HL-LHC to the ILC500, the HL-LHC performs better (worse) than the ILC500 for κ λ < ∼ ( > ∼ ) 0.5. Combining the ILC500 measurements with the final stage of the ILC1000, the Linear Collider shows a substantially better result than the HL-LHC for all the values of λ hhh in the allowed intervals that we have found.
The production of two light Higgs bosons can also depend on λ hhH in the 2HDM. In this case, the prediction in the alignment limit is λ hhH = 0, but here we reach the maximum (minimum) value around c β−α = ±0.05. We find that the total allowed interval of this coupling is [−1.4, 1.5] for type I and [−1.6, 1.8] for type II.
Concerning the triple Higgs couplings involving two heavy 2HDM Higgs bosons, we find large allowed values for both 2HDM type I and II. For λ hHH , λ hAA and λ hH + H − we find maximum values of up to 15, 16 and 32, respectively. These triple Higgs couplings can have a very strong impact on the heavy di-Higgs production at pp and e + e − colliders. Large coupling values can possibly facilitate the discovery of such heavier Higgs bosons. However, it must be kept in mind that the larger values of triple Higgs couplings of h with two heavy Higgs bosons are realized for larger values of the respective heavy Higgs-boson mass. Consequently, the effects of the large coupling and the heavy mass always go in opposite directions. A detailed analysis of the various production cross sections will be analyzed elsewhere [70]. λ i basis: = −iv 3λ 1 c α c β s 2 α + 3λ 2 c 2 α s α s β + (λ 3 + λ 4 + λ 5 ) c 3 α c β − 2c 2 α s α s β − 2c α c β s 2 α + s 3 α s β .
Physical basis: / a c 1 u 0 7 d + / t t w / u n 6 q s k A w n L E s y e R 6 A w o Q L n G i u E z z P J U I a J H g W L F 5 X / N k F S s U z 8 V Y v c / R T m A s e c Q b a Q L M D q z 8 N c M 5 F q f n i Y 8 6 Z L i S u W s 4 G j H A p U h A r 7 9 8 R 4 y 7 0 / J b j T I 2 n x g / U j L Q 8 j A 9 X r y 6 B n u T z W I + z i I J P u 0 G v w U G Q X S D d K g K j Y J X H 8 b u n O y 4 B J t n 7 p j J c K / t r p e O Y R U M O c w n p k 9 I A 1 T b 9 P v V C U D G G 9 d n N m f W e V b J d O K z h y n O K I t z e u l V P j W B m 7 c 5 w M K y L X m 3 c T d M h m z q Z t f 9 M w 4 w V K Q r N E l D K c 4 e 5 9 k u Q m r P E G E 4 L h T m w B c z R M 6 2 A F J V f 1 v 9 1 R R 9 H m a Q 6 R l r P l 7 U l p E o t 0 8 B o U t C x 2 u U q 8 H + c V + j o p V 9 y k R c a B T M S w 0 V F Q n V G q + d B Q y 6 R 6 W R p G m C S m y 0 p i 0 E C 0 + Y R N Z x Y l u a g x + 7 A 3 N q v U + p v g z M x u b u h X G 1 O R w P 3 + W D 0 Z t Q 5 6 m 4 C 2 y M P y S P S J S 5 5 Q Y 7 I M T k h E 8 K s T 9 Z n 6 6 v 1 z f 5 i / 7 B / 2 r / W U t v a f P O A N M r + / R d 2 e P k 7 < / l a t e x i t > λ i basis: = iv λ 1 c β s α s 2 β − λ 2 c α c 2 β s β + λ 3 c 3 β s α − c α s 3 β + (λ 4 + λ 5 ) c α c 2 β s β − c β s α s 2 β .
Physical basis: