Putting a Stop to di-Higgs Modifications

Pair production of Higgs bosons at hadron colliders is an enticing channel to search for new physics. New colored particles that couple strongly to the Higgs, such as those most often called upon to address the hierarchy problem, provide well motivated examples in which large enhancements of the di-Higgs rate are possible, at least in principle. However, in such scenarios the di-Higgs production rate is tightly correlated with the single Higgs production rate and, since the latter is observed to be SM-like, one generally expects that only modest enhancements in di-Higgs production are allowed by the LHC Run 1 data. We examine the contribution of top squarks (stops) in a simplified supersymmetry model to di-Higgs production and find that this general expectation is indeed borne out. In particular, the allowed deviations are typically small, but there are tuned regions of parameter space where expectations based on EFT arguments break down in which order 100% enhancements to the di-Higgs production rate are possible and are simultaneously consistent with the observed single Higgs production rates. These effects are potentially observable with the high luminosity run of the LHC or at a future hadron collider.


Introduction
A comprehensive experimental program to characterize the 125 GeV Higgs boson [1,2] and determine the underlying nature of electroweak symmetry breaking is underway at the LHC. Based on the complete Run 1 data set, significant progress has been made through the study of final states with a single Higgs particle. The largest and best measured single Higgs production channel is the one loop gluon fusion process, which is in good agreement with the predictions of the Standard Model (SM) [3,4]. In addition, an important long term goal of this program is to observe and study final states with two Higgs bosons. The di-Higgs channel is sensitive to the trilinear self coupling of the Higgs particle, which in turn gives information about the shape of the scalar potential, and can furthermore provide a sensitive probe of physics beyond the SM (BSM). Like single Higgs production, the dominant di-Higgs production channel at the LHC is gluon fusion, which is depicted in Fig. 1. In the SM and its extensions, di-Higgs production probes a different combination of couplings and masses than other loop processes such as single Higgs production via gluon fusion. One could then imagine that, even if for some reason the hGG coupling were SM-like, there could be large deviations in di-Higgs production. This expectation is further motivated by the fact that in the SM the two diagrams 1 of Fig. 1 interfere destructively making the SM di-Higgs production cross section smaller than the naive expectation [5][6][7]. Thus, typical BSM scenarios provide ample opportunity for significant modifications of di-Higgs production at hadron colliders when this cancellation is spoiled. Indeed this is the case for models with modified electroweak sectors and/or models where di-Higgs production may be resonantly enhanced through the production of new heavy fields which decay to Higgs pairs.
In this work we instead focus on another potential source of modifications: new colored fields that couple to the Higgs. We investigate how much these scenarios may modify di-Higgs production at the LHC and future hadron colliders through their impact on the momentum-dependent hGG and h 2 GG vertices (shaded green in Fig. 1) while keeping the Higgs quartic coupling λ at its SM value. Throughout we refer to these as 'non-resonant' corrections.
As a first step, consider the effective field theory (EFT) below some cutoff Λ for the Higgs-gluon couplings hGG and h 2 GG. In general, if new heavy colored fields that couple to the Higgs and have mass m ∼ Λ are integrated out, they generate operators of the form where H is the Higgs doublet in the unbroken theory. In the broken theory, we can write the operators in terms of the physical Higgs field, h, and work to quadratic order in h. If we also include the SM contribution, which we will denote with the coefficient c SM α s /12π, we obtain the effective operators 2) where v = 174 GeV and the sign flip between single and double Higgs couplings in the SM has been included.
We now introduce a core observation from the first run of the LHC: modifications to the total single Higgs production rate are small. For a model in which the only BSM physics is new colored fields coupled to the Higgs, the cross section modifications must be O(20%) [3,4], implying modifications to the hGG coupling of O(10%). We may understand the implications of this observation for non-resonant contributions to di-Higgs production by studying Eq. (1.2) more closely.
If the new physics is heavy and respects decoupling, the usual rules of EFT apply. In particular, small corrections to single Higgs production imply c 1 v 2 /Λ 2 c SM and we can safely ignore the higher order terms. Then, Eq. (1.2) implies that the magnitude of corrections to the h 2 GG coupling must also be small if corrections to hGG are small, as the magnitude of both are controlled by the same parameter combination c 1 v 2 /Λ 2 c SM . Thus, we should expect non-resonant contributions in both diagrams for the di-Higgs production amplitude of Fig. 1 to be similarly suppressed.
It is worth noting that the impact of non-resonant new physics generically exhibits constructive interference between the triangle and box diagrams, unlike the top contribution in the SM. This implies that non-resonant corrections to di-Higgs production may spoil the cancellation in the SM and be larger than corrections to single Higgs production, but this will not be a very large effect.
Quite generally then, the constraint that the hGG coupling be SM-like implies that models with only colored, non-resonant, BSM states will have fairly SM-like di-Higgs rates for regions where the EFT is valid and current single Higgs constraints are taken into account. Clearly, the best chance for large deviations in the SM di-Higgs rate in this scenario is that the new particles are somewhat light so that an EFT analysis is inapplicable. In this case models must be checked on a case by case basis. In this work we explore this possibility in the context of scalar top partners (stops) in a simplified model as a supersymmetric extension of the SM.
Supersymmetry is attractive because it provides a solution to the hierarchy problem. In a natural SUSY model one expects stops with masses below the TeV scale. Such stops have been searched for directly at colliders, but these searches depend strongly on the superpartner spectrum and specific decay modes of the stop. The bounds on stops decaying to a top and neutral LSP are approaching the TeV scale when the LSP is light [8][9][10][11], and are expected to get stronger with future LHC data [12,13]. The bounds on very light stops, with masses in the 100 -200 GeV range are much more difficult to evade. One possibility is that the stop could decay in a way that makes it much harder to discover at a collider. For example, it could be stealthy and nearly degenerate with the top [14][15][16][17][18][19], or part of a compressed spectrum such that it is heavy but approximately degenerate with the particle it decays to [20][21][22][23][24][25][26][27], 2 or decay into other light MSSM particles (e.g. staus [30,31]), or decay via baryon number R-parity violation [32][33][34] where LHC searches are just starting to become sensitive [35]. Because stops can be hidden in various exotic decay modes, complementary indirect bounds on top squarks are a crucial tool in the exploration of weak scale SUSY.
Indirect probes of stops include modifications to the W mass [36,37], corrections to Higgs production rates and branching ratios [38,39] in loop processes, Higgs kinematic distributions [40,41] especially at high p T , effects on Higgs wavefunction renormalization [42,43], and stop-onium resonances [44][45][46][47]. Stronger constraints could be obtained with future colliders [48,49]. Because these probes of new physics are indirect, if a deviation is found it will be difficult to solve the inverse problem: what is the nature of the new physics that modifies a particular observable? Therefore, it is very important to explore as many different complementary probes as possible.
Higgs pair production has been studied in the Minimal Supersymmetric Standard Model (MSSM) [50,51], with [52,53] exploring the effects of scalars in loops. In this paper we show, using stops as a concrete and well motivated example, that the absence of large deviations in single Higgs gluon fusion makes it very difficult to generate large enhancements in double Higgs production from non-resonant contributions alone. We show this in the context of an effective field theory and also with stops using low energy theorems [54][55][56][57] as well as with a full loop calculation [52,53]. Despite these considerations, we do find that current Higgs data allow small, tuned, regions of parameter space with O(1) deviations in the di-Higgs total cross section.
In the following section, we survey the experimental and phenomenological literature on di-Higgs production at hadron colliders. While there is still significant uncertainty, we use it to select sensitivity benchmarks that we will use in this study. In Sec. 3, we analyze generic (and decoupling) heavy physics contributions to di-Higgs production using effective field theory, while in Sec. 4 we analyze heavy stops in the non-decoupling regime using low energy theorems. Finally in Sec. 5 we do a full loop calculation which is necessary for the case of light stops, and we find regions of parameter space where di-Higgs production has potentially observable modifications which are nonetheless consistent with single Higgs production constraints from Run 1. We conclude in Sec. 6, and we give results for a 100 TeV collider in the appendix.

Collider Phenomenology
We begin by reviewing the prospects to measure the di-Higgs channel at the LHC and future hadron colliders. Due to its importance in understanding electroweak symmetry breaking, di-Higgs production is a well studied channel. In the SM the di-Higgs production rate was calculated long ago [58,59], and at LHC energies the gluon fusion channel (see Fig. 1) dominates [60]. This process was computed at leading order (LO) [50,60] and next-toleading order (NLO) in the heavy top limit [61], with more recent computations including higher orders in 1/m t [62][63][64], parton shower effects [65], and virtual corrections [66]. There are also computations of di-Higgs plus one jet [67][68][69] and vector boson fusion (di-Higgs plus two jets) [70]. The computations continue to improve, but due to the difficulty of the final state, the uncertainty in projecting the collider reach in this channel is dominated by experimental challenges.
With Run 1 data, ATLAS has released a search for non-resonant di-Higgs in the bbγγ channel [71] setting a limit three orders of magnitude above the SM prediction. 3 There are also resonant searches in the 4b channel from CMS [72] and ATLAS [73], and in the bbγγ [74] and the multi-lepton/photon channel [75] from CMS, all of which have cross section limits that are O(pb), while the pair production cross section in the SM at 8 TeV is O(fb). Future projections depend very strongly on the projections for experimental efficiencies and systematics. Preliminary studies for high luminosity LHC at ATLAS [76] and CMS [77] in the bbγγ channel and CMS in the bbW W [77] show a marginal sensitivity to observing pair production with 3,000 fb −1 at 14 TeV, but further studies are ongoing.
There are also phenomenological studies that are more optimistic about the reach, but their sensitivity estimates vary greatly, even among those considering the same channels. For the most studied channel, bbγγ [78][79][80][81][82][83] significance estimates span from about 2σ to 6σ. Other channels, including bbτ τ [67,79,84], bbW W [67,79,85], and 4b [67,86,87] have similar qualitative variance in the observability of these channels. Therefore, we take uncertainty benchmarks of 30% and 60% for observing deviations from the total SM rate, but ultimately more study will be needed to determine the true sensitivity of future searches.
It is important to note, however, that di-Higgs modifications from stops will also lead to a modified spectrum in the di-Higgs invariant mass m hh or p T . Thus, to obtain the strongest possible limit one would ideally perform an analysis which is sensitive to not only the total cross section but also the spectrum, especially features at higher center of mass energies. Such an analysis would depend heavily on the final state which is being observed. Therefore, instead of a full shape analysis for a specific final state we consider two invariant mass bins to demonstrate the importance of considering the spectrum.
If loops of new particles such as stops are responsible for a modification to the di-Higgs total rate, then other di-Higgs production channels will have SM-like rates and can be used to disentangle new physics scenarios. Vector boson fusion is a large component of di-Higgs plus two jets. This channel has been studied [70,79,88] but because of the small cross section, it is quite challenging at the LHC. Higgs pair production in association with tt is another challenging channel [89,90], but perhaps a combination of these channels in conjunction with improvements in collider analysis could yield sensitivity in the future. Di-Higgs production has also been explored for physics beyond the SM, both in the context of effective field theory [82,83,[91][92][93][94], as well as for various specific new physics models [7, 50-53, 57, 95-113].
Planning is underway for higher energy hadron colliders where the cross section for Higgs pair production increases and prospects for measurements are potentially dramatically improved. The details of any putative collider and detector are still largely uncertain, but there have been several phenomenological studies of this process. The bbγγ [82,114], 4W [115], and bb + leptons (and possibly also photon or missing energy) [116] all appear to be promising ways to measure di-Higgs production at a 100 TeV collider. In App. A we consider modifications to di-Higgs production due to stops for a 100 TeV proton-proton collider, taking precision benchmarks of 10% and 20% on the rate.
Finally, as we have emphasized, it is useful to compare the process gg → hh to gg → h. The fitted rates for single Higgs production in gluon fusion, normalized to the SM value, are 0.85 +0. 19 −0.16 at CMS [3] and 1.23 +0.23 −0.20 at ATLAS [4], so we take the current bound to be 20%. These bounds will improve in the future, but ultimately will be systematics limited because of uncertainties in the SM prediction as well as experimental complications. With 3,000 fb −1 , the expected error on the coupling is 3-5% [117], so we take the ultimate expected error on the rate (twice the error on the coupling) to be 10%.

EFT Modifications to di-Higgs Production
In this section, we consider the generic effects of new heavy colored particles on di-Higgs production from an EFT perspective. When integrated out, these states will induce the effective operators presented in the introduction in Eq. (1.2). We can then write the relevant couplings contributing to di-Higgs production as Here we have defined the relative coupling shifts induced by the higher dimension operators defined in Eq.
We would like to understand the extent to which these coupling shifts can modify the di-Higgs production rate while being consistent with the observed SM-like single Higgs production.
The total di-Higgs production cross section can be written as Here the gluon parton luminosity is defined as where f g (x, Q) is the gluon parton distribution function, with factorization scale Q. Throughout this paper we use the MSTW [118][119][120] parton distribution functions when calculating the hadronic differential cross sections, with renormalization and factorization scales set to the invariant mass of the di-Higgs system. The partonic cross section in Eq. (3.2) is given byσ With the couplings in Eq. (3.1), the function A(ŝ) is given by Let us consider the case in which the new heavy colored states decouple from the Higgs as their mass is raised. This will happen if these states primarily obtain their mass from sources other than electroweak symmetry breaking. In this case, there is a separation of scales, v Λ, and the EFT expansion in Eq. (1.2) is a useful one. The leading dimension 6 operator dominates over the dimension 8 (and higher) operators, c 1 v 2 /Λ 2 c 2 v 4 /Λ 4 and there is a well-defined relation between the single and double Higgs production rate via gluon fusion in terms of the parameter κ h 1 , which is κ h 1 = −κ hh 1 . In Fig. 2 we plot the ratio of the di-Higgs production cross section to the SM prediction as a function of the hGG coupling shift κ h 1 arising in the EFT. As the Run 1 Higgs results restricts |κ h 1 | < 10%, we observe that an enhancement or suppression of the di-Higgs production rate of order 30% is still allowed by the data within the context of the EFT. In this case, one can easily understand the origin of the enhancement (suppression) when κ h 1 is negative (positive) by examining the interference between the box and triangle diagrams (see Fig. 1) via the function A(ŝ) in Eq. (3.4). For instance, when κ h 1 is negative, the smaller triangle amplitude is suppressed, while κ hh 1 = −κ h 1 is positive and the dominant box amplitude is enhanced. This implies that the interference between the amplitudes is reduced in comparison to the SM and the di-Higgs rate is enhanced.
There are other qualitatively distinct cases to consider. The first is when the new heavy colored states do not decouple from the Higgs as their mass is raised. This will occur if the new states obtain a substantial portion of their mass from electroweak symmetry breaking. In the language of the EFT, each operator in Eq. (1.2) is of similar size and thus the expansion is not useful from a practical point of view. This type of non-decoupling behavior is of course very familiar from the top quark contribution to the hGG and hhGG couplings. In this case it is instead necessary to specify the model for the new heavy colored states and apply the low energy theorems [54][55][56][57], as seen for light stops in Sec. 4.
Finally, the last case to consider is when the new states are light such that neither the EFT nor LET descriptions are valid. In the case of di-Higgs production, this occurs when the masses of the new states in the loop are similar to the characteristic invariant mass of the di-Higgs system under consideration. In this situation it is necessary to specify the model under consideration and compute the full one loop contribution to di-Higgs production. This is carried out for light stops in Sec. 5.
-7 -For the remainder of the paper we specialize to the case of stops in supersymmetry, which provides a well-motivated, concrete example of new colored particles with significant couplings to the Higgs. As is well known, the MSSM requires two Higgs doublets. However, motivated by the lack of evidence for new scalars and the fact that the Higgs production and decay rates are measured to be near the SM value, we take the 125 GeV Higgs to be the lightest neutral scalar boson and work in the decoupling limit. For the light stops that we consider in this work, we will typically not be able to obtain the 125 GeV Higgs mass in the MSSM. However, there are many possible scenarios that raise the Higgs mass including, for example, the NMSSM (for a review see [121,122]) or non-decoupling D-terms [123,124]. Therefore, we take the Higgs potential, particularly the triple Higgs coupling, to be that of the Standard Model in order to focus on the stop contributions.
We begin by describing our conventions for the stop sector. The stop mass matrix is given by where we have defined with∆ Q = 1 2 ( 1 2 g 2 − 1 6 g 2 ),∆ U = 1 2 ( 2 3 g 2 ). We also take v 2 u + v 2 d = v = 174 GeV and define tan β ≡ v u /v d . This matrix can be diagonalized, with eigenvalues m 1 and m 2 satisfying m 2 > m 1 , by performing a rotation of the basis by the angle θ defined by In this section we examine the generic corrections to the di-Higgs production rate in the limit that the stops are heavy in comparison to the typical di-Higgs invariant mass. As alluded to in the previous section, the stops can in general exhibit non-decoupling behavior as their masses are raised if the X t parameter is also raised in a correlated fashion. This is analogous to the case of the top quark in the SM. Because of this potential non-decoupling behavior, we we apply the Low Energy Theorem (LET) [54][55][56][57] to derive the couplings of the Higgs to gluons induced by stops. The starting point is the stop threshold contribution to the running of α s . After canonical normalization of the gluon field, we obtain the following effective Lagrangian: where b c 0 = 1 6 is the QCD beta function coefficient for stops.
Using Eq. (4.4) we determine the couplings of the Higgs h to gluons generated from stops by Taylor expanding around v u and v d in the Higgs fluctuations. Including the dominant SM top quark contribution, we arrive at the following effective Lagrangian describing the Higgs couplings to gluons: The coefficients κ h t , κ hh t (κ h t , κ hh t ) encode the top quark (stop) contributions to the hGG and h 2 GG couplings. In particular, for the stop contribution we have Here α is the mixing angle between the light and heavy CP-even Higgs bosons. Neglecting the small contributions from D-terms (g, g → 0) and taking the decoupling limit (α → β − π/2) we obtain where in the final step we have written κ hh t in terms of κ h t . These stop-induced contributions are to be compared to the top quark contributions, which in the decoupling limit are κ h t = κ hh t = 1. Therefore, the parameters κ h t and κ hh t measure the relative coupling shift from the SM values in an analogous way to the EFT coupling shifts defined in the previous section. We see from the last line in Eq. (4.8) that a definite correlation exists between the hhGG and the hGG couplings, and in the limit of heavy stops, m 1,2 m t , the hhGG coupling shift is fully determined by κ h t . As emphasized above, the current Run 1 data probe deviations in the hGG coupling at the 10% level, i.e., |κ h t | 10%. One can use this constraint to estimate the allowed size of the corrections to the di-Higgs rate from heavy stops by using Eq. (4.8). This is shown in Fig. 2, where we observe that O(50%) corrections are possible when the hGG coupling is smaller than its SM value by about 10%. The behavior can be easily understood by examining the couplings κ h t and κ hh t and accounting for the interference between the two diagrams depicted in Fig. 1. For instance, when κ h t is negative the s-channel Higgs exchange amplitude is slightly suppressed compared to its SM value, while the larger-inmagnitude contact diagram is instead mildly enhanced (since κ hh t is positive when κ h t is negative, assuming the stops are heavy). Therefore, the interference between the diagrams is less effective leading to the enhanced rate in this region, as shown in Fig. 2.
In Fig. 2 we can also see the importance of the non-decoupling behavior by comparing the EFT to the LET calculation. Because A-terms can cause the stops to get a large fraction of their mass from electroweak symmetry breaking even if they are relatively heavy, different and potentially larger effects in di-Higgs can be induced. Therefore, if a deviation is observed but no on-shell states are discovered, the size of the deviation could disentangle different types of decoupling vs non-decoupling new physics scenarios.

Light Stop Modifications: Full Loop Calculation
Finally, we consider the effects of light stops on the di-Higgs rate, which requires a full one loop analysis. To calculate the parton-level single Higgs and di-Higgs production cross sections we implemented the SM+Stops model described above into the FeynArts package [125,126] and employed the FeynArts, FormCalc, and LoopTools suite of packages [125,126] to calculate the amplitudes and evaluate loop functions. We used the MSTW [118][119][120] parton distribution functions when calculating the hadronic differential cross sections, with renormalization and factorization scales set to the invariant mass of the di-Higgs system. For the spectra in Fig. 3 we use constant K-factors to normalize our LO result to the NLO results in [79]. However, these K-factors cancel out in all other plots as only ratios of the BSM rate with the SM rate are shown. We have also cross checked our results using the full one-loop MSSM computations of Refs. [52,53], finding good agreement. 4 We begin by examining some benchmark models and their effect on the di-Higgs invariant mass spectra. In the SM, the amplitude for di-Higgs production vanishes at threshold because of a cancellation between the top box diagram and a triangle diagram that utilizes the triple Higgs coupling [5][6][7], and this is true for any field content as long all masses are acquired via the Higgs vacuum expectation value. Therefore, the invariant mass distribution in the SM is very small near threshold and grows to a peak near m hh ∼ 2m t , as we see in Fig. 3. Generic new physics that mediates one-loop di-Higgs production will spoil this cancellation, so light colored particles can lead to large deviations near threshold. We demonstrate this for some benchmark cases in Fig. 3.
Benchmark A has m 1 = 325 GeV, m 2 = 500 GeV, sin θ = 0.4: it has both stops light but the mixing angle is such that the the rate of gg → h is only enhanced by ∼ 15% and the h → γγ rate is within 5% of the SM value. This is a typical case where even having light stops the di-Higgs spectrum looks SM-like, and the total rate is ∼ 86% of the SM; a modification unobservable at the LHC. This also illustrates the effect found in Fig. 2 that 4 We differ in the writing of the function F3 defined in equation (B.2) of [53]: The total di-Higgs rate is enhanced by ∼ 70%, the single Higgs rate is reduced by ∼ 20%, and the di-photon modification is small. In benchmark C we show the generic but excluded case with one light stop: m 1 = 150 GeV, m 2 = 1000 GeV, sin θ = 0. Here the cancellation in the matrix element at threshold discussed in the introduction is spoiled and there is a large cross section enhancement at low invariant mass. The total cross section is enhanced by ∼ 90%, but the single Higgs rate is also enhanced by ∼ 80%.
We now discuss the expected modifications to the di-Higgs production rate as a function of more general stop sector parameters. Throughout we consider corrections to single Higgs and di-Higgs production. We will also consider two bins of di-Higgs invariant mass: 260 < m hh < 350 GeV and 260 < m hh < 2000 GeV. The first region is motivated because for light stops, the di-Higgs invariant mass spectrum can deviate significantly from the SM prediction for m hh < 2m t , as was illustrated in Fig. 3. Thus, although the total number of signal events may be smaller, when constraining new non-resonant contributions to di-Higgs production it may help to focus on di-Higgs invariant mass bins close to the threshold for production as this is where corrections are likely to be greatest. We also consider the full invariant mass regime to make contact with previous phenomenological studies that also do so.
In this section we consider corrections only at 14 TeV and provide contours for 100 TeV in App. A. The total di-Higgs production cross section increases substantially when going from 14 to 100 TeV, which is essentially due to the increased gluon luminosity. This is the main reason that sensitivity to di-Higgs production improves significantly with a 100 TeV proton-proton collider when compared to the LHC. However, for light stops the ratio of cross section modifications to the SM cross section remains roughly the same for both colliders. The reason for this is that although the total gluon luminosity in both cases is significantly different, the gradient of the gluon luminosity with respect to parton center of mass energy is not significantly different in the region of interest for di-Higgs production. Thus, when integrating over the parton distribution functions the increased gluon luminosity is roughly a constant factor, especially in the low invariant mass bin, and hence when the ratio of total cross section with stops to the total cross section in the SM is taken this factor essentially drops out. Therefore, the fractional corrections are very similar at 14 and 100 TeV. This does not persist whenever the stops are heavy and features in the invariant mass distribution appear at large m hh where the gluon luminosity between 14 and 100 TeV is significantly different, but in this case the corrections are typically smaller than the expected sensitivity. Thus, for the fractional corrections to the total cross section the 14 TeV results are also roughly illustrative of the 100 TeV result, although the expected sensitivity is increased at higher center of mass energy, so it should be kept in mind that contours of different di-Higgs cross section are appropriate in this case.
In general the stop parameter space can be described by three physical parameters, such as the two stop mass eigenvalues m 1 , m 2 , and the mixing angle, or alternatively the two soft masses m L , m R and the mixing parameter X t . To plot the corrections a projection down to a two-dimensional subspace is necessary. Results for a variety of projections for the full loop calculation are shown in Fig. 4, Fig. 5, and Fig. 6. In Fig. 4 the stop mixing X t -terms are set to zero and only the physical mass eigenvalues are varied. In Fig. 5 the two soft masses are set equal, m L = m R , and varied and the X t -term is also varied. The results are shown in the basis of physical masses. In Fig. 6 we fix the mass eigenvalue of the heavy stop to a benchmark value and then vary the light stop mass and the stop mixing angle.
In all of the figures a consistent picture emerges. Cases which lead to an observable deviation in the di-Higgs production rate also typically have observable deviations in the single Higgs production rate. Furthermore, in the 'blind spot' region where the single Higgs corrections are small the di-Higgs corrections are also typically suppressed unless both stops are quite light, which is consistent with our understanding based on the EFT arguments in Sec. 1. Thus, in order to indirectly constrain the existence of light stops which may have evaded direct detection at the LHC, the single Higgs production and di-Higgs production processes are highly complementary indirect probes and the strongest indirect constraints would arise from the combination of the two. Furthermore, taking into account current constraints, Fig. 4, Fig. 5, and Fig. 6 suggest that tuned regions of parameter space may remain after LHC8 in which observable non-resonant contributions to di-Higgs deviations may still arise at LHC14 from stops.   Fig. 4 with the exception that both stop soft masses are set equal and the Aterms are varied. In both cases regions which lead to a ∼ −20% change in the single production rate typically imply a ∼ 30% change in the pair production rate. The approximate color breaking vacuum constraint shown in blue is relevant for large mass splittings due to the large X t -terms.
It is also interesting that, as advertised previously, in Fig. 5 and Fig. 6 it is clear that deviations relative to the SM may be significantly larger in low invariant mass bins than they are for the total cross section. However, due to the smaller signal rate, the statistics   Fig. 4 with the exception that the heavy stop mass is fixed at 1000 GeV (upper panels) and 500 GeV (lower panels) and the light stop mass and mixing angle are varied. will be lower in the low mass bin than for the total cross section. Thus in a collider analysis aimed at indirectly constraining stop squarks it may be necessary to study a number of invariant mass cuts to determine the optimal constraint.
Finally, in Fig. 6 it is clear that if both stops are light the standard 'blind spot' in stop contributions to single Higgs productions may be closed by constraining di-Higgs production. This is consistent with our EFT discussion in Sec. 1 as even when the stop loop contributions to the hGG coupling have been tuned to precisely zero there will remain contributions to the h 2 GG coupling coming from a dimension-8 operator. Thus the h 2 GG coupling in the blind spot will typically be O(m 4 t /m 2 1 m 2 2 ). Hence, if we face the unfortunate situation that both stops are light and hGG deviations are absent due to a pernicious cancellation between stop loop contributions to the hGG coupling, it may still be possible to indirectly constrain this scenario through di-Higgs production measurements at the LHC.

Additional Indirect Constraints
As mentioned in Sec. 1, there are other indirect constraints on stops, and here we briefly comment on how those compare to the constraints and predictions considered here. The most relevant of these constraints comes from the observation that much of the parameter space considered has very large A-terms, and this can generate charge-or color-breaking vacua in the scalar potential that are deeper than the electroweak vacuum [127][128][129][130][131][132][133][134][135][136][137]. One can approximate the maximum allowed A-term by considering a D-flat direction in field space where H u = t L = t R , and requiring that all minima in that direction have positive vacuum energy. This leads to the condition [127][128][129][130][131] In the decoupling limit, m 2 Hu + |µ| 2 = −m 2 h /2 where m h is the physical Higgs mass. We can take the small µ or large tan β limit which sets A t = X t . This allows us to plot the bound from Eq. (5.2) in Fig. 4, Fig. 5, and Fig. 6.
We stress that Eq. (5.2) is a very crude approximation for the stability bound on the A-terms. In order to properly compute the bound, one must take into account loop contributions [138,139], tunneling effects [140][141][142], and properly account for cosmological history [143]. There now exist sophisticated computer codes [144] which can compute bounds in various different supersymmetric models [145][146][147]. Other groups have recently considered the effect of the 125 GeV Higgs on these bounds [148][149][150][151]. A full computation of the vacuum stability of our scenario is beyond the scope of this work, so we use Eq. (5.2) to give a rough sense of where those bounds would lie.
Precision electroweak observables can also be used to constrain this scenario [152]. One particularly important constraint comes from ρ-parameter which measures the splitting of electroweak multiplets. This constraint depends sensitively on the mass of the right handed sbottom as well as the mixing in the sbottom sector, and because of this additional model dependance we do not show the constraint on our figures. We find that generically the constraints from the ρ-parameter are weaker than the vacuum stability constraints.

Conclusions and Outlook
The observation and study of the di-Higgs channel is a primary objective of the LHC as well as future hadron colliders and it is a promising place to look for signatures of BSM physics. In this paper we have explored the impact of new colored states coupled to the Higgs particle on the production of Higgs boson pairs. Such states are well motivated by naturalness, with prime examples being top-partners. This class of non-resonant new physics can in principle lead to significant modifications to di-Higgs production. In most cases, however, the current experimental constraints on single Higgs production in the gluon fusion channel limit the extent to which the di-Higgs rate can deviate from the SM prediction. This can easily be seen in the case of heavy new colored states from an EFT analysis. The case of new light colored states requires a more detailed specification of the model and a full one loop calculation of the di-Higgs rate. We have performed such an analysis for the case of stops in supersymmetry, finding that that modifications are typically small, but that tuned regions with O(1) enhancements to the cross sections exist. This result demonstrates that future di-Higgs measurements could be used to place indirect constraints on the presence of light stops if they have somehow otherwise evaded detection at the LHC. However, these modifications are likely to be modest given the present constraints on single Higgs production. Thus, if large modifications in the di-Higgs production rate were observed this work would suggest that they are more likely to come from resonant new physics, or modifications of the weak sector and/or Higgs self-coupling, rather than from non-resonant contributions from new colored fields coupled to the Higgs.