Sensitivity to triple Higgs couplings via di-Higgs production in the 2HDM at the (HL-)LHC

An important task of the LHC is the investigation of the Higgs-boson sector. Of particular interest is the reconstruction of the Higgs potential, i.e. the measurement of the Higgs self-couplings. Based on previous analyses, within the 2-Higgs-Doublet Model (2HDM) type I, we 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 the (HL-)LHC with a center-of-mass energy of 14TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$14 \,\, \textrm{TeV}$$\end{document} at next-to-leading order in the heavy top-quark limit with the code HPAIR. We investigate in particular the process gg→hh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$gg \rightarrow hh$$\end{document}, with h being the Higgs boson discovered at the LHC with a mass of about 125GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$125 \,\, \textrm{GeV}$$\end{document}. The top box diagram of the loop-mediated gluon fusion process into Higgs pairs interferes with the s-channel exchange of the two CP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{C}\mathcal{P}}$$\end{document}-even 2HDM Higgs bosons h and H. The latter two involve the triple Higgs couplings (THCs) λhhh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{hhh}$$\end{document} and λhhH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{hhH}$$\end{document}, respectively, possibly making them accessible at the HL-LHC. Depending on the size of the involved top-Yukawa and THCs as well as on the mass of H, the contribution of the s-channel H diagram can be dominating or be highly suppressed. We find regions of the allowed parameter space in which the di-Higgs production cross section can differ by many standard deviations from its SM prediction, indicating possible access to deviations in λhhh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{hhh}$$\end{document} from the SM value λSM\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{\textrm{SM}}$$\end{document} and/or contributions involving λhhH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{hhH}$$\end{document}. The sensitivity to the beyond-the-SM (BSM) THC λhhH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{hhH}$$\end{document} is further analyzed employing the mhh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{hh}$$\end{document} distributions. We demonstrate how a possible measurement of λhhH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{hhH}$$\end{document} depends on the various experimental uncertainties. Depending on the underlying parameter space, the HL-LHC may have the option not only to detect BSM THCs, but also to provide a first rough measurement of their sizes.


Introduction
The discovery of a new scalar particle with a mass of ∼ 125 GeV by ATLAS and CMS [1][2][3] -within the experimental and theoretical uncertainties -is in agreement with the properties of the Standard Model (SM) Higgs boson.On the other hand, no conclusive sign of Higgs bosons beyond the SM (BSM) has been observed so far.However, the experimental results about the state at ∼ 125 GeV, whose couplings are known up to now to an experimental precision of roughly ∼ 10 − 20%, leave ample room for interpretations in BSM models.Many BSM models feature extended Higgs-boson sectors with correspondingly extended Higgs potentials.Consequently, one of the main tasks of present and future colliders will be to determine whether the observed scalar boson forms part of the Higgs sector of an extended model, or not.
In contrast to the Higgs couplings to the SM fermions and gauge bosons, the trilinear Higgs self-coupling λ hhh remains to be determined.So far it has been constrained by ATLAS [4] to be inside the range −0.4 < λ hhh /λ SM < 6.3 at the 95% C.L. and −1.24 < λ hhh /λ SM < 6.49 at the 95% C.L. by CMS [5], both assuming a SM-like top-Yukawa coupling of the light Higgs.Many BSM models can still induce significant deviations in the trilinear coupling λ hhh of the SM-like Higgs boson with respect to the SM value, see, e.g., Ref. [6] for an up-to-date investigation.For recent reviews on the measurement of the triple Higgs couplings at future colliders see for instance Refs.[7,8].In case a BSM Higgs sector manifests itself, it will be a prime task to measure also the BSM trilinear Higgs self-couplings.
One of the simplest extensions of the SM Higgs sector is the 2-Higgs-Doublet Model (2HDM) [9][10][11][12], where a second Higgs doublet is added to the SM Higgs sector.After electroweak symmetry breaking this leads to five physical Higgs bosons, two CP-even bosons h and H, where by convention m h < m H , one CP-odd boson A and two charged Higgs bosons H ± .The ratio of the two vacuum expectation values of the neutral components of the two Higgs doublets is defined as tan β ≡ v 2 /v 1 .To avoid flavor-changing neutral currents at tree level, a Z 2 symmetry is imposed [13], possibly softly broken by the parameter m 212 .Depending on how this symmetry is extended to the fermion sector, four types of the 2HDM can be realized: type I and II, flipped (type III) and lepton specific (type IV) [11].
In Ref. [14] an analysis was presented of the possible size of triple Higgs couplings (THCs) in the 2HDM type I and II taking into account all relevant experimental and theoretical constraints. 1 For that analysis it was assumed that the lightest 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.(An update and extension to type III and IV was presented in Ref. [15].)Future e + e − linear colliders, like the ILC [16] and CLIC [17], will play a key role for the measurement of the Higgs potential and in detecting possible deviations from the SM with high precision [7,8,[18][19][20][21].Employing the results of Ref. [14], in Ref. [22] the sensitivity of the ILC and CLIC to various 2HDM THCs (including BSM THCs) was analyzed.Further analyses of THCs at e + e − colliders were presented in Refs.[23,24].Recent reviews on triple Higgs couplings at e + e − colliders can be found in Refs.[7,8,20,21].
In this paper, based on the results of Ref. [14], we complement the above results with an analysis of the senstivity to BSM triple Higgs couplings at the LHC, and in particular its high-luminosity phase, the HL-LHC.Further analyzes involving BSM triple Higgs couplings can be found in Refs.[6,8,[25][26][27].However, while these papers took the effects of BSM THCs into account, to our knowledge no analysis for the (HL-)LHC exists attempting to quantify the potential sensitivity to BSM triple Higgs couplings.
The main Higgs pair production process at the LHC is gluon fusion into Higgs pairs [28].Here we investige in particular gg → hh in the 2HDM type I and II, with h corresponding to the state discovered at the LHC at ∼ 125 GeV.The process is loop-mediated already at leading order and consists of a triangle and a box top-loop contribution.For small values of tan β the bottom loop plays only a minor role.In the SM, the box diagram interferes destructively with the triangle contribution.In the 2HDM, we have both the h and H s-channel exchange in the triangle contribution, where a resonantly produced H with subsequent decay into hh can lead to a significantly enhanced cross section.For our analysis, we take into account the next-to-leading order QCD corrections to the process in the heavy top-quark limit [29] by making use of the accordingly modified [6,30] program HPAIR.We find regions of the allowed parameter space in which the di-Higgs production cross section can differ by several standard deviations from its SM prediction, indicating possible access to deviations in λ hhh from λ SM and/or contributions involving λ hhH .The sensitivity to λ hhH is further analyzed employing the m hh distributions in the gg → hh production cross section.We investigate how a possible measurement of λ hhH depends on the assumed experimental uncertainties in m hh , such as smearing, bin width, as well as on the position of the bins.We demonstrate that, depending on the underlying parameter space, the HL-LHC may have the potential not only to detect BSM triple Higgs couplings, but also to provide a first rough measurement of their size.
Our paper is organized as follows.In Sect. 2 we briefly review the 2HDM, fix our notation, define the benchmark planes used later for our investigation and summarize the constraints that we apply (which are the same as in Refs.[14,15]).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 in the contribution from the s-channel H exchange.In Sect.3.3 we analyze a possible sensitivity of the di-Higgs production cross section at the (HL-)LHC to λ hhh and in particular to λ hhH .Finally, in Sect. 4 we present the possible HL-LHC sensitivity to λ hhH via the m hh distribution, and in Sect. 5 also assess its dependence on smearing, bin width and position of the bins.Our conclusions are given in Sect.6.

The Model and the constraints
In this section we give a short description of the 2HDM to fix our notation.We briefly review the theoretical and experimental constraints, which are taken over from Refs.[14,15].Finally we will define the benchmark planes for our analysis of the gg → hh production cross section.

The 2HDM
We assume the CP-conserving 2HDM [9][10][11][12], where the potential can be written as, After electroweak symmetry breaking (EWSB) the two SU (2) L doublets Φ 1 and Φ 2 can be expanded around their two vacuum expectation values (vevs) v 1 and v 2 , respectively, as with the vev ratio given by tan where v 246 GeV is the SM vev.The eight (scalar) degrees of freedom, φ ± 1,2 , ρ 1,2 and η 1,2 , give rise to three Goldstone bosons, G 0 and G ± , and five physical scalar fields, two CP-even scalar fields, h and H, where by convention m h < m H , one CP-odd field, A, and one charged Higgs pair, H ± .The mixing angles α and β diagonalize the CP-even and CP-odd/charged Higgs mass matrices, respectively.
The occurrence of tree-level flavor-changing neutral currents (FCNC) is avoided by imposing a Z 2 symmetry, only softly broken by the m 2  12 term in the Lagrangian.The extension of the Z 2 symmetry to the Yukawa sector prohibits tree-level FCNCs.This results in four variants of the 2HDM, depending on the Z 2 parities of the fermion types.In this article we focus on the Yukawa type I and II.The assignment of the fermions to the Higgs doublets are listed in Tab. 1.Here we work in the physical basis of the 2HDM, where most of the free parameters in Eq. ( 1) are expressed in terms of a set of "physical" parameters given by where we use 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 the LHC at ∼ 125 GeV.The couplings of the Higgs bosons to SM particles are modified w.r.t. the SM Higgscoupling predictions because of the mixing in the Higgs-boson sector.The couplings of the neutral Higgs bosons to fermions and to gauge bosons are given by, Here m f , m W and m Z are the fermion, W boson and Z boson masses, respectively.The modification factors in the couplings to fermions and gauge bosons, ξ f h,H,A and ξ V h,H,A , for the 2HDM of type I and II are given in Tab. 2. The generic triple Higgs coupling λ hh i h j involving at least one SM-like Higgs boson h is defined such that the Feynman rules are given by h h i h j < l a t e x i t s h a 1 _ b a s e 6 4 = " a z g o f

type I type II
where n is the number of identical particles in the vertex.Relevant for our analysis here are λ hhh and λ hhH .We adopt this convention in Eq. ( 5) so that the light Higgs triple coupling λ hhh has the same normalisation as λ SM in the SM, which is given by −6ivλ SM with λ SM = m 2 h /2v 0.13.We furthermore define κ λ ≡ λ hhh /λ SM .The explicit expressions of the two triple Higgs couplings are given by where m2 , derived from m 2 12 , is given by: The triple Higgs couplings depend on c β−α .In particular, in the "alignment limit", c β−α → 0, where the light CP-even Higgs couplings to the SM particles recover SM values, and the triple Higgs couplings approach the values λ hhh = λ SM and λ hhH = 0, respectively.

Theoretical and experimental constraints
In this subsection we briefly summarize the various theoretical and experimental constraints considered in our analysis (more details can be found in Refs.[14,15]).Note, that we did not check for constraints arising from di-Higgs measurements at the LHC.The analysis performed in [6] showed that non-resonant and resonant di-Higgs searches start to cut in the parameter spaces of extended Higgs sector models.However, the parameter spaces of the CP-conserving 2HDM investigated here are not affected yet.

• Theoretical constraints
The important theoretical constraints come from tree-level perturbartive unitarity and stability of the electroweak vacuum.They are ensured by an explicit test of the underlying Lagrangian parameters (details of our approach can be found in Ref. [14]).The parameter space allowed by these two constraints can be enlarged, if we allow for a mass term breaking the imposed Z 2 symmetry softly, i.e. we choose a non-zero m 2 12 .In some of the sample scenarios that we will investigate later, we chose m 2  12 as • Constraints from electroweak precision data For SM extensions based solely on extensions of the Higgs sector the constraints from the electroweak precision observables (EWPO) can be expressed in terms of the oblique parameters S, T and U [31,32].Most constraining in the 2HDM is the T parameter, requiring either m H ± ≈ m A or m H ± ≈ m H .In Ref. [14] three scenarios were defined to meet this constraint: Here it should be kept in mind that the EWPO used to set these constraints do not take into account the recent measurement of m W at CDF [33], which deviates from the SM prediction by ∼ 7 σ.After this result for m W was published, many articles appeared to describe the CDF value in BSM models, including analyses in the 2HDM, see Refs.[34][35][36] for the first papers.It was shown that the large m W value can be accomodated by introducing a certain amount of splitting between the masses of the heavy 2HDM Higgs bosons.This also holds (albeit with a smaller amount of splitting) if a possible new world average, see Ref. [37], is taken into account [38].However, we will not include this possibility into our analysis.

• Constraints from direct Higgs-boson searches at colliders
The exclusion limits at the 95% confidence level of all relevant BSM Higgs boson searches (including Run 2 data from the LHC) are included in the public code HiggsBounds v.5.9 [39][40][41][42][43]. 2 For a parameter point in a particular model, HiggsBounds determines on the basis of expected limits which is the most sensitive channel to test each BSM Higgs boson.Then, based on this most sensitive channel, HiggsBounds determines whether the point is allowed or not at the 95% CL.As input HiggsBounds requires some specific predictions from the model, like branching ratios or Higgs couplings, that were computed with the code 2HDMC-1.8.0 [45]. 3 Constraints from the properties of the ∼ 125 GeV Higgs boson Any model beyond the SM has to accommodate a Higgs boson with mass and signal strengths as they were measured at the LHC.In the parameter points used the compatibility of the CP-even scalar h with a mass of 125.09GeV, h 125 , with the measurements of signal strengths at the LHC is tested with the code HiggsSignals v.2.6 [49][50][51].
The code provides a statistical χ 2 analysis of the h 125 predictions of a certain model compared to the measurement of Higgs-boson signal rates and masses from the LHC.As for the BSM Higgs boson searches, the predictions of the 2HDM have been obtained with 2HDMC [45].For a 2HDM parameter point to be allowed it was required [15] that the corresponding χ 2 is within 2 σ (∆χ 2 = 6.18) from the SM fit: χ 2 SM = 84.73 with 107 observables. 4

• Constraints from flavor physics
Constraints from flavor physics can be significant in the 2HDM, in particular because of the presence of the charged Higgs boson.Flavor observables like rare B decays, mixing parameters of B mesons, and LEP constraints on Z decay partial widths etc. are sensitive to charged Higgs boson contributions [52,53].To test the parameter space, taking into account the most constraining decays B → X s γ and B s → µ + µ − , the code SuperIso [54,55] was used, again with the model input given by 2HDMC (see Ref. [14] for more details).

Benchmark planes
Based on the analysis in Ref. [14] we define four benchmark planes that exhibit an interesting phenomenology w.r.t. the di-Higgs production cross sections in the gluon fusion channel, gg → hh.The four planes are all in the 2HDM Yukawa type I, as this type allows for larger deviations from the alignment limit taking into account the experimental constraints.The parameters are chosen as: 12 fixed via Eq.( 9), free parameters: c β−α , tan β

Cross section results
In this section we start our numerical analysis with the evaluation of the di-Higgs production cross sections in the benchmark planes defined in Sect.2.3.We first discuss details of the calculation and then present the results, where we analyze the impact of a possible heavy Higgs, H, in the s-channel.We perform this analysis in all benchmark planes listed above to give a broad overview about the possible phenomenology of di-Higgs production in the 2HDM.In the following sections we will discuss selected planes and points to further examine the effects of the various triple Higgs couplings and the properties of the heavy CP-even Higgs boson.

Calculation of gg → hh
The main di-Higgs production process at the LHC is given by gluon fusion.The diagrams contributing at leading order are shown in Fig. 1.They both involve a heavy quark loop (top or bottom), where for small tan β the bottom-quark loop only plays a minor role.In the SM the triangle diagram, Fig. 1a, gives access to the trilinear Higgs coupling, λ hhh , with the SM Higgs exchange in the s-channel.The box diagram, Fig. 1b, interferes destructively with the triangle diagram, resulting in a small cross section.In BSM theories additional diagrams can contribute.In particular in the case of the 2HDM the second CP-even Higgs can be exchanged in the s-channel, involving λ hhH .This diagram will usually be referred to as the "resonant diagram"5 , whereas the SM-like contributions will be referred to as "the continuum".Note that the Yukawa coupling and the trilinear Higgs self-coupling λ hhh of the SM-like Higgs boson h can deviate from the SM values so that the observed destructive interference between the triangle and box diagram in the SM may not be effective.
For our numerical evaluation we use the code HPAIR [6,29,30,56], adapted to the 2HDM.The original code evaluates the cross section of the production of two neutral Higgs bosons through gluon fusion at the LHC for the SM and the MSSM.The calculation is done at leading order (LO, see Fig. 1), and includes next-to-leading order (NLO) QCD corrections in the heavy top-mass limit.In this limit, it is assumed that the contribution of the bottom quark is negligible (it would introduce modifications of less than 1% in the SM) and then the  top mass is taken to infinity.This assumption becomes less accurate at high values of tan β because the bottom quark loop contribution gets larger.
At LO the calculation includes top-and bottom-quark loops with full mass dependence.It is equivalent to the calculation done in the Minimal Supersymmetric Standard Model (MSSM) since its Higgs sector is equivalent to the 2HDM Type II.The only changes that are implemented in the case of the 2HDM are the modification of the Yukawa couplings of the MSSM according to Tab. 2, for the corresponding type and the change of the triple Higgs couplings. 6As the QCD corrections in the heavy top-quark limit only involve couplings between coloured particles, they can straightforwardly be taken over from the MSSM to the 2HDM.For further details, we refer to Refs.[6,29,30,56].

Analysis of the cross section predictions
In Figs. 2 -5 we show the results for the di-Higgs production cross section in the 2HDM normalized to the SM value calculated at NLO for the gluon fusion process.The SM prediction was obtained with HPAIR assuming the alignment limit and coincides with the values given in Ref.
The results for all the benchmark planes are shown as follows.In the upper left plot of each figure, we present the NLO 2HDM cross sections normalized to the NLO SM value, as indicated by the color coding.The red line shows where the ratio is one, and the inner part of the solid black line is allowed by all theoretical and experimental constraints (as evaluated in Refs.[14,15]), see Sect.2.2.The upper right plot shows the K-factor, K = σ NLO /σ LO of the 2HDM hh cross sections.It should be noted that for the determination of the K-factor we consistently evaluated the LO cross section with LO pdfs and the strong coupling α s at LO and the NLO cross section is evaluatied with NLO pdfs and NLO α s Refs.[70,71].The lower left plot indicates the total width of the heavy CP-even Higgs boson, which contributes via the s-channel diagram (Fig. 1a with h i = H).This quantity will be relevant for the discussion of the dependence of the cross section on the underlying parameter space.Finally, the lower right plot shows the ratio of the full cross section, devided by the cross section obtained omitting the diagram with H in the s-channel, both evaluated at leading order.Large deviations from unity indicate an important contribution from the resonant H diagram.The analysis for the benchmark plane 1 is presented in Fig. 2, where tan β is plotted against cos(β − α) for the four quantities described in the previous paragraph.We observe that the maximum deviation from the SM prediction within the allowed area occurs precisely at the "tip" that is furthest away from the alignment limit cos(β − α) = 0 (see upper left plot).Here the enhancement factor in the cross section is ∼ 3.This corresponds to the minimum size of the triple Higgs coupling λ hhh that was obtained in the allowed region of this benchmark plane (κ λ ∼ −0.4).If a deviation between the SM prediction and experiment is observed, this could point to a deviation of the κ λ coupling.The K-factor in this region results in a factor close to 2, more specifically in the allowed region of ∼ 1.92 to ∼ 1.97 (upper right plot).One can see that the decay width of H in this region can amount to ∼ 25 GeV (lower left plot).A large total width Γ tot H has a suppressing effect on the H contribution to the cross section due to its appearance in the denominator of the s-channel propagator, as will be discussed below.We see that the H contribution has no enhancing effect on hh production within the allowed area given by the black lines (lower right plot).Including the resonance either slightly suppresses the cross section or leaves it unchanged (we find the ratio 1 exactly at the "tip", where the hh cross section is maximal within the allowed region).In conclusion, in this plane the maximum enhancement of the cross section is precisely due to the deviation of the triple Higgs coupling from the expected SM value.In the case of the benchmark plane 2, shown in Fig. 3, where we now plot m 2 12 versus cos(β − α), it can be observed that the cross section does not have a significant enhancement in the allowed region, as it does not even reach a factor of ∼ 2 times the SM value.Also for this benchmark plane we observe that the largest value of the cross section falls in the region of the minimum value of κ λ , i.e.where the destructive interference between box diagram and h exchange is minimal.The K factor in the allowed region is roughly ∼ 1.91 to ∼ 1.95, again close to 2. In the evaluation of the effect of the heavy Higgs H, we observe that neither the decay widths nor the resonant enhancement are significant in the allowed region.The H contribution almost has no effect on the hh production cross section.In the benchmark plane 3, shown in Fig. 4, where we plot m H = m A = m H ± against cos(β − α), we allow for a variation of the heavier Higgs masses by fixing the value of tan β and the mass m 2 12 according to the Eq. ( 9).Concerning the normalized NLO cross section there is no enhancement below 250 GeV since the heavy Higgs is not produced on-shell.Above this threshold there is resonant enhancement due to the contribution of the heavy Higgs in the s-channel.In particular, the cross section is enhanced by up to a factor of ∼ 8 in the "tail" at c β−α ∼ −0.1.In this region one finds κ λ close to 1, so that the enhancement of the cross section w.r.t. the SM is given by the diagram with H in the s-channel that resonantly decays into hh.This can also be seen in the lower right plot of Fig. 4, where a ratio of resonant over non-resonant cross section of up to ∼ 8 is found.The K-factor in this region is slightly above 2 (up to ∼ 2.06).
Since in this plane we allow for a variation of the mass of the heavy Higgs in the propagator one can observe the enhancement of the cross section around m H ∼ 350 − 400 GeV that is expected from single Higgs production above the di-top threshold, m H ∼ 2m t .This enhancement is clearly visible in the lower right plot of Fig. 4 (particularly for negative c β−α within the allowed region).The final plane we investigate in detail is plane 4 as shown in Fig. 5.Here the two free parameters, that we plot against each other, are m H and m A = m H ± .Correspondingly, one expects a large enhancement of the 2HDM cross section for m H ∼ 400 GeV where m H > 2m h and above the di-top threshold.Exactly this can be observed in the upper left plot (total cross section) and the lower right plot (relative enhancement from the resonant diagram).The cross section can be up to 60% larger than the SM di-Higgs production cross section, with a K factor again close to 2. The total width of the heavy CP-even Higgs ranges from very small values up to larger than 200 GeV within the allowed region (found for low m A = m H ± and large m H ). For large total widths the resonance enhancement is not effective.The largest enhancements of the cross section are found for relatively small values of Γ tot H < ∼ 10 GeV.

Dependence on triple Higgs couplings
In this subsection we analyze the cross section with respect to the triple Higgs couplings involved in the di-Higgs production process.In particular, we will show in which part of the parameter space the total di-Higgs production cross section has a relevant dependence on λ hhh and/or λ hhH .
Here we focus on a statistical treatment of the errors of the total cross section measurement, which is assumed to be Gaussian, neglecting systematic uncertainties.It was found in Refs.[8,72] that the statistical uncertainty of the total di-Higgs cross section measurement, assuming SM values, will reach a level of 4.5 σ at the end of the HL-LHC, combining ATLAS and CMS.(Taking into account systematic effects could lower this value to ∼ 4 σ.)Consequently, we will approximate the corresponding error in the measurement as δxs = xs/4.5. 8The significance of the deviation of the (to be measured) 2HDM cross section w.r.t. the SM value can then be expressed as It should be noted that this approximation becomes worse for larger deviations of xs 2HDM from xs SM , since the precision of the measurments, δxs, has been evaluated assuming the SM value.For higher cross sections a more precise measurement can be expected.We present our results within the four benchmark planes discussed above in Figs. 6 -10.For each plane in the upper left and middle plot we will show the predictions of κ λ and λ hhH , as obtained in Refs.[14,15].The upper right plot, for better comparison, repeats the results of σ 2HDM /σ SM at NLO QCD as presented in the upper left plots in Figs. 2 -5, where we here also show the maximum (minimum) value of the coupling that is realized within the allowed region as red (blue) dots.The three upper plots are always given in the plane of the two free parameters involved in the respective benchmark choice.The lower left plot shows which combinations of κ λ and λ hhH can be reached in each plane, where the points inside the area allowed by theoretical and experimental constraints are marked in red (and indicated with a red arrow).The lower middle plot, focusing on the allowed region in the κ λ -λ hhH plane, presents the values of σ 2HDM /σ SM at NLO QCD in this plane (with the SM point κ λ = 1, λ hhH = 0 marked by a red star).This indicates the dependence of the total 2HDM di-Higgs production cross section on the two involved triple Higgs couplings.The black points represent the values of the THCs that are reached in this plane.The lower right plot shows the same area in the κ λ -λ hhH plane, now indicating the expected number of ∆σ SM , see Eq. (11), that the (to be measured) 2HDM result differs from the SM prediction, i.e. with which significance such a deviation can be measured experimentally.
In benchmark plane 1, Fig. 6, one can observe from the comparison of the upper left and right plots that within the allowed range, as discussed above, the smallest κ λ value gives rise to the largest value of σ 2HDM .As can be inferred from the lower middle plot, in this benchmark plane the cross section depends strongly on κ λ , but effectively not on λ hhH .This is due to the fact, as discussed above, that the heavy CP-even Higgs is too heavy to give a sizable s-channel contribution.Overall, one can observe that for the smallest allowed κ λ values, κ λ ∼ −0.4, a cross section enhancement of up to ∼ 3 can be found.
Finally, we see from the lower right plot that for the smallest κ λ , corresponding to the largest σ 2HDM , a deviation of up to ∆σ SM ∼ 9 can be expected.This indicates that within this 2HDM benchmark plane a clear distinction between the 2HDM and the SM via the di-Higgs production cross section can be possible.Deviations of more than 2σ can be expected for κ λ < ∼ 0.6.In benchmark plane 2 a similar result as in plane 1 can be observed, as shown in Fig. 7.The largest cross sections are found for the smallest κ λ values, and the predicted 2HDM di-Higgs cross section depends only mildly on λ hhH .The latter can again be understood because of the relatively large value of m H = 650 GeV in this benchmark plane.The maximum significance of the 2HDM deviation w.r.t. the SM value is less than for plane 1 with a value of at most 3.5 σ, reached for κ λ ∼ 0.9 and λ hhH ∼ −0.5.
The situation is more involved in plane 3, which we show in Fig. 8.As discussed in the previous subsection, very large enhancements can be reached in this parameter plane, and larger allowed regions are found for both signs of c β−α .The lower middle plot, showing the cross section enhancement w.r.t. the SM seems to show a relatively small enhancement of up to ∼ 3. The larger effects that actually occur (enhancements of up to ∼ 8) are found in comparably small regions and are thus not well visible in this figure (but will be shown clearly below).The lower right plot shows the ∆σ SM , and for large parts of the parameter space we find the same feature as in the two benchmark planes above: the largest deviations are found for the smallest κ λ values, and independent of λ hhH , reaching about 5 σ.However, In Fig. 9 we present the results where we have divided the allowed region of benchmark plane 3 and the obtained cross sections and their sensitivity in the couplings plane in positive and negative values of c β−α .The left (right) column in Fig. 9 shows the results for c β−α negative (positive).The upper row indicates the combination of κ λ and λ hhH that can be reached in the full benchmark plane.The middle row shows σ 2HDM /σ SM , whereas the lower row indicates the level of ∆σ SM that can be reached, where we have zoomed further into the interesting region.The middle left plot demonstrates that for κ λ close to 1 the cross section can be strongly enhanced with the enhancement strongly depending on the BSM THC λ hhH .This behavior can be traced back to the H contribution in the s-channel for relatively small values of m H , as discussed in the previous subsections.Looking into the (zoomed in) result for ∆σ SM (lower left plot) one can observe that for the smallest values of λ hhH ∼ −0.2 a deviation from the SM by up to 35 σ can be expected.More importantly, the size of the deviation may give an indication of the value of λ hhH .Going to positive values of c β−α as presented in the right column, one can observe that for large parts of the paramerter space a dependence solely on κ λ is found, as in the previous benchmark planes.However, again for κ λ ∼ 1 strong enhancements are found due to the presence of the heavy CP-even Higgs in the resonance.This is better visible in the lower right plot, demonstrating that for small, but positive values of λ hhH deviations of up to 6 σ can be seen, whereas for negative values of λ hhH ∼ −0.2 even deviations of up to 9 σ can be found.Here it should be kept in mind that this analysis only demonstrates the possible dependences and effects of the THCs on the di-Higgs production cross section.As will be discussed below, an actual possibility for a determination of λ hhH is not implied, as it depends on the precise knowledge of the other (free) parameters.
Our final analysis in this section is done for benchmark plane 4, as shown in Fig. 10.In this m A = m H ± -m H plane the value of κ λ is always close to 1, varying only by about ∼ 10%, so that large variations of the di-Higgs production cross section can only be produced by resonant enhancement.The coupling λ hhH varies between 0 to ∼ −1.5 in the allowed region (and even down to ∼ −4.5 for the largest m H values).The cross section, as discussed already in the previous subsection shows an interesting enhancement of up to ∼ 60%, where the heavy Higgs is resonant and not too heavy, m H ∼ 400 GeV.We will use the behavior of σ 2HDM in this benchmark plane for a more detailed analysis of the invariant m hh distribution in the next section.
The projection into the κ λ -λ hhH plane shows only a line, which can be understood as follows.Looking at Eq. ( 6) and Eq. ( 7) and discarding all the terms proportional to constants (the angles α and β) we find the following relations: where the c 1,2,3,4 are constant terms, and it is taken into account that according to Eq. ( 9) m2 ∝ m 2 12 ∝ m 2 H . Consequently, one finds resulting in the linear dependence that is observed in the lower plots of Fig. 10.Now we proceed to analyze the values of the cross section that are possible for the different values of the THCs, even though we do not have a truly 2-dimensional plot in these cases.
For the benchmark plane 4, in the lower middle plot of Fig. 10 the cross section is badly defined in the sense that more than one value of the cross section corresponds to a particular value of the THCs.This happens when we allow for a change in the masses but fix the angles, as discussed above.The THCs change in a coherent way for different masses m H (see upper left and middle plots in Fig. 10), while the cross section has different possible values.As an example, for m H in the range ∼ 220 GeV to ∼ 800 GeV it can vary within ∼ (0.8 − 1.6) × σ SM , as can be seen in the upper right plot.Therefore, in the lower middle plot we represented the mean value of the cross section as a circle for a particular combination of (κ λ , λ hhH ).We show maximum (upper triangle "slightly displaced above the circles") and minimum (lower triangle "slightly displaced below the circles") values of the normalized cross section for this same combination of (κ λ , λ hhH ).One can observe that the highest cross section is realized for a κ λ ∼ 0.985 and λ hhH ∼ −0.2.In this point the value of the cross section varies roughly from 1 to 1.6 and the sensitivity that can be reached in the most optimistic scenario (i.e. the largest deviation from the SM that is realized) is almost 2.5 σ, as can be seen in the lower right plot.This enhancement is relatively small, but it demonstrates that in this plane the relevant role of the coupling λ hhH is more significant than κ λ , which is very close to 1, where these effects are found.The size of the deviation clearly depends on λ hhH in this case.

Analysis of m hh
In the next step of our analysis we will analyze the influence of the THCs on the di-Higgs production cross sections by evaluating the di-Higgs invariant mass distribution, m hh .We first demonstrate in a toy example the effect of the characteristic properties of the resonant Higgs boson: its mass, its width, and the sign of the coupling combination (λ hhH × ξ t H ). Subsequently, the analysis will be performed for several benchmark points located in the planes discussed in the previous section, where the effects of the characteristics of the resonance will be demonstrated in a real model.
The invariant mass distributions, dσ/dm hh , are also obtained with the code HPAIR.We will use a grid of values for the invariant mass m hh that range from 250 GeV to 1250 GeV.As a default value we will use a bin size of 20 GeV (where experimentally a bin size of ∼ 50 GeV appears more realistic, see the discussion below).This bin size is used for demonstrative purposes.In a later step we will analyze the effect of different bin sizes and other experimental effects to obtain a more realistic picture of m hh distributions.

General analysis of the effects
In this subsection we will analyze a toy model for the resonance to demonstrate the effects of the mass, width and couplings of the resonant Higgs boson in the s-channel exchange.
The effect of the total decay width of the heavy Higgs boson, Γ tot H , is important whenever the resonant diagram gains significance in the calculation of the cross section.This happens close to the resonance at m H ∼ m hh , as discussed in the previous section.For a correct treatment close to the resonance the total width has to be included into the propagator, From this expression one can clearly see that the dominant effect of Γ tot H appears when the intermediate Higgs boson mass is equal to the (reduced) center of mass energy Q 2 .In the Higgs pair production process, the total decay width of the heavy Higgs becomes relevant near the resonant region where the behavior of the cross section can be dominated by the interference between the resonant and the non-resonant contributions, which is proportional to We use this expression to investigate the resonant behavior of the m hh distribution.In Fig. 11 (left) we show σ interf as a function of m hh = Q.We have chosen m H = 300 GeV and three exemplary values of Γ tot H : 0.1 GeV (red), 10 GeV (dark blue) and 50 GeV (light blue).In all cases σ interf shows a peak-dip structure, with the change exactly at m hh = m H , as expected from Eq. ( 15).Furthermore, one observes that the highest (smallest) peak-dip structure is obtained for the smallest (highest) value of Γ tot H , following the analytical behavior of Eq. ( 15).We furthermore observe that the "total width of the effect", given by the width of the peak at half of its maximum value, increases with increasing Γ tot H , as expected.The features observed for σ interf are also found in the calculation of the m hh distribution of the complete cross section, i.e. the result from taking into account the complete resonant and the non-resonant contributions, as shown in the right plot of Fig. 11.Here we depict dσ/dm hh as a function of m hh for one benchmark point of benchmark plane 4 with m A = m H ± = 544.72GeV and m H = 515.5 GeV9 .For the total width of H we find Γ tot H ∼ 3 GeV, resulting in the red curve.In order to illustrate the effects of the size of Γ tot H , as seen in the left plot, we also show the results for ad-hoc set values of Γ tot H = 10 GeV (dark blue) and 50 GeV (light blue).The main features of Γ tot H (height of the peak-dip structure and the "width of the effect") are found in the full calculation exactly as in σ interf .However, there is one important difference between the results for σ interf and the full invariant mass distribution results which can be observed in Fig. 11.While for σ interf a peak-dip structure is found, in the full calculation for our chosen benchmark point a dip-peak structure can be observed.This difference can be traced back to the sign of (λ hhH × ξ t H ), which enters as prefactor in the the resonant diagram.In the left plot of Fig. 12 the resonant H diagram is shown with the two coupling factors entering the amplitude: the top-Yukawa coupling modification factor ξ t H of the heavy Higgs H and the THC λ hhH .The right plot in Fig. 12 demonstrates the effect of sign(λ hhH × ξ t H ). The red curve is identical to the red curve in the right plot of Fig. 11, as obtained for the value of λ hhH = −0.3975,with λ hhH × ξ t H < 0. The blue curve shows the m hh distribution for the (ad hoc) flipped sign, i.e. with λ hhH = +0.3975(normalized to the corresponding value of the full cross section obtained for this trilinear Higgs coupling).As can be expected, the flip of the sign(λ hhH × ξ t H ) also flips the dip-peak structure to a peak-dip structure (σ interf shown above corresponds to a positive sign).The question arises whether an experimental analysis can be sensitive to the difference between peak-dip and dip-peak, and thus provide a handle on the sign of (λ hhH × ξ t H ).This question will be analyzed below.

Model based analysis: benchmark plane 3
We start our model based m hh analysis for several points given in the benchmark plane 3. First we will investigate the points that present the largest enhancement of the cross section w.r.t. the SM within the allowed region, as listed in Tab. 3. Second we will look at points with c β−α ∼ 0.1, as listed in Tab. 5. Finally we will analyze points with c β−α ∼ 0.2, i.e. a relatively large deviation from the alignment limit, as listed in Tab. 6.The aim of this study is to extract the general behavior and the influence of specific parameters on the experimental measurement of the cross section.This will allow us to track variations of the parameters that we are mostly interested in (m H , Γ tot H and λ hhH ).

Benchmark plane 3: large di-Higgs production cross sections
We first analyze three points in benchmark plane 3 with large enhancements of the di-Higgs production cross section w.r.t. the SM.They are located close to the alignment limit and can be seen in the left part of Fig. 13 as red, blue and black dots.In the first step of the analyses we choose a bin size of 5 GeV to make the large resonant enhancement, which is very narrow, clearly visible.The values of the parameters of each point are listed in Tab. 3. The di-Higgs production process is kinematically forbidden for m hh < 250 GeV = 2 m h .Once this threshold is surpassed, one can observe a resonant enhancement for m hh ∼ m H .This is clearly seen at the location of the resonant peaks in the invariant mass distribution in Fig. 13 (right).For the red and the black points the resonant peak is found around ∼ 265 GeV, while for the blue point it is located at ∼ 292 GeV, corresponding to the respective m H value.The "height" and "width" of the peaks is related to the total decay width Γ tot H of the heavy Higgs boson, which is largest for the red point (σ/σ SM ∼ 8) and smallest for the black point (σ/σ SM ∼ 2.5), as shown in Tab. 3. Furthermore, the resonant heavy Higgs contribution yields the already observed typical pattern, a peak-dip or dip-peak structure, depending on the parameter point.The peak-dip structure is observed in the blue and red points, whereas in the case of the black point one can only see the peak, see the discussion below.Moreover, we observe for all three points an enhancement at an invariant mass of 350 − 400 GeV, which is related to the top pair production threshold in the resonant diagram 10 .
The three points have different values of c β−α , which change the Yukawa coupling of the top quark according to the expression: resulting in a positive (negative) sign of (λ hhH × ξ t H ) for the red and blue (black) points, as shown in Tab. 2. This in turn results in a peak-dip (dip-peak) structure, cf.Tab. 4. The pattern of the differential distribution changes according to the sign of (λ hhH × ξ t H ). For the black curve only a peak is visible, because the dip before it cannot be produced for masses below 250 GeV.
A more detailed analysis of the three points is presented in Fig. 14, where we analyze the contribution of individual (groups of) diagrams.In the upper left plot we show the total, i.e. including all diagrams, distribution for the three points (as in the right plot of Fig. 13)   but changing the bin size to 20 GeV in order to represent a more realistic experimental set up.One can observe already in the example of the black point that the bin size (and location) is important for the observation or non-observation of an enhancement.
In the upper right, and lower row of Fig. 14 we disentangle the contributions of the individual diagrams for the red, black and blue points, respectively.We have calculated the total differential distribution including all three diagrams for di-Higgs production (red), the SM-like cross section, called continuum, -including only the box and the SM-like Higgs boson in the s-channel diagrams -(black), the contribution of the diagram with no triple Higgs couplings involved, i.e. the box diagram (yellow), the contribution of the two diagrams with the triple Higgs couplings, which includes the h and H in the s-channel (purple), and the contribution with only h (light blue) and with only H (dark blue), respectively, in the s−channel.For all three points the same pattern can be observed.The 2HDM cross section (red) follows largely the SM-like distribution (black), apart from the strong resonant enhancement at m hh ∼ m H .This is caused by the H s-channel contribution (blue), potentially providing sensitivity to λ hhH .Furthermore, the destructive interference of the box-diagram (yellow) and the SM-like h-exchange contribution (light blue) is clearly visible in the continuum, i.e.SM-like, contribution of the distribution (black).

Plane 3:
Next we proceed to study a set of points that are all located at the same value of c β−α ∼ 0.1.The exact value of c β−α = 0.1015 results from the grid used to scan this plane.In this case the only change between the benchmark points, as listed in Tab. 5, is the common mass m H of the heavy Higgses, with correspondingly modified couplings and decay widths.The points are also shown as orange, yellow, purple, garnet and green dots (in ascending mass order) in the upper left plot of Fig. 15 (repeating the upper left plot of Fig. 8).The upper right plot of Fig. 15 (repeating the middle right plot of Fig. 9) indicates the location of these points in the κ λ -λ hhH plane.One can observe that with increasing m H the points are decreasing in κ λ (from κ λ ∼ 1 down to zero) and are increasing in λ hhH (from λ hhH close to zero to λ hhH ∼ 1).
The m hh distributions for the five benchmark points are presented in the lower plot or Fig. 15, with a bin size of 5 GeV.The color indicates the value of m H (as defined in the upper plots).The red arrows indicate the location of the "resonant peaks" (for m H > 250 GeV).The resonant enhancement is found to be tiny, despite the non-negligible values of λ hhH .The reason for the unrealistically small bin size of 5 GeV is to see any peak at all.It is clear that for these points the enhancement of σ 2HDM w.r.t.σ SM is caused purely by the reduced κ λ values which alleviates the destructive interference between triangle and box contributions present in the SM.The resonant H exchange hardly contributes to the total cross section.The structure of the enhancement in this case is hard to infer from the plot, but looking closely one can see a peak-dip structure.The reason for this small resonant contribution can be found in the top Yukawa coupling value of the heavy CP-even Higgs.Following Eq. ( 16) we obtain a value of ξ t H = 5 × 10 −4 > 0, thus rendering the resonant contribution negligible.The triple Higgs couplings λ hhH listed in Tab. 5 are all positive, so that the overall sign for the coupling factor of the triangle contribution is positive and the (hardly visible) structure is a peak-dip one as expected.Furthermore, we have seen in Fig. 14 that the largest contribution in the lower mass spectrum comes from the diagram with a light Higgs h exchange, i.e. from the diagram involving κ λ .This diagram hence drives the behavior of the distribution.In the lower plot of Fig. 15 this trend is clearly visible in the lower part of the spectrum, m hh < ∼ 350 GeV.The smaller the value of κ λ for a particular point (as seen in upper right plot), the larger the enhancement in the invariant mass distribution at lower m hh .The most extreme point is the green one, for which κ λ is close to zero, and the m hh spectrum shows a clear bump at m hh ∼ 350 GeV.It should be recalled that the gluon fusion di-Higgs production cross section has a minimum for κ λ ∼ 2.5 and is higher for small (or very large values) of κ λ .

Plane 3:
We finish our analysis of benchmark points in plane 3 with five points with a relatively large value of c β−α ∼ 0.2, i.e. relatively far away from the alignment limit, as given in Tab. 6.As above, the exact value of c β−α = 0.203 is given by the scanned grid.The points are shown in the upper left and upper right plot of Fig. 16 as colored dots in orange, yellow, purple, garnet and green dots (in ascending mass order), color coding indicating σ/σ SM in the c β−α -m H and κ λ -λ hhH plane in the upper left and right plot, respectively.As can be observed in the upper right plot, all points have κ λ ∼ 1, i.e. no relevant change in the cross section can be expected from the contribution of the h-exchange diagram, so the lower part of the m hh spectrum is very similar for all the points.On the other hand, the values of λ hhH decrease from around zero down to λ hhH ∼ −0.5.However, as can be seen in the two upper plots of Fig. 16, the variation of the total cross section is relatively small.The largest cross sections are found for m H = 312.0GeV (yellow) and m H = 399.75GeV (purple), i.e. around ∼ 350 GeV, see the discussion in Sect.4.2.1.In the lower plot of Fig. 16 we show the m hh distribution for the five points.For four of them with m H > 2m h a clear resonance dip-peak structure can be observed at m H ∼ m hh , as expected.In Fig. 17 we present a more detailed analysis of the sensitivity to the triple Higgs couplings in the invariant mass distribution following the same notation as in Fig. 14.We show the invariant mass distribution from all the diagrams in the upper right plot and then split the individual contributions for each particular mass point in the rest of the plots.The first one for m H = 244.5 GeV is shown in the upper middle plot, which has m H below the di-Higgs production threshold.Consequently, no enhancement due to the diagram containing λ hhH can be observed, and the total cross section is almost indistinguishable from the SM-like result in this case.For the other masses we find a similar result as in Fig. 14.One can observe that the s-channel contribution involving the heavy Higgs with its trilinear coupling λ hhH is responsible for the enhancement close to the mass of the intermediate Higgs boson, while the effect of the κ λ is mostly significant in the low mass region of the plot.The contribution of the diagrams involving THCs (purple) interferes with the continuum (box diagram) shown in yellow and creates the dip-peak structure that can be observed in the total distributions, see Fig. 17 (upper left).
The top Yukawa coupling for this value of c β−α is ξ t H = 0.1020 > 0, and thus the sign of (λ hhH × ξ t H ) is negative, resulting in the dip-peak structure observed.In principle, this type of distributions can yield a handle on the size and sign of λ hhH , as will be discussed in more detail below.
In Fig. 18 we further analyze the interference contributions of the s-channel diagrams in Fig. 1.In the left plot we show the interference term of the diagram with the s-channel h exchange (A) and diagram with the s-channel H exchange (B).The interference term in this case is defined as |A + B| 2 − |A| 2 − |B| 2 .The solid (dashed) lines indicate positive (negative)  interference.One can observe a similar behavior to the above discussed interference between resonant and box contributions, i.e. that these two diagrams interfere constructively up to m hh ≤ m H , and destructively for larger m hh values, i.e. the interference term changes its sign.The diagram A and the box diagram (C) interfere negatively accross the whole invariant mass range as shown in the middle plot of Fig. 18.This behavior corresponds to the result found for the SM di-Higgs production, where only these two diagrams are present Finally, the interference of B and C, shown in the right plot of Fig. 18, has two sign changes.Up to m hh ≤ m H the interference is negative, leading to the dip in the total m hh distribution.For larger values it turns positive, leading the subsequent peak in the total distribution, see the discussion in Sect.4.1.The second sign change happens because the interference approaches zero at an m hh value not related to m H .In the plot the interference lines in principle go down to zero, which, however, is not visible due to the log scale and the finite bin width.

Model based analysis: benchmark plane 4
To complete our m hh analysis we investigate benchmark points in plane 4. In this plane the heavy Higgs boson masses are free parameters and are in agreement with the applied constraints over a relatively large interval, i.e. the effects of a mass variation (as well as a variation of λ hhH ) can be readily analyzed.The selected points are listed in Tab. 7.They are all located at the same mass of m A = m H ± = 544.75GeV.As was demonstated in Sect.3.3, κ λ and λ hhH are proportional to each other, and allowed values of λ hhh are within ∼ 2% of κ λ = 1, whereas λ hhH is found for our benchmark points in the interval [∼ 0, −0 The location of the points in the benchmark plane is shown in the upper left plot of Fig. 19, and in the κ λ -λ hhH plane in the upper right plot.The lower row shows the corresponding m hh distributions with the color indicating the m H value.In the left plot, for better visibility, we show the idealized case of a bin size of 5 GeV, whereas the right plot shows the more realistic case of 20 GeV.
One can observe that the resonant structure in this case is dip-peak since (as in the previous case) the top Yukawa coupling is positive and the BSM THC is negative.For the purple point no enhancement can be observed, since its mass is below the di-Higgs production threshold, m purple H = 222.5 GeV.The remaining points show a resonant enhancement as the invariant mass approaches the mass of the heavy Higgs, i.e. one expects to have the highest sensitivity to the THC λ hhH in this region.The difference between the heights of the peaks is found to be largest for small m H , but naturally also decreases with increasing bin size, which will be further discussed below.The increasing bin size, in particular, decreases the visibility of the dip, an effect that is in particular visible for low m H .This demonstrates the relevance of the bin size in a realistic experimental analysis (and will also be further discussed below).

Impact of experimental uncertainties
In this section we will analyze the impact of experimental uncertainties on the possible sensitivity to λ hhH .These effects are the experimental smearing, i.e. the uncertainty in the m hh measurement, the experimental resolution, i.e. the size of the bin width, as well as the arbitrary location of the bin.We will also demonstrate how the experimental results for λ hhH change with a variation of sign(λ hhH × ξ t H ). In order to estimate the sensitivity to the BSM coupling λ hhH , following Ref.[22], we define a theoretical quantity that aims to quantify the "sensitivity to λ hhH " (but is not meant as a determination of its precision, which requires a detailed experimental analysis, which is beyond the scope of our paper).We define where N R is the number of events of the resonant contribution, and N C is the number of events of the continuum.The window in which the events are counted is defined by The sum over i in Eq. ( 17) runs over all the bins that fulfill this condition.The chosen condition in Eq. ( 18) starts with a minimum of 1000 excess events due to the resonance when the bin size is 50 GeV and 200 events when the bin size is 10 GeV, i.e. smaller bin sizes are not "punished".Using the absolute value in the definition of R in Eq. ( 17), as well as in the definition of the window in Eq. ( 18) effectively makes use of both the dip and the peak of the smeared distribution.This constitutes a simplified theory definition, where in a realistic experimental analysis the dip-peak structure would be taken into account via a template fitting, see e.g. the analysis in Ref. [73].The numbers of events are in turn obtained using the relation between the cross section and the integrated luminosity of the collider, where we have used L = 6000 fb −1 , i.e. the sum of the anticipated luminosity of ATLAS and CMS combined at the end of the HL-LHC run.This constitutes the most optimistic case.

Smearing
Differential cross section measurements are affected by the finite resolution of the detectors.This translates into a blurred or "smeared" spectrum that can be observed in such experiments.We try to mimic this effect by artificially smearing the theoretical prediction for the invariant mass distributions of the chosen benchmark points.To do this we introduce a statistical error to our prediction of the invariant mass.We apply the uncertainties in m hh by allowing the value of an event to shift to the left or to the right in the spectrum according to a Gaussian probability distribution.The amount of smearing is defined in terms of a percentage of smearing p that indicates that the predicted value x should stay within the interval [(1 − p) • x, (1 + p) • x] with a 78 % probability, which corresponds to the definition of the full width half maximum of the Gaussian distribution (as given for the experimental analyses).We illustrate this effect in Fig. 20 for one particular example of a benchmark point taken from the benchmark plane 4 with the masses fixed to m A = m H ± = 544.72 GeV and m H = 515.5 GeV11 .In this figure we show in blue the m hh distribution without smearing (the ideal case).The solid line depicts the full distribution, whereas the dashed line shows the result for the continuum (non-resonant) diagrams.The red lines demonstrate the effect of applying a 10% (left plot) and 15% (right plot) smearing on the theoretical prediction of the m hh distributions, where the solid (dashed) lines indicates the full (continuum) result.
While a 15% smearing was given as a realistic future estimate, the 10% smearing indicates a potential optimistic improvement.One can observe that from the original dip-peak structure as seen in the solid blue line effectively only a peak or bump around the original peak remains.The original dip is visible only as a very small reduction of the unsmeared distribution, as the relative weight of the points below the continuum is smaller than those above the continuum (note the logarithmic scale).

Bin width
As a further step in the evaluation of the experimental challenges, we analyze the effect of the bin width.The binning means that the data in a particular interval in m hh is presented as the mean value of the differential cross section of all the points that fall in that interval.Assuming that at least one of the Higgs bosons analyzed will decay in a b b pair 12 , the bin size will eventually be determined by the b-jet mass resolution from the reconstruction of the h → b b decay mode.This affects the visualization of the results in a realistic experimental set up, but also the counting of events for the evaluation of the experimental sensitivity, see Eq. ( 17).The binning is applied after the smearing discussed in the previous subsection.
The analysis is performed for the same benchmark point as in Sect.5.1.In Figs.21 and 22 we show the same spectrum but for a different bin size in the m hh variable: 10 GeV (upper left), 20 GeV (upper right), 40 GeV (lower left) and 50 GeV (lower right).Fig. 21 assumes a 10% smearing, whereas in Fig. 22 we show the more realistic result with 15% smearing.The red lines show the true (smeared and binned) prediction, whereas the other colors indicate the unsmeared, but binned results for comparison.One can observe that the effect of the smearing becomes less significant in the region of resonant production for a larger bin size.The resonance is already partially diluted by the smearing, and the effect of the binning becomes less visible, as can be observed best in the lower right plots of Figs.21 and 22.The effect of the binning is less important once the smearing of the experimental data is taken into account.After the binning the "dip" is effectively indistinguishable from the continuum contribution.The peak is still persistent and for larger bin size approaches the same height as the bump at ∼ 400 GeV before binning.
In the most conservative result the expected experimental resolution should have a bin size of 50 GeV and a smearing of ∼ 15%.The expected results in this case would possibly give access to the location of the resonance (the mass of the CP-even H should be know via single production by the time the di-Higgs cross section is measured) and partially to the height, and thus possibly to the size of λ hhH .In order to make a quantitative estimate of the sensitivity of the signal produced by the resonant diagram we have calculated the value of the variable R defined in Eq. ( 17) that is obtained from Figs. 21 -22, shown in Tab. 8. Overall, one can see that the values of R are significant, i.e. the signal could possibly be distinguished at the HL-LHC if our assumptions on the experimental uncertainties are met.It should be noted that we are not taking into account the efficiency of the particle detectors, which could reduce signifficantly the estimate of R. Comparing the two columns in Tab. 8 one observes R is roughly 10% worse as the assumed percentage of smearing increases by 5%.However, R is is somewhat more stable w.r.t. the bin size, where deviations within ∼ 5% are found.This would constitute a rather positive feature for an experimental set-up.

Bin location
The next part of the analysis concerns the arbitrary choice of the location of the bin.This choice can also affect the pattern of the invariant mass distribution.The concrete value of m min hh (the value of m hh at the bin start) and m max hh (the value of m hh at the bin end) affects the number of events that fall into that bin and thus can have an impact on the evaluation of the sensitivity R. For the previously used benchmark point we change the location of the bin for a 10% and 15% smeared distribution and a 50 GeV bin size.In Fig. 23 we show the difference in the invariant mass distribution created by a change in the location of the first bin by 10 or 20 GeV, i.e. we start the distribution at 251, 261, 271 GeV as orange, yellow and purple lines, respectively.In both plots we show the difference between the total differential cross section (solid lines) and the continuum contribution (dashed lines).The left (right) plot uses a smearing of 10% (15%).One can observe that for all three choices of bin locations the peak structure remains similarly visible (the dip is strongly diluted from the smearing and the binning as discussed in the previous subsections).To quantitatively evaluate the significance of the signal of the resonant enhancement we list the values of R for the two plots discussed above in Tab. 9.In Tab. 9 one can observe that the variation in R stays within 5% when we modify the location of the bins.That means that the uncertainties associated to the location of the bin are smaller than the ones associated to the smearing and about the same as for the bin size.Therefore, overall we find that the experimental resolution of the particle detector, which we tried to mimic by smearing the data, has a larger impact on the resonance, and the width and location of the binning has a smaller effect in diluting the resonance.

Effect of sign(λ
Finally, taking into account the above discussed experimental uncertainties, we analyze the possibility to access experimentally the relative sign of the involved couplings in the resonant diagram.For this analysis we will use the same benchmark point as above, where the effect of an ad-hoc change of sign(λ hhH × ξ t H ) has been shown already in the right plot of Fig. 12 and is reproduced for completeness of this analysis in Fig. 24.On this "ideal" result we will apply a smearing and a binning in order to see whether the difference in the peak-dip vs. dip-peak structure, and thus in the sign(λ hhH × ξ t H ), persists in the experimental analysis.This will demonstrate whether the access to this parameter is a realistic goal at the HL-LHC.The invariant mass distribution of the benchmark point is shown in red for the original value of λ hhH = −0.3975and a positive sign of (λ hhH × ξ t H ), whereas the blue line gives the result for an ad-hoc changed sign of λ hhH = +0.3975and a negative sign of (λ hhH × ξ t H ). Both distributions are normalized to the corresponding value of the total cross section, which is indicated in the legend of Fig. 24.The change in the structure is clearly visible in this plot.We remind that the value of the top Yukawa coefficient in this case was, as obtained from Eq. ( 16), ξ t H = 0.104 > 0. Therefore we see that for an overall minus (plus) sign the imprinted structure of the resonance in the m hh spectrum is a peak-dip (dip-peak) one.
We now impose the various experimental uncertainties on the "ideal" result shown in Fig. 24.We start with the application of smearing as presented in Fig. 25.In the left (right) plot a 10% (15%) smearing is applied.The red and blue curves correspond, as before, to the negative and positive sign of λ hhH .The black curves indicate the continuum distribution.One can observe that the blue curve is above (below) the red one right before (after) the bump.However, the overall structure cannot be resolved once the data is smeared.The final step is the application of binning on top of a 15% smearing, which is the most realistic experimental set up.The results are shown in Fig. 26.One can conclude that effectively the dip-peak vs. peak-dip structure becomes unresolvable once both smearing and binning is applied.

Conclusions
In this work we have analyzed the impact of triple Higgs couplings on the production cross section of two ∼ 125 GeV Higgs bosons at the HL-LHC in the framework of the 2HDM.The first goal was to analyze the impact of λ hhh and λ hhH on the di-Higgs production cross section.In a second step we analyzed the potential sensitivity of the HL-LHC on the BSM THC λ hhH .This sensitivity has been analyzed w.r.t.various experimental uncertainties.
We have chosen several observables that allow us to trace the impact of the THCs.The first one is the Higgs pair production cross section that is affected by the diagram containing this coupling.Using the code HPAIR we have evaluated the cross sections of di-Higgs production including NLO QCD corrections (in the heavy top limit) in the specific benchmark planes and analyzed the particular impact of the THCs, the contribution of the heavy Higgs boson in the propagator and its decay width.
Differences of the di-Higgs production cross sections can originate from a changed value of λ hhh , as well as from additional "resonant" contributions, given by the s-channel exchange of the heavy CP-even Higgs boson, H.In the analyzed benchmark planes we have found both types of effects, which can yield a strong enhancement of the di-Higgs production cross section w.r.t.its SM result.While these results indicate that the effect of the resonant s-channel H contribution may be visible, from the cross section alone it is not possible to disentangle the two sources.
In order to gain more direct access to λ hhH we analyzed the invariant di-Higgs mass distributions, dσ/dm hh as a function of m hh .In a toy example we analyzed the interference of the resonant H exchange with the non-resonant (continuum) diagrams.The interference term changes sign for Q = m hh = m H , resulting in a peak-dip (or dip-peak) structure of the m hh distribution.We demonstrated that the interference effect becomes smaller with the total decay width of the H and that the "total width of the effect", given by the width of the peak at half of its maximum value, increases with increasing Γ tot H . Furthermore, the structure of the interference (peak-dip vs. dip-peak) depends on the sign of (λ hhH × ξ t H ), which enters as a prefactor of the interference contribution (ξ t H denotes the top Yukawa coupling of the heavy CP-even Higgs boson).
The effects in the toy example were reproduced for the complete m hh distribution for several points in the before chosen benchmark planes.We demonstrated the dependences of the size of the interference effects on λ hhH and ξ t H for large di-Higgs production cross sections, as well as for selected points closer to the alignment limit (c β−α ∼ 0.1) and further away from the alignment limit (c β−α ∼ 0.2).These results indicate that the HL-LHC may have sensitivity to see effects of the BSM THC λ hhH .As a theoretical variable quantifying the relative sensitivity to λ hhH , the variable R was defined, see Eq. ( 17).However, it is not meant as a determination of the experimental λ hhH precision that requires a detailed experimental analysis, which is beyond the scope of our paper In order to further analyze a possible sensitivity, we applied various experimental uncertainties on R. The first one is the smearing due to the detector resolution, which realistically can go down to ∼ 15%, but we also analyzed a potential (optimistic) improvement down to 10%.It is found that from the dip-peak structure mostly a bump survives, with a very small reduction due to the dip w.r.t. the unsmeared result.The next effect is the binning of the result, mostly given by the b-jet mass resolution from the reconstruction of the h → b b decay of one of the Higgs bosons.A realistic value of a 50 GeV bin size was compared to more optimistic sizes down to 10 GeV.While the smearing has a visible effect on R, the binning hardly reduces its value.Similarly, the location of the bin, which is partially arbitrary, has a smaller impact on R once we take into account the finite resolution of the detector (smearing).As a last step we showed that the sign(λ hhH × ξ t H ) does not leave any measurable imprint on the m hh distributions, once the experimental uncertainties are taken into account.
We conclude that, depending on the values of the underlying Lagrangian parameters, a sizable resonant H contribution to the di-Higgs production cross section can leave possibly visible effects in the m hh distribution.This would pave the way for a first determination of a BSM THC, a step that is crucial for the reconstruction of the Higgs potential of the underlying BSM model.

Figure 1 :
Figure 1: Leading-order diagrams to SM-like Higgs pair production in gluon fusion processes at hadron colliders.The red dot shows the triple Higgs coupling, and h i = h, H.

Figure 2 :
Figure 2: Plane 1. 2HDM type I, tan β versus cos(β − α).Upper left: Cross section prediction for di-Higgs production in 2HDM normalized to the SM value, both evalued at NLO QCD in the heavy top-quark limit.Allowed area inside the black contour.Red lines indicate the values at with the ratio is 1 (i.e.σ = σ SM ).The color coding indicates the size of this ratio.Upper right: K-factor, defined as the ratio of the NLO and LO cross sections.Lower left: Total decay width of the heavy Higgs H. Lower right: Ratio of the cross section with and without resonant enhancement, both evaluated at LO.

Figure 6 :
Figure 6: Plane 1. Upper line: Triple Higgs coupling predictions in the 2HDM and value of the normalized cross section (w.r.t. the SM value) evaluated at NLO QCD in the same parameter space.Red (blue) dots show the maximum (minimum) value of the trilinears that is realized in the allowed region.Lower left: Points realized in the couplings plane, the red area around κ λ = 1 and λ hhH = 0 represents the points that fall into the allowed region, indicated by a red arrow.Lower middle: Zoom into the previous plot, color code indicates the normalized cross section within the allowed region for different values of triple Higgs couplings.Black dots indicate the existent scan parameter points in the κ λ -λ hhH plane within and outside the allowed region that was obtained from the figures in the upper row.The red star indicates the SM limit (κ λ = 1 and λ hhH = 0).Lower right: Expected sensitivity to the deviation of the cross section from the SM value.Red star indicates the SM.

Figure 9 :
Figure 9: Plane 3. Split for negative (left) and positive (right) values of c β−α .Upper line: Allowed region.Middle line.Total di-Higgs production cross section at NLO QCD w.r.t the SM.Lower line.Expected sensitivity to the cross section deviation (zoomed into the interesting region).

Figure 11 :
Figure 11: Left: Partial result σ interf Eq. (15) for three different decay widths.Right: Effect of Γ tot H on the the invariant mass distribution for one benchmark point and three values of Γ tot H (see text).

Figure 13 :
Figure 13: Analysis of points with large di-Higgs production cross sections in the benchmark plane 3. Left.Location of the benchmark points in the total production cross section plot.Right.Invariant mass distribution for selected points (with a bin size of 5 GeV).The colors of the points in the two figures are matching.The values of σ tot in the legend indicate the LO inclusive cross section prediction for each point.

Figure 14 :
Figure 14: m hh distributions for selected points in plane 3. Upper lef t: Total m hh distributions for the points of Tab.3; upper right (red), lower left (blue), lower right (black): Individual contributions to the m hh distribution for the three points: total cross section (red), continuum cross section (black), m hh involving only λ hhh (λ hhH ) in light (dark) blue, m hh involving both (no) THCs in purple (yellow).

Figure 15 :
Figure 15: Sensitivity to triple Higgs couplings for points with c β−α = 0.1015 in benchmark plane 3. Upper left: The common mass m ≡ m H = m A = m H ± and correspondingly m 2 12 change, as marked with orange, yellow, purple, garnet and green dots (in ascending mass order, color coding indicating σ/σ SM ).Upper right: Location of the benchmark points in the κ λ -λ hhH plane (color coding indicating σ/σ SM ).Lower plot: m hh distribution for the five benchmark points with a 5 GeV bin size.The colors indicate m H , the red arrows show the location of the "resonant peak".

Figure 16 :
Figure 16: Sensitivity to triple Higgs couplings for points with c β−α = 0.203 in benchmark plane 3. Upper left: m ≡ m H = m A = m H ± and correspondingly m 2 12 changes, as marked with orange, yellow, purple, garnet and green dots (in ascending mass order, color coding indicating σ/σ SM ).Upper right: Location of the benchmark points in the κ λ -λ hhH plane (color coding indicating σ/σ SM ).Lower plot: m hh distribution for the five benchmark points.The colors indicate m H .

Figure 17 :
Figure 17: m hh distributions for selected points in plane 3 with c β−α ∼ 0.2.Upper left: Total m hh distributions for the points of Tab.6; upper middle, left and lower left, middle and right: Individual contributions to the m hh distribution for the five points: total SM-like cross section in black, total cross section in red, m hh involving only λ hhh (λ hhH ) in light (dark) blue, m hh involving both (no) THCs in purple (yellow).

Figure 18 :
Figure 18: Interference contributions for the benchmark points in plane 3 with c β−α = 0.203.Left: Interference of the two triangle diagrams.Middle: Interference of the SM-like diagrams.Right: Interference between the box and the resonant diagram.Solid lines indicate that the interference is positive and dashed lines indicate that it is negative.

Figure 19 :
Figure 19: Sensitivity to triple Higgs couplings for points within benchmark plane 4. Upper left: The points have fixed m A = m H ± = 544.75GeV, c β−α = 0.2, whereas m H and correspondingly m 2 12 change, as marked with purple, orange, blue, light blue, green, yellow and garnet dots (in ascending mass order, color coding indicating σ/σ SM ).Upper right: Location of the benchmark points in the κ λ -λ hhH plane.Lower line: m hh distribution for the seven benchmark points with a bin size of 5 (20) GeV in the left (right) plot.The colors indicate m H .

Figure 20 :
Figure 20: Theoretical (blue) and smeared (red) invariant mass distributions for the selected benchmark point (see text).Solid (dashed) lines show the contribution of the total (continuum) differential cross section.Left (right) plot has a 10% (15%) smearing.

Figure 21 :
Figure 21: Different bin sizes for a 10% smeared distribution in the example benchmark point (10, 20, 40 and 50 GeV).The red lines correspond to the true (smeared and binned) prediction of the m hh distribution.The other color indicates the corresponding binned, but unsmeared distribution.Solid (dashed) lines represent the total (continuum) contribution to the cross section.The grey region represents the region that falls into the window defined to compute the variable R. The black vertical line indicates the value of the resonant mass, i.e. 512.5 GeV.

Figure 23 :
Figure 23: Invariant mass distributions for different bin locations assuming a bin size of 50 GeV and a smearing of 10% (left) and 15% (right).Solid (dashed) curves show the full (continuum) result.

Figure 24 :
Figure 24: Invariant mass distributions for different signs of (λ hhH × ξ t H ) (see the right plot of Fig. 12).

Figure 25 :
Figure 25: Invariant mass distributions for different signs of (λ hhH × ξ t H ) with a 10% (left) and 15% (right) smearing applied.The black line shows the continuum distribution.

Table 1 :
Fermion couplings in the 2HDM type I and II.

Table 2 :
Factors appearing in the couplings of the neutral Higgs bosons to fermions, ξ f h,H,A , and to gauge bosons, ξ V h,H,A , in the 2HDM of type I and type II.

Table 3 :
Selected points in benchmark plane 3 with large di-Higgs production cross section.

Table 4 :
Structure of the resonance.

Table 7 :
.7].Similarly, also Γ tot H grows with increasing m H .All corresponding numerical values can be found in Tab. 7. Selected points in benchmark plane 4 for m A = m H ± ∼ 545 GeV.

Table 8 :
Values of the variable R for the significance of the signal for different bin sizes for a 10% and 15% smeared distribution, see Figs.21-22.

Table 9 :
Values of the variable R for the sensitivity of the signal for different bin locations for a 10% and 15% smeared distribution and a bin size of 50 GeV.