Corrections to di-Higgs boson production with light stops and modified Higgs couplings

The Higgs pair production in gluon fusion is a sensitive probe of beyond-Standard Model (BSM) phenomena and its detection is a major goal for the LHC and higher energy hadron collider experiments. In this work we reanalyze the possible modifications of the Higgs pair production cross section within low energy supersymmetry models. We show that the supersymmetric contributions to the Higgs pair production cross section are strongly correlated with the ones of the single Higgs production in the gluon fusion channel. Motivated by the analysis of ATLAS and CMS Higgs production data, we show that the scalar superpartners' contributions may lead to significant modification of the di-Higgs production rate and invariant mass distribution with respect to the SM predictions. We also analyze the combined effects on the di-Higgs production rate of a modification of the Higgs trilinear and top-quark Yukawa couplings in the presence of light stops. In particular, we show that due to the destructive interference of the triangle and box amplitude contributions to the di-Higgs production cross section, even a small modification of the top-quark Yukawa coupling can lead to a significant increase of the di-Higgs production rate.


INTRODUCTION
A scalar resonance with mass of approximately 125 GeV has been detected at the run-I of the LHC [1,2]. Since its discovery, there has been a lot of effort in studying its properties. In particular, the main production rates and decay modes at the LHC have been analyzed, leading to results that are close to the ones predicted for the Higgs boson in the SM. The accuracy of each of these measurements is low and hence at present a departure of the SM properties may only be obtained by a combined analysis of all production and decay channels. A recent combined analysis of the Higgs data collected at run I by ATLAS and CMS [3,4] in different production channels was used to determine the best fit to κ i = g hii /g SM hii , the ratio of the Higgs couplings with respect to the SM predicted values. All relevant ratios κ i are consistent with unity at the 2σ level, although errors are still large and moderate deviations of the Higgs couplings with respect to the SM values are possible. In fact, the best fit values of κ i present moderate deviations with respect to the SM predictions, which allows the presence of BSM effects in Higgs physics.
The double Higgs production provides a probe for new physics. In the SM, at the leading order, the Higgs pair production process in gluon fusion, gg → hh, receives contributions from two different quark-loop induced amplitudes, corresponding to a triangle (gg → h * → hh), and a box diagram, (gg → hh), shown in Fig. 1, with top-quarks giving the main contribution. The amplitudes associated with the two diagrams interfere with each other destructively. When the di-Higgs invariant mass, m hh , is below the threshold for the top quarks in the QCD loop to be produced on-shell, m hh ≤ 2 m t , the amplitudes associated with these two diagrams only contains real parts, and the destructive interference leads to an exact cancellation of the total one-loop amplitude at the di-Higgs production threshold m hh = 2 m h . The resulting cross section in the SM is small and the statistical significance of the Higgs pair production process becomes very low, making this process very sensitive to possible deviations of these amplitudes from their SM values.
The triangle and the box diagrams are both very sensitive to the Higgs couplings of the colored particles running in the QCD loop. This means that even a small deviation of the Higgs coupling to the top quark with respect to its SM value may lead to considerable impact on their contributions to the di-Higgs production. In addition, the triangle diagram, where a single off-shell Higgs is produced and transforms into a pair of Higgs bosons, is correlated to the diagram of single Higgs production via gluon fusion, gg → h, and it is also proportional to the triple Higgs coupling, λ 3 , which provides an important information in probing the Higgs potential. Moreover, the box and the triangle loop amplitudes become very sensitive to new heavy colored particles running in the QCD loop. Therefore, the di-Higgs production detection becomes a very promising channel in probing new physics, being sensitive to various kinds of new effects .
In this article, we shall analyze the possible modifications of the di-Higgs production rate within low energy supersymmetry models. These models allow for the presence of new light colored particles coupled strongly to the Higgs, namely the stops. Moreover, the Higgs sector in these models is extended to include an extra Higgs doublet, in the Minimal Supersymmetric extension of the SM (MSSM), and an additional extra singlet within the Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM). This implies, in general, that the Higgs boson couplings will depend strongly on the mixing of this particle with the additional neutral Higgs states and may present small deviations with respect to the SM ones. In particular, as emphasized above, the departures of the Higgs coupling to top quarks and the triple Higgs couplings from the SM values, may have an important impact on the di-Higgs production rate. This can also provide a probe to the nature of the electroweak phase transition due to its close connection to the triple Higgs coupling as pointed out in [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].
The correlation of the di-Higgs production amplitude with the single Higgs production amplitude implies that the new physics contributions to the di-Higgs production channel are restricted by the Higgs production rate, which, as mentioned above, is bounded to be close to the SM prediction. This was stressed in Ref. [6], where they quantify the enhancements in the di-Higgs production due to stops while keeping gluon fusion single Higgs production rate, trilinear Higgs coupling (λ 3 ) and modification of the top quark Yukawa coupling with Higgs at or close to the SM values. Their results confirm that it is very difficult to achieve large deviations of the di-Higgs production at the LHC from the SM value under these constraints. In this work, we vary all three quantities within 2σ of the experimentally allowed ranges allowed by the combined ATLAS-CMS analysis [3,4], taking also into account the constraints coming from run-2 data. We also extremize X t at This work is organized as follows. In Sec. II, we develop a general understanding of heavy colored scalar particles contribution to the di-Higgs production cross section. As an example of heavy colored particle, we take the stop contribution to the di-Higgs production process for the range of value of the couplings allowed by the results of the combined ATLAS and CMS Higgs data at the run-I of the LHC, and comment on possible modifications induced by run-II data. In Sec. III, we discuss the light stop one-loop contribution to 2h production. We show that when the gluon, as well as the bottom and top Higgs couplings are allowed to vary within the range consistent with the best fit values of κ i , the 2h production may be greatly enhanced compared to the SM predictions. Moreover, we study the di-Higgs invariant mass distribution, showing that its study may lead to relevant information on the new particles contributing to the di-Higgs production rate. We reserve Sec. IV for our conclusions. In the Appendix A, we present the form factors we use to perform the full one loop calculation.

Modifications of the gluon fusion Higgs and di-Higgs production cross sections
In this section, we consider the modification to double Higgs production cross section in the presence of light stops. We begin by writing down the stop mass matrix.
The parameters m Q and m U are soft supersymmetry (SUSY) breaking mass terms of the left-handed and right-handed stops respectively, X t = A t − µ cot β is the stop mixing parameter, Q t = 2/3 is the top quark charge, T t 3 = 1/2, s W is the sine of the weak mixing angle and tan β is the ratio of Higgs vacuum expectation values (VEVs). Neglecting the small contributions from D-terms (m 2 z /3 m 2 t , m 2 Q,U ) we obtain: In the presence of light stops, in addition to the triangle and box diagrams with top-quarks in the loop, shown in Fig. 1, there are new diagrams contributing to the double Higgs production at the leading order, shown in Fig. 2. Diagram (1) and (2) is the SM contribution, which may be modified by departures of the top-quark Yukawa coupling and the trilinear coupling with respect to the SM values. Diagrams (3) to (8) represent the stop contributions. While the dimensionful trilinear coupling of the Higgs to the stops has a strong dependence on the Higgs mixing parameter X t , which can be larger than the stop masses, the quartic coupling (bilinear in both the Higgs and stop fields) is fixed by the square of the top-quark Yukawa coupling. As the LHC has determined that the stops are significantly heavier than the top-quarks, diagrams (5) and (6) lead to only small contributions to the di-Higgs production rate, since there is no source of parametric enhancement associated with them. Diagrams (3), (4), (7) and (8) depend relevantly on the stop mixing parameter and tend to give the most relevant light stop contributions to the Higgs production rate.
We have computed the leading order (LO) double Higgs production amplitudes depicted in Fig. 2, finding agreement with the results presented in Ref. [5]. The Standard Model next to leading order (NLO) corrections have been computed in the literature in different approximations [33][34][35][36][37][38][39][40]. Recently, the full NLO QCD corrections have become available [41]. Far less is known about the supersymmetric corrections. The NLO QCD corrections to the supersymmetric contributions have been calculated in the limit of vanishing external momenta [42], and they can be sizable with stops below the TeV scale.
Due to the lack of a full NLO calculation in the supersymmetric case, in this work we will compute the cross section in MSSM at LO, then present its ratio to the SM LO value. According to [42], the K-factor in the MSSM can be larger than the K-factor in the SM, and therefore the ratio of the cross sections can be further enhanced by the NLO corrections. In light of these results, our computation could be viewed as a conservative estimate of the possible enhancement due to supersymmetric particles. In order to numerically compute the supersymmetry effects on di-Higgs production we have modified the calculations presented in the public program MCFM-8.0 [43], by including the stop contributions and allowing also for possible modifications of the top-quark Yukawa coupling y t and the trilinear Higgs coupling λ 3 .
In order to generate relevant modifications to the double Higgs production cross section, the lightest stop mass should not be too far above from the weak scale. Therefore, the constraints on the stop masses coming from LHC searches put strong restrictions on the possible size of the supersymmetric contributions to the di-Higgs production process. Considering the direct production of stops at the LHC, the constraints on the stop masses depend strongly on the how stops decay, and on the masses of other particles in the decay chain. The most relevant constraints come from the region in which the difference between the stop and the lightest neutralino masses are larger than the top-quark mass. In such a case, one would expect to find a significant decay of the stop into a top-quark and a neutralino.
In the simplified models LHC considers, a stop decays one hundred percent either to a top and the lightest neutralino, or to a bottom and the lightest chargino, which then decay to a W + and the lightest neutralino. The final state is therefore a bottom, a W + and missing energy in both cases, but the kinematic distributions and efficiencies are different. The current constraints on the stop mass in this region of parameters are at least 500 GeV, and depend on the mass of the lightest neutralino, becoming stronger for larger mass difference of the stop with the lightest neutralino [44][45][46][47][48][49][50][51]. Moreover, the stop bound in the compressed regions is also of the order of 500 GeV [44,45,52,53].
With more complicated decay chains, the constraints on stops could be weaker than the 500 GeV limit reported by the LHC. For example, in the presence of a light stau [54,55], the decay chain of the stop ist → bχ ± 1 → b ντ → b ν τχ 0 1 . The final state is two b-jets, two τ 's and missing energy. As τ 's are difficult to detect at the LHC, the stop constraints in this scenario could be weakened significantly compared to that in the simple models described above. Another way to evade the constraints in the compressed region is to consider gauge mediation models [56]. In gauge mediation models, the lightest neutralino can decay to a photon and a gravitino. Then in the compressed region, all other decay products other than the photons are too soft, and the final states are two photons and missing energy. The current diphoton plus missing energy search is only focused on the high mass region of the squarks, and the limit for stops around 500 GeV or below is weak [57,58].
Long lived stops which dominantly decay through a RPV coupling λ ijkūidjdk can have weaker constraints from the LHC as well [59]. In such scenario, the long lived stops can decay into a pair of down-type quarks, which lead to a displaced dijet final state. Then by recasting the 8 TeV data [59], for cτ ∼ 0.1 mm, stops lighter than 200 GeV can be allowed, and for cτ ∼ 0.4 mm, stops around 400 GeV can be allowed. The heavy stable charge particle (HSCP) search in this scenario is weaker compared to the displaced dijet final state for low values of cτ [59][60][61].
Of course, the exact limit in the above three scenarios can be only obtained by doing a detailed recast of the current LHC data. The current recasting tools only include the data up to 2.3 f b −1 of the 13 TeV Run [62,63], and it is beyond the scope of this work to analyze the exact stop limit in those scenarios. In this analysis, we are going to consider stops as light as 300 GeV as a reference value, showing how the effects on the di-Higgs production cross section depend on the exact stop mass bound.
In addition to direct constraints, light stops also modify the single Higgs production cross section in the gluon fusion channel, which is well measured at the LHC [3,4,64,65]. The effective gluon Higgs coupling in the presence of light stops, and modified top Yukawa can be calculated from the low energy effective theory (EFT) approach. The leading contribution to the gg → h process can be obtained from the one-loop QCD beta functions of the heavy particles [66][67][68][69][70]. The effects may be understood from the contribution of heavy particles to the gluon kinetic term, namely where g s is the strong coupling constant and β i is the contribution of the particle of mass m i to the QCD β function. If the masses of the loop particles depend on the Higgs VEV, one can obtain the couplings of the SM-like Higgs by replacing the dependence on the VEV. by a dependence on the Higgs field, i.e. m i → m i (h). If the particles are much heavier than the relevant energy scale m h , then we can integrate out those particles and write the single Higgs production process, gg → h using a 1/m i expansion of the effective Lagrangian. By substituting h → h + v, the Taylor expansion of the QCD effective Lagrangian on the Higgs field leads to Therefore, one can get an approximation to the single and di-Higgs couplings to the gluon field strength. The effective Lagrangian for multi-Higgs couplings to gluons has been extended to N 4 LO [71]. However, these effective couplings are computed at zero momentum transfer and hence they lose validity, if the typical momenta are of the order or larger than the masses of the particles running in the loop.
Let us stress that in single Higgs production the relevant scale is given by the Higgs mass and therefore this expansion leads to a good description of the effective coupling of the single Higgs to the gluon field strength for loop particle masses of the order of or above the weak scale. In the case of double Higgs production, however, the relevant di-Higgs invariant mass scale may be much higher than the particle masses, and hence the effective field theory tends to fail for relatively light particles, leading to large corrections to the di-Higgs production process [6,8], as we discuss in the Appendix B. Hence, in our analysis, we shall use the full one-loop computation of the di-Higgs production rate.
From Eq. (2.4), the contribution of the new particles to the linear coupling of Higgs to gluons may be obtained by the first order expansion of QCD effective Lagrangian can be written as: where we have grouped particles with the same quantum numbers in a singlet mass matrix M i and β i denotes their common contribution to the QCD beta function.
The contribution from stops to the single Higgs production rate has been previously considered in the literature. Including possible modifications in the Higgs coupling to tops, the modifications in κ g is given by [80][81][82][83][84] (see Appendix B) The value of κ t is governed, at tree-level, by the mixing between the CP -even Higgs bosons.
In the MSSM, for instance, κ t cos α/ sin β, where α is the CP -even Higgs mixing angle. Close to the alignment (or decoupling) limit sin α − cos β and cos α sin β, and for moderate values of |µ| and tan β,X 2 t is very well approximated by X 2 t . Let us stress, however, that in the MSSM relevant deviations of κ t from one can only be obtained at low values of tan β, for which the loop corrections are insufficient to bring the Higgs mass in agreement with observations for stops at the TeV scale [72][73][74][75][76][77]. In general, even for κ t = 1, we will assume that the Higgs mass is not given by the MSSM relations, but fixed by additional D-terms, that could arise in gauge extensions of the MSSM [78], or F-terms, like happens in the NMSSM [79]. In the NMSSM, similar relations betweenX 2 t and X 2 t are obtained, and small deviation of κ t from one may be obtained for both heavy as well as light scalar singlets [83,84].
In addition, let us comment that κ t = 1.1 in general implies the presence of additional, relatively light, non-standard Higgs bosons. The additional CP -even Higgs could lead to a resonant double Higgs production if their masses are larger than 250 GeV (for an analysis of the resonant production in singlet extensions of the SM, see Ref. [87,88]). For instance, in the MSSM (with additional D-terms to fix the Higgs mass and induce a significant Higgs mixing), for a 350 GeV heavy CP -even Higgs, with low tan β 1, and relatively heavy stops, the gluon fusion cross section is of the order of several pb [89]. The branching ratio for such heavy Higgs to a pair of SM Higgs bosons, in the absence of light charginos or neutralinos, is around a few tens of percent, which leads to a cross section for pp → H → hh higher than the nonresonant double Higgs production. In this case, we expect to see a resonance in the m hh distribution.
When the heavy Higgs mass is increased to about 500 GeV, the gluon fusion cross section is around a pb, and the branching ratio of H → hh decreases to about a few percent as the tt channel opens up. Then the resonant production rate is comparable to the nonresonant production rate. Because of the destructive interference between the resonant double Higgs production diagram and the box diagram, we expect a dip-peak structure in the m hh distribution at the parton level. Whether this structure is visible at the LHC depends on how large the cross section is, which depends on model parameters as tan β, and mt, and the detector resolution.
We should also stress that these results are strongly model dependent. In the NMSSM, when light singlets are present, or when there are large splittings between the CP -even and CP -odd Higgs bosons, the dominant Higgs decay of the heavy CP -even Higgs boson is into non-standard Higgs states, and the resonant production of pairs of SM-like Higgs bosons is highly suppressed [84][85][86]. Moreover, even with the MSSM Higgs sector, in the presence of light neutralinos or charginos, the decay branching ratio into pair of SM-like Higgs bosons could be highly suppressed [83]. In this article, we shall concentrate on the nonresonant production of SM-like Higgs bosons, and analyze the impact of the mixing with additional Higgs bosons via the modifications of the top-quark Yukawa and trilinear Higgs self couplings.
As can be seen from Eq. (2.6), a small enhancement in the Higgs coupling to tops, as currently allowed by data [3,4], would not only enhance the top-quark Yukawa coupling, but also the stop contribution to the gluon fusion cross section. Let us stress however, that the run I indications of a high value of κ t [3,4] have not been confirmed by the current run II data [90][91][92] and hence in the following we shall consider only small variations of this coupling. Moreover, the stop effects may be significantly enhanced for large values of X t , which could also lead to a reduction of the Higgs coupling to gluons κ g , within the range allowed by the run I best fit values. However, a very large X t might also affect the Higgs vacuum stability [93][94][95][96][97][98]. In this paper, following the results of Ref. [98], we shall use the approximate bound when µ is small, A t X t . Here r = m 2 Q /m 2 U and the last term represents the impact of the CP -odd Higgs mass. In general, we will take a conservative approach and neglect the second term, containing the dependence on m A and m z , and consider the above bound on A t as a bound on X t , as would approximately arise from Eq. (2.7) at large values of tan β and moderate values of µ and m A . However, we will also discuss the impact of considering the above bound on A t for low values of tan β and moderate values of m A and µ.
In summary, we have calculated the modifications to double Higgs production with a modified version of MCFM in the presence of a light stop, modified top Yukawa and Higgs trilinear couplings. The mass of the light stop, the stop mixing angle, and the modifications to the couplings are subject to direct stop searches, precision Higgs signal strength measurements, and the vacuum stability. The numerical results for the double Higgs production cross section and the current experimental constraints will be discussed in the next section.

Collider Phenomenology
The collider phenomenology depends strongly on the precise stop masses, the stop mixing angle and the values of the top-quark Yukawa and trilinear Higgs couplings. In Fig. 3, we compute the variation of the di-Higgs production cross section in the absence of stops. As it is clear from this figure, even a mild variation of the top quark Yukawa coupling, κ t = 1.1, can lead to an increases of the cross section by 50 percent. The reason for this is that the contribution to the SM amplitude associated with the box diagram, which increases quadratically with κ t , is about a factor 2.5 larger than the one associated with the triangle diagram at the 2 m t threshold, which increases only linearly with κ t , and interferes destructively with the box amplitude.
We also show the variation of the cross section with a modification of the trilinear Higgs coupling. For κ t = 1, the value of λ 3 = 2.5 λ SM 3 maximizes the destructive interference between the box and triangle diagram amplitudes, and hence leads to a general reduction of the di-Higgs production cross section. On the contrary, for small values of λ 3 0, only the box diagram contributes, and hence the cross section is not only enhanced with respect to the SM case, but depends quartically on the top quark coupling κ t . Di-Higgs production cross section values of the order of 4 times the SM value may be obtained for the maximal variations of κ t and λ 3 considered in Fig. 3.
In the Figs. 4, 5, 6 and 7, we show the results for the double Higgs cross section in the presence of light stops. For each values of m Q and m U , we calculated the largest value of Figure 3: Di-Higgs production cross section in the absence of stops, as a function of the top-quark Yukawa coupling, κ t , for different values of the Higgs trilinear coupling λ 3 . Here, we have κ t = κ g .
|X t | that can be allowed by a lower bound on stop mass and a stable Higgs vacuum, with a Higgs vacuum expectation value of v = 246 GeV. The lower bound on the stop masses used in Figs. 4, 5, 6 and 7 are 400 GeV, 300 GeV,500 GeV, and 400 GeV respectively. Then we use the previously mentioned modified version of MCFM-8.0 to calculate the double Higgs production cross section, which is normalized to the SM value, as shown by the green dashed contours. For the stability condition, we decided to be conservative and ignore the m A and m z dependence in Eq. (2.7). The dependence on m A of the vacuum stability bound on X t , and of the resulting double Higgs production cross section, will be discussed later.
We also calculate the single Higgs production cross section in the gluon fusion channel, as shown in the orange regions. The left panels in all three figures correspond to a value of the top-quark Yukawa coupling normalized to the SM value, κ t = 1.0, while the right panel corresponds to κ t = 1.1. The modification of the triple Higgs coupling is defined as The first and last row in each of the Figs. 4, 5 and 6 corresponds to δ 3 = 0 and δ 3 = −1. The latter is used as an example to demonstrate the effect on the di-Higgs production in the case of no destructive interference between the triangle and the box diagrams. Such values of δ 3 can be realized in a scalar singlet extension of the SM Higgs sector, like the NMSSM, as demonstrated in [30]. Finally, Fig. 7 shows the case δ 3 = 1.5, in which, as stressed above, there is large destructive interference between the triangle and box diagram amplitudes.
Assuming the Higgs top-quark and triple Higgs couplings acquire SM values, the double Higgs cross section can be as large as 1.8 times the SM production cross section. When a small enhancement of the Higgs coupling to top-quarks is allowed, for instance κ t = 1.1, the double Higgs production cross section can be as large as about 3 times the SM cross section with a 400 GeV light stop. Allowing the modifications of the Higgs trilinear coupling, for δ 3 = −1 and κ t = 1, similar enhancement of the di-Higgs production rate is obtained. When both modifications are considered together, i.e., κ t = 1.1 and δ 3 = −1, then the values of more than 4 times the SM di-Higgs production rate may be obtained as shown in the bottom right panel of Fig. 4.
An important property that may be also extracted from Fig. 4, is the strong correlation between the modification of the di-Higgs production cross section and the one of the gluon fusion single Higgs production rate. A large modification in the di-Higgs production indicates a large value of the stop mixing parameter X t , which decreases the gluon fusion single Higgs production rate, Eq. (2.6). Values of the Higgs coupling to gluons somewhat lower than 0.9 times the SM value are required to obtain the largest corrections to the di-Higgs production rate for κ t = 1, while for κ t = 1.1, the required values are 0.9 ≤ κ g ≤ 1. These values of κ g are consistent with the combined fit to the run-I Higgs data, that leads to a preference for lower values of κ g 0.81 +0. 13 −0.11 [4].  As emphasized before, in the SM case, this value corresponds to a maximal reduction, due to the destructive interference between the triangle and the box diagrams, of the di-Higgs production [35,37]. This value can also be realized in singlet scalar extensions of the SM Higgs sector, like the one present in the NMSSM, and is strongly correlated with the obtention of strongly first order phase transition, which is of particular interest from the perspective of EW baryogenesis, as shown in [20,22,25,30,99,100]. In the absence of stops and κ t = 1, the double Higgs cross section falls to ∼ 0.4 times that of the SM value [35,37], making statistical significance too low even at the end of the LHC run. In Fig. 7, the left panel shows that the addition of light stops, with the lighter one having a mass of ∼ 400 GeV, can more than double this value to make it about 0.9 times the SM cross section. In the right panel, we show that in addition to a ∼ 400 GeV stop, small modifications (∼ 10%) of κ t can increase the double Higgs cross section to about 40% above the SM value.  cross sections obtained for light stops and modified couplings are approximately given by the ones obtained for light stops and SM values of the couplings times the ones obtained by modifying the couplings in the SM case, shown in Fig. 3. Also, for Figs. 3,4,5, 6 and 7, as we go to the region of high values for both m Q and m U (of about or more than 1 TeV), the vacuum stability constraints on X t constrain the lightest stop mass to values larger than the minimum allowed one, and therefore the di-Higgs production cross section results become independent of the allowed lightest stop mass. Thus, the mass of the lighter stop increases for higher m Q and m U , leading to the decoupling of the stop effects. We also compare the full one-loop calculation(solid lines) with the EFT calculation (dashed lines) in Fig. 8 (for the detailed EFT calculation, see Appendix B). We chose κ t = 1 for orange, red and green lines, and κ t = 1.1 for the blue lines. The value of κ g increases monotonously along each line except the orange lines with increasing mass of the lightest stop. The range of κ g values are from 0.85 to 0.98 for the red, 0.77 to 0.98 for the green and 0.90 to 1.08 for the blue line. For a given lightest stop mass and κ t , the κ g value is same for EFT and one-loop case. The orange plots have κ g = 1 by definition. All lines are plotted taking m Q = m U and δ 3 = 0.
In all cases, when the stops are heavy enough, i.e., above ∼ 1 TeV, the EFT calculation and the full one-loop calculation agree well, as one would expect. Also, with heavy stops, the cross section ratio approaches one for κ t = 1 cases, and 1.6 for the κ t = 1.1 case (blue line), which is in agreement with the results shown in Fig 3. The orange line depicts results for κ g = 1, namely when X 2 t is chosen to be m 2 (see Eq. (2.6)). In the red and blue lines, instead, X 2 t is chosen to saturate the vacuum stability condition, Eq. (2.7), in a conservative way by neglecting the m A and m z terms, and the results are in agreement with the ones shown in the top panels of Figs. 4,5 and 6, where a similar vacuum stability constraint was considered. Finally, the green lines show the dependence of the double Higgs production cross section on the stability bound on X t . In order to show this dependence, we use m A = 350 GeV, µ = 400 GeV, and tan β = 1 to calculate the saturation of the stability condition, Eq. (2.7). Therefore a larger X t can be allowed, and larger modification to the di-Higgs production cross section can be achieved. For instance, for a lightest stop mass of 500 GeV, the bound on X t 3 m Q , instead of X t 2.6 m Q that is obtained when one neglects the m A dependence of the vacuum stability bound, and the modification of the di-Higgs production cross section can be as large as 60 %, compared to about 30 % when no m A dependence is considered.
In Fig. 9, we show the effect of stop mixing parameter X t on the di-Higgs production cross section for a fixed value of the mass of the lighter stop in the case of m Q = m U and δ 3 = 0. Red, green and blue represent fixed lighter stop mass of 300, 400 and 500 GeV respectively. Solid lines correspond to κ t = 1, while dashed lines correspond to κ t = 1.1. The maximum value of X t for green and blue lines correspond to the condition X t ∼ 3m Q as set by taking m Q = m U and suitably high experimentally allowed values of m A in Eq. (2.7). The maximum value of X t used for the red line is less than 3m Q in order to increase the readability of the plot (X t < 2.88m Q and X t < 2.78m Q for the solid and dashed lines, respectively).
For lower values of X t and κ t = 1, the di-Higgs production cross section becomes smaller than the SM one and then starts increasing owing to the sign flip for the trilinear coupling of lighter stops with the Higgs boson, which becomes linearly dependent on X t . Therefore, for large values of X t the contribution to the amplitude that grows quadratically with X t (last row of Fig. 2) becomes dominant and the cross section grows proportional to X 4 t . This behavior is clearly seen in the red line in Fig. 9, corresponding to a lightest stop mass of 300 GeV. For the green and blue lines the vacuum stability condition of X t < 3m Q puts an upper bound on X t /mt 1 , and makes other contribution to the amplitude competitive with 0 400 800 1200 1600 2000 2400 0 1.

7.
8. Figure 9: X t dependence of the Di-Higgs production cross section for δ 3 = 0 and m Q = m U for lighter stop masses of 300 GeV (Red), 400 GeV (Green) and 500 GeV (Blue) GeV after neglecting the D-terms. Solid lines correspond to κ t = 1, while dashed lines correspond to κ t = 1.1. The condition for vacuum stability described in Eq. (2.7) is taken to be X t < 3 m Q , which can always be achieved for a suitably high value of m A allowed experimentally. For presentational reasons, for the red line is cut at values of X t smaller than 3 m Q (see the text).
those that depend quadratically on X t , preventing the X 4 t behavior to develop.
As can be seen from Fig. 9, the increase in X t leads to significant enhancements in the cross section. For lighter stop masses as low as 300 GeV, large enhancements by a factor of order ten can be obtained before the vacuum breaking condition X t < 3m Q is met. Even for experimentally more viable values of the lighter stop mass such as 500 GeV, considerable enhancements of 60% and 140% are possible for κ t = 1 and κ t = 1.1, respectively.

Di-Higgs Search Channel
The general strategy in the search for double Higgs is to require one Higgs to decay to a pair of bottoms for enough statistics, as the total rate for double Higgs production is about three orders of magnitude smaller compared to single Higgs production. Then, we can consider the other Higgs decay to a pair of photons, bottoms, W ± 's, or τ 's. In this work, we are going to discuss the modifications to distributions in the presence of light stops, and we will focus on the bbγγ channel, as this channel provides the best resolution.
The cross section for bbγγ final state depends not only on the di-Higgs production cross section, but also on the Higgs decay branching ratios to bb and γγ. These decay branching ratios depend strongly on the Higgs coupling to W ± gauge bosons and bottom quarks, called κ w and κ b respectively. It is then important to see what are the values of κ w and κ b allowed by Higgs data, for the values of κ t and κ g considered in this work. In order to do this, we recall that, the gluon fusion production rate is modified by a factor κ 2 g , while the vector boson fusion and associated production with vector boson channels are modified by κ 2 w and the tth channels are modified by a factor κ 2 t . Moreover, the modified branching ratios are given by where BR(h → XX) is the branching ratio of the Higgs decay into a pair of X particles.
In Fig. 10, we fix κ t at 1 or 1.1, κ g at 0.80, 0.90 or 1, which are representative values for the gluon Higgs couplings necessary to obtain sizable modifications of the di-Higgs production cross section. Having fixed these values, we fit for the preferred values of κ b and κ w . We include all the Higgs data from Run I [3,4], except h → ZZ * and h → τ τ , as they mostly depend on κ z and κ τ , which are beyond the scope of discussion in this study. The production for VBF also depends on κ z , which we fix at the run I best fit value κ z = 1. Due to the small value of the BR(h → ZZ), fixing κ z = κ w makes no difference in our results. The value of κ γ is considered to be consistent with the values induced by the presence of light stops and modifications of κ t and κ w . Using effective field theory to evaluate the top and stop contributions, one obtains, approximately where we used Eq. (2.6) and the fact that the relation between the top and stop contributions to κ g and κ γ are the same.
The region within 1 σ of the best fit value for κ b and κ w is shown in blue, and the region within 2 σ of the best fit value is shown in light blue. Then, for given values of κ w and κ b , we calculate the Higgs decay branching ratios to bb and γγ, and we show the contours of BR(h → bb) × BR(h → γγ), normalized to the SM value. We also show the Run 2 results for gluon fusion, h → γγ in orange(ATLAS) [101], and green(CMS) [65]. The solid lines are the central values, and the dashed line show the 1 σ range. The region above the dotted line is consistent with the Run 2 measurement of associated production of Higgs with vector bosons, V h, with h → bb within 1 σ [102]. It can be seen from the top two panels that κ t does not change the fit, as the tth channel has large uncertainties. κ t does not change the branching ratios either, because by allowing new particles in the loop, or considering κ g as an independent parameter, κ t does not change the Higgs decay. Then for κ g = 0.9 and κ g = 1 we only consider κ t = 1.
Our results are roughly consistent with the ones obtained by the combined ATLAS and CMS Higgs data [3,4]. As can be seen from these contours, some small modifications to BR(h → bb) × BR(h → γγ) are expected, which would modify the hh → bbγγ rate. However, the largest modification is about ±20%. Let us stress that the inclusion of run II data is likely to move κ b towards larger values. However, as is apparent from Fig. 10, this modification is unlikely to modify the above conclusion. Therefore, only mild variations are expected in the product of the bb and γγ decay branching ratios and the hh → bbγγ rate is mainly controlled by the modifications of the di-Higgs production rate with respect to the SM value.

Modifications of the di-Higgs invariant mass distribution
As pointed out in [30,103], a modification in λ 3 can lead to a drastic change in the kinematic distributions for double Higgs production. The m hh distribution shifts significantly to lower values for δ 3 2. In this section, we study the possible modifications to the m hh distribution with modified κ t and the presence of a light stop.
As we emphasized before, in the SM, there are two diagrams contributing to the double Higgs production, the box diagram and the triangle diagram. The two diagrams interfere with each other destructively. The cross section for the triangle diagram scales with κ t as κ 2 t , and the cross section of the box diagram scales with κ t as κ 4 t . Therefore, a modification in the top Yukawa can change not only the di-Higgs production cross section but also the m hh distribution. However, without modifications in λ 3 , the box diagram dominates over the triangle diagram, and as only a few tens of percent deviation is allowed in the top Yukawa, we do not expect that this change can modify the m hh distributions in any relevant way, as we have checked in our numerical simulations and can be seen from the blue dashed line in Fig 11. The modification of the invariant mass distribution could become relevant when the cancellation between the two diagrams become strong. This occurs for values of λ 3 ∼ 2.5, for which the amplitudes associated with the two diagrams become comparable in size. In this case, a cancellation of the production rate appears at some value of m hh . Then modifications of κ t would change the relative weight of the triangle and box diagrams, inducing a more relevant change in the invariant mass distribution, m hh . This can be seen from the green dot dashed line and the magenta solid line in Fig. 11, where λ 3 is 2.5λ SM 3 in both lines. When κ t is 1 (green dot dashed line), the cancellation is at m hh about 2m t . As κ t increases to 1.1 (magenta solid line), the box diagram increases more than the triangle diagram, and the exact cancellation occurs for smaller values of m hh , at about 330 GeV.
Furthermore, in the presence of a light stop, the amplitudes for diagram (3) -(8) in Fig. 2 develop imaginary parts when the invariant mass m hh crosses the 2 mt threshold, inducing a second peak in the m hh distribution a little above 2 mt. We selected benchmarks to study the distributions of m hh with a light stop, and possible modifications in λ 3 and κ t . The benchmarks are listed in Table 1, in which we are neglecting D-terms when calculating    Table 1: Benchmarks points for light stops giving a sizable correction to the di-Higgs production cross section at hadron colliders and therefore the second peak is around 640 GeV and 1 TeV, respectively.
As can be seen from Fig. 12, the kinematic distributions are similar to the ones in the SM and hence the m hh cut efficiency in this benchmark points is similar to the SM case.
Other kinematic variables that have been used at the LHC, including the invariant mass distributions of the bottom quarks, m bb , and of the diphotons, m γγ , as well as cuts on the p T of the b-jets, and the photons are expected to have a similar behavior as in the SM. Therefore, the projected sensitivity scales approximately with the signal rate and therefore we use the ATLAS SM results to estimate the projected sensitivity for our benchmarks at the High Luminosity run of the LHC (HL-LHC) [104], with a projected luminosity of 3 ab −1 , in Table 2. CMS shows a similar sensitivity in this channel [105].
A recent work proposes to use the log-likelihood ratio to identify kinematic regions and shows an improved sensitivity [106]. As one can see from Table 2, only using the bbγγ channel, the HL-LHC will be sensitive to light stops with a large mixing, which can be a indirect probe for light stops regardless how the stops may decay. For stops as heavy as 500 GeV, the LHC sensitivity is limited to the cases of a large mixing, a negative correction to the Higgs trilinear coupling, which is well motivated by a strong first order phase transition,  and/or a small positive correction to the top-quark Higgs coupling, such as it appears in benchmark point D.
In summary, the presence of a light stop, modifications of the top Yukawa and trilinear Higgs couplings can lead to sizable contribution to double Higgs production. The stop contributions are summarized in Fig. 8, the contributions from a modified top Yukawa coupling is summarized in Figs. 3. We present some benchmarks and their projected sensitivities in Table 1 and Table 2.

Conclusions
The search for di-Higgs production is one of the main goals at hadron colliders. This is due to the sensitivity of this channel to new physics and its dependence on the Higgs potential parameters. The sensitivity of the LHC experiments to this channel is limited by the small rate and large backgrounds in the main final state channels. It is therefore very important to study under which conditions the di-Higgs production rate may be enhanced, allowing for its study at a high luminosity LHC. Barring the possibility of resonance di-Higgs production via the presence of heavy scalars decaying into pairs of SM-like Higgs bosons, it is known that this can be done in the presence of negative corrections to the trilinear Higgs coupling and/or positive corrections to the top quark coupling to the Higgs. In this work we emphasized the strong dependence of the di-Higgs production cross section to small, positive corrections of the top-quark coupling to the Higgs, which are still allowed by the current LHC Higgs data.
Furthermore, we studied the additional effects of light stops on the di-Higgs production cross section. We computed the one-loop corrections associated with light stops, finding agreement with previous expressions in the literature. We then incorporated these corrections into a modified version of the program MCFM-8.0, including the possibility of light stops together with possible modifications of the top-quark and trilinear Higgs couplings. We found out that large corrections to the di-Higgs production rate are possible in the case of relatively light stops, with a large stop mixing parameter. The effect of light stops may become even stronger under modifications of the top-quark or trilinear Higgs couplings.
In general, we found that the modifications of the di-Higgs production rate are strongly correlated with similar modifications of the gluon fusion Higgs production rate, and can significantly enhance the LHC sensitivity to this production channel. Moreover, we also found that the precise constraints on the trilinear Higgs-stop coupling coming from the requirement of vacuum stability, have a major impact on the size of the possible stop corrections for lightest stop masses above 500 GeV, that is the current stop bound on the stop mass in standard decay channels.

A Form Factors
In this section, we present scalar form factors we use. The cross section of gg → hh is determined by two Lorentz and gauge invariant structure functions given in in [5,[107][108][109][110][111].
The differential cross section reads, Here the matrix-element M are separate by four parts according to the helicity of the incoming glouns, where + and − denotes the right-and left helicity gluons. Thus [112]: The amplitude can also be written in terms of form factors of triangle F , and box diagrams, and quartic coupling C 2 normalized to unity [113], where the two terms inside the bracket in Eq. (A.3) are in one to one correspondence with the ones given in Eq. (A.2) and G F = 1 √ 2v 2 is the Fermi constant. We will discuss the value of the above form factors below.
The above amplitudes depend on the couplings g hab of the mass eigenstates of quarks and squarks to the Higgs field. The stop mass eigenstate couplings may be obtained from the corresponding coupling of the weak eigenstates, namely where R is the rotational matrix which rotates the left-and right-handed squark fields to the mass eigenstates. [114]. The mixing angle θ of the 2 × 2 rotation matrix : Considering only small deviations of the Higgs couplings to fermions, close to the decoupling limit α β − π/2 and neglecting the small contribution of the D terms, the linear and bilinear couplings of the stop mass eigenstate to the Higgs are given by The form factors are associated with the diagrams given in Fig. 2, where the Higgs triple coupling constant × T C h0 ijj (t) + C h0 jii (t) + U C h0 ijj (û) + C h0 jii (û) + 2ŝC hh iji (ŝ) ijji (t,û) + 2ŝm 2 t j D h0h0 ijji (t,û) + D hh00 jijj (ŝ,t) + D hh00 jijj (ŝ,û) (A.20) × D h0h0 jijj (t,û) + D hh00 ijii (ŝ,t) + D hh00 ijii (ŝ,û) Here p 1 , p 2 are the momentum of incoming glouns, k 1 , k 2 outgoing Higgs, with p 2 1 = p 2 2 = 0, k 2 1 = k 2 2 = m 2 h . U = (m 2 h −û), T = (m 2 h −t) and S = (m 2 h −ŝ). C ab ijk (α) and D abcd ijkl (α,β) are defined in terms of the Passarino-Veltman functions C and D, which are given in [115]

B Effective Field Theory Analysis
In the limit that the colored particles in QCD loop is much heavier than the relevant energy scale in the theory, m Q m hh , the form factors in (A.3) may be computed using effective field theory techniques for the effective vertices shown in Fig. 13. In Sec. 2, we discussed the leading contribution to the single Higgs production. For the di-Higgs production, according to Eq.(2.4), the second order coupling of the Higgs, necessary for the computation of the gg → hh amplitude, can be written as: We can now define dimensionless couplings g In the limit of vanishing soft supersymmetry breaking terms, these form factors become constants, proportional to the particle contributions to the QCD β function [8], In the large mass approximation, it is instructive to use EFT to consider the new particles implication to gg → hh process.
If we now consider BSM modification to the di-Higgs production process, we have the parton cross section: h /(β t g t h ) is the coupling of the Higgs to a pair of gluons, normalized by the SM induced one. Once the set of couplings, κ, is introduced to parametrize deviations from the SM couplings of the Higgs bosons couplings to SM bosons and fermions [116], one could compute the Higgs cross section obtain by the annihilation of a particle i by κ 2 i = σ i /σ SM i . Therefore, using Eq. (2.1) and (2.2), the couplings gt h and gt hh are given by : The coupling of single Higgs boson production via gluon fusion process gg → h, κ g , is defined to be 1 in the SM based on the contribution form top and bottom quarks, and if we consider the contribution of top squarks running in the QCD loop, it would be given by : This could still be taken as a good approximation, and we use this equation to constrain the values of κ g we used in this paper.