NLO-QCD corrections to Higgs pair production in the MSSM

We take a step towards a complete NLO-QCD determination of the production of a pair of Higgs scalars in the MSSM. Exploiting a low-energy theorem that connects the Higgs-gluon interactions to the derivatives of the gluon self-energy, we obtain analytic results for the one- and two-loop squark contributions to Higgs pair production in the limit of vanishing external momenta. We find that the two-loop squark contributions can have non-negligible effects in MSSM scenarios with stop masses below the TeV scale. We also show how our results can be adapted to the case of Higgs pair production in the NMSSM.


Introduction
After the discovery of a Higgs boson in Run 1 of the LHC [1,2], one of the major goals of Run 2 is the experimental exploration of its properties. In Run 1, the couplings of the Higgs boson to fermions and to gauge bosons have already been measured, and found to be compatible with the predictions of the Standard Model (SM) within an experimental accuracy of (10 -20)% [3]. On the other hand, the self-couplings of the Higgs boson, which are accessible in multi-Higgs production processes, have not been probed yet. While a measurement of the quartic Higgs self-coupling lies beyond the reach of the LHC [4,5], previous studies showed that the Higgs pair production process, and hence the trilinear Higgs self-coupling, might be accessible for high integrated luminosities in the bbγγ [6][7][8][9][10][11], bbτ τ [7,12], bbW + W − [13] and bbbb [14][15][16] final states.
Not only is Higgs pair production interesting as a probe of the trilinear Higgs self-coupling in the SM, but it also can help constrain the SM extensions. First limits on scenarios with strongly increased cross section, which occurs, e.g., in models with novel hhtt coupling [17][18][19], or if the Higgs boson pair is produced through the decay of a heavy new resonance, have been given in refs. [20][21][22][23][24].
In the minimal supersymmetric extension of the SM (MSSM) the Higgs sector consists of two SU(2) doublets, H 1 and H 2 , whose relative contribution to electroweak (EW) symmetry breaking is determined by the ratio of vacuum expectation values (VEVs) of their neutral components, tan β ≡ v 2 /v 1 . The spectrum of physical Higgs bosons is richer than in the SM, consisting of two neutral scalars, h and H, one neutral pseudoscalar, A, and two charged scalars, H ± . The couplings of the scalars to matter fermions and gauge bosons, as well as their self-couplings, differ in general from the SM ones. However, in the so-called decoupling limit of the MSSM Higgs sector, m A m Z , the lightest scalar h has SM-like couplings and can be identified with the particle discovered at the LHC, with m h ≈ 125 GeV [25].
The dominant mechanism for Higgs pair production in the MSSM is gluon fusion, 1 mediated by loops involving the top and bottom quarks and their superpartners, the stop and sbottom squarks. Only for relatively light squarks, with masses below the TeV scale, do the squark contributions lead to sizeable effects on the cross section for the production of SM-like Higgs pairs [28]. Direct searches leave several corners of parameter space for light stops open, e.g. for reduced branching ratios or difficult kinematic configurations [29][30][31]. However, the measured value of m h implies either stop masses in the multi-TeV range or a large and somewhat tuned left-right mixing in the stop mass matrix. Scenarios allowing for light stops are thus restricted to the latter possibility.
Due to the extended Higgs spectrum of the MSSM, a pair of light scalars can also be produced resonantly through the s-channel exchange of a heavy scalar, leading to a sizeable increase in the cross section [32][33][34][35]. In addition, mixed scalar/pseudoscalar pairs, and pairs of pseudoscalars, can as well be produced in gluon fusion. In this paper, however, we will restrict our attention to the production of scalar pairs.
In the SM, the leading-order (LO) cross section for Higgs pair production via gluon fusion, fully known since the late eighties [36], is subject to large radiative corrections. The next-to-leading order (NLO) QCD contributions of diagrams involving top quarks were computed in the late nineties in the limit of infinite top mass m t , or, equivalently, of vanishing external momenta [33]. Whereas this approximation was shown to work quite well for single Higgs production [37], it can be expected to be less effective for pair production, due to the larger energy scale that characterizes the latter process. Unfortunately, an exact two-loop calculation of the "box" form factor that contributes to Higgs pair production at the NLO in QCD is currently not available. In contrast, the "triangle" form factor entering diagrams where a single (s-channel) Higgs boson splits into a Higgs pair can be borrowed from the the calculation of single Higgs production [37][38][39][40].
In order to improve the NLO result for Higgs pair production in the SM, ref. [33] factored out the LO cross section with the full top-mass dependence. The uncertainties of this approach were estimated to be of O(10 %) in refs. [41][42][43][44]. The next-to-next-to leading order (NNLO) contributions in the heavy-top limit were computed in refs. [45][46][47]. Soft gluon resummation at next-to-next-leading logarithmic (NNLL) order was performed in refs. [48,49]. Furthermore, NLO contributions in the heavy-top limit have been computed for the SM extended with dimension six operators [50], for an additional scalar singlet [51] and for the two-Higgs-doublet model [52].
In the case of the MSSM, the triangle form factor 2 that contributes to the production of a scalar pair at the NLO can again be borrowed from the calculation of single-scalar production. In particular, the contributions of two-loop diagrams involving only quarks and gluons can be adapted from the corresponding SM results [37][38][39][40] via a rescaling of the Higgs-quark couplings. The contributions of two-loop diagrams involving only squarks and gluons are fully known [39,40,53,54]. In contrast, an exact calculation of the two-loop diagrams involving quarks, squarks and gluinos -which can involve up to five different masses -is still missing. Calculations based on a combination of numerical and analytic methods were presented in refs. [55,56], but neither explicit formulae nor computer codes implementing the results of those calculations have been made available so far. Approximate results for the quarksquark-gluino contributions can however be obtained in the presence of some hierarchy between the relevant masses. The top-stop-gluino contributions were computed in the vanishing Higgs-mass limit (VHML) in refs. [57][58][59], and both the top-stop-gluino and bottom-sbottom-gluino contributions were computed in the limit of heavy superparticles -but without assuming a hierarchy between the Higgs mass and the quark mass -in refs. [60][61][62]. In particular, the calculation in ref. [59] relied on a low-energy theorem (LET) [63][64][65], connecting the amplitude for Higgs-gluon-gluon interaction to the derivatives of the gluon self-energy with respect to the Higgs fields, to provide explicit and compact analytic formulae for the top-stop-gluino contributions to the triangle form factor in the VHML. 2 In the MSSM, loop topologies other than triangle and box contribute to scalar pair production, due to the existence of quartic interactions involving squarks. With a slight abuse of language, in the following we denote as "triangle" all diagrams that involve the s-channel exchange of a single scalar, and as "box" all of the remaining diagrams.
For what concerns the box form factor, in the MSSM the contributions of one-loop diagrams involving quarks differ from their SM counterparts by a rescaling of the Higgs-quark couplings, and their calculation must be extended to account for the possibility of two different scalars in the final state [32]. The contributions of one-loop diagrams involving squarks have been computed in refs. [66,67] (see also ref. [28]). Going beyond the LO calculation, the contributions of two-loop diagrams involving top quarks and gluons in the heavy-top limit can be adapted from the corresponding SM results via a rescaling of the Higgs-top couplings [33]. On the other hand, the diagrams involving bottom quarks -whose effect is negligible in the SM, but can become relevant in the MSSM where at least one of the scalars has tan β-enhanced couplings to down-type quarks -are known only at one loop, because the heavy-quark limit adopted in the existing NLO calculations cannot, of course, be applied to them. Finally, no calculation of the contributions to the box form factor from two-loop diagrams involving squarks has, to our knowledge, been presented so far.
In this paper we take a step towards a complete NLO-QCD determination of the production of a pair of Higgs scalars in the MSSM. Relying on the same LET as in ref. [59], we obtain analytic results for the contributions to the box form factor from one-and two-loop diagrams involving top quarks and stop squarks in the limit of vanishing external momenta. We also obtain, by direct calculation of the relevant two-loop diagrams, the subset of bottom/sbottom contributions that involve the D-terminduced EW Higgs-squark coupling and survive in the limit of vanishing bottom mass. To assess the importance of the newly-computed corrections, we include the squark contributions to both triangle and box form factors in a private version of the public code HPAIR [68], which computes the NLO-QCD cross section for Higgs pair production in the SM and in the MSSM. We find that the two-loop squark contributions can have a non-negligible effect in scenarios with stop masses below the TeV scale. We conclude by discussing the limitations of the approximation of vanishing external momenta. Finally, in the appendices we collect some analytic formulae for the two-loop box form factors, and we show how our results can be adapted to the case of Higgs pair production in the next-to-minimal supersymmetric extension of the SM (NMSSM).

Higgs pair production via gluon fusion at NLO in the MSSM
In this section we summarize some general results on the gluon-fusion production of a pair of neutral Higgs scalars, denoted as φ and χ (each of them can be either h or H). The hadronic cross section for the process h 1 + h 2 → φ + χ + X at center-of-mass energy √ s can be written as where: M 2 φχ is the invariant mass of the φ + χ system; f a,h i (x, µ F ) is the density for the parton of type a (with a = g, q, q) in the colliding hadron h i ; µ F is the factorization scale;ŝ = s x 1 x 2 is the partonic center-of-mass energy;σ ab is the cross section for the partonic subprocess ab → φ + χ + X. The partonic cross section can be written in terms of the LO contribution σ The LO cross section is where: G F is the Fermi constant; α s (µ R ) is the strong gauge coupling expressed in the MS renormalization scheme at the scale µ R ; the Mandelstam variables of the partonic process,t and (for later convenience)û, are defined aŝ with θ the scattering angle in the partonic center-of-mass system, and The integration limits in eq. (3) are given bŷ corresponding to cos θ = ±1. Finally, in eq. (3) F φχ, 1 and G φχ, 1 represent the one-loop parts of the spin-zero and spin-two form factors for the process gg → φχ, respectively. While the spin-two form factor G φχ receives only contributions from box diagrams, the spin-zero form factor F φχ can be decomposed in box and triangle contributions as: In particular, F φχ contains the spin-zero part of the box diagrams, while F h ∆ (F H ∆ ) contains the contribution of the triangle diagrams for the production of an off-shell scalar h (H) which subsequently decays into the pair φχ through the factor C hφχ ∆ (C Hφχ ∆ ), defined as where λ hφχ is the trilinear scalar coupling 3 and Γ h is the width of the scalar h (in turn, C Hφχ ∆ is obtained from eq. (9) with the replacement h → H). The form factor F φ ∆ is decomposed in one-and two-loop parts as and analogous decompositions hold for F φχ , F φχ and G φχ . 3 We normalize all trilinear Higgs couplings to λ0 = m 2 Z /v, with v = ( The coefficient function G ab (z) in eq. (2) can in turn be decomposed, up to NLO terms, as with the LO contribution given only by the gluon-fusion channel: The NLO terms include, besides the gg channel, also the one-loop induced processes gq → qφχ and qq → gφχ. The gg-channel contribution, involving two-loop virtual corrections to gg → φχ and one-loop real corrections from gg → φχg, can be written as where In eq. (13), C A = N c (N c being the number of colors), β 0 = (11 C A − 2 N f )/6 (N f being the number of active flavors) is the one-loop β-function of the strong coupling in the SM, P gg is the LO Altarelli-Parisi splitting function and The first line of eq. (13) displays the two-loop virtual contribution regularized by the infrared singular part of the real-emission cross section. The second line contains the non-singular contribution from the real gluon emission in the gluon-fusion process. The function R gg is obtained from one-loop diagrams where only quarks or squarks circulate into the loop, and in the limit of vanishing external momenta it becomes R gg → −11(1 − z) 3 /(6z). The form factors F φχ ∆∆ and G φχ ∆∆ in eq. (14) represent the contributions of two-loop double-triangle diagrams with t/u-channel gluon exchange. In the limit of vanishing external momenta, the double-triangle form factors can be expressed in terms of the one-loop triangle form factors: with Finally, the contributions of the gq → qφχ and qq → gφχ channels are given by: where with The functions R qq and R qg in (19) are obtained from one-loop quark and squark diagrams, and in the limit of vanishing external momenta become

Box form factors in the limit of vanishing external momenta
As mentioned in section 1, exact results for the one-loop form factors F φχ, 1 and G φχ, 1 which determine the cross section for Higgs pair production at the LO have been known for a long time, both for the SM [36] and for the MSSM [32,66,67]. At two loops, the triangle contributions to the form factors can be borrowed from the calculation of the cross section for single Higgs production. However, explicit formulae for the contributions of triangle diagrams involving quarks, squarks and gluinos are available only in approximate form, assuming the existence of some hierarchy among the relevant masses and momenta [57][58][59][60][61][62]. Two-loop results for the box contributions to the form factors are known only for the diagrams involving top quarks and gluons, and only in the heavy-top limit [33].
In this section we present a novel calculation of the contributions of diagrams involving top quarks and stop squarks to the box component F φχ of the spin-zero form factor F φχ , up to the two-loop order. We restrict our calculation to the limit of vanishing external momenta, which, for the topgluon contribution alone, corresponds to the heavy-top limit. Note that the corresponding triangle component F φ ∆ can be extracted from ref. [59], and that the spin-two form factor G φχ vanishes in the zero-momentum limit. We also present results for the contributions of the diagrams involving sbottom squarks, under the additional approximation of vanishing bottom mass. Finally, we show how the formulae for the two-loop part of the form factors are affected by a change in the renormalization scheme of the parameters entering the one-loop part.
It is convenient to decompose the triangle and box form factors for the production of scalar mass eigenstates as where T F = 1/2 is a color factor (we make it explicit to follow the notation of ref. [59]), the angle α relates the scalar mass eigenstates, h and H, to the real parts of the neutral components of the two MSSM Higgs doublets, S 1 and S 2 , and H i and H ij , with i, j = (1, 2), are form factors in the interaction basis. As mentioned above, the form factors H i were computed in refs. [59,60,62] for single Higgs production. Finally, we further decompose the form factors H ij into top/stop and bottom/sbottom contributions,

Top/stop contributions via the low-energy theorem
In our derivation of the top/stop contributions to the box form factors we rely on the same LET for Higgs interactions [63][64][65] that was employed in ref. [59] for the calculation of the top/stop contribution to the triangle form factors. In our case, the LET connects the form factor for the interactions of two gluons with two Higgs scalars at vanishing external momenta to the second derivatives of the gluon self-energy with respect to the Higgs scalars. In particular, we can write the top/stop contributions to the form factors in the interaction basis as where Π t (q 2 ) denotes the top/stop contribution to the transverse part of the dimensionless (i.e., divided by q 2 ) self-energy of the gluon. In analogy with the effective-potential calculation of the MSSM Higgs masses in ref. [69] and with the LET calculation of single Higgs production in ref. [59], the dependence of the gluon self-energy on the Higgs fields S i can be identified through the field dependence of the top mass m t , the stop masses m 2 t 1 and m 2 t 2 and the stop mixing angle θ t , defined as A lengthy but straightforward application of the chain rule for the derivatives allows us to express the form factors as where A t is the trilinear soft-SUSY breaking Higgs-stop coupling, µ is the Higgs/higgsino mass term in the superpotential (with the sign convention of refs. [59,69]), and we define s 2θt ≡ sin 2θ t and, for later convenience, c 2θt ≡ cos 2θ t . We note that the first line of each equation and those in the second lines read where θ W being the Weinberg angle.
In appendix A we provide the explicit definitions of the two-loop functions F 2 i , F 2 , G 2 , F 2 i and D 2 in terms of the derivatives of the gluon self-energy. For the latter, we define the shortcut Z ≡ (2/T F ) Π 2 , t (0), after decomposing the gluon self-energy in one-and two-loop parts as Analytic formulae for the first derivatives of Z, computed under the assumption that the one-loop part of the gluon self-energy is expressed in terms of DR-renormalized top/stop parameters, were given in ref. [59]. Indeed, the functions F , G and D entering eqs. (29)-(31) coincide with those defined in that paper for the case of single Higgs production. Analytic formulae for the second derivatives of Z, which enter the functions F i and F i , can be easily obtained from those for the first derivatives, using the recursive relations for the derivatives of the two-loop function Φ(m 2 1 , m 2 2 , m 2 3 ) given e.g. in appendix A of ref. [70]. However, those formulae are too lengthy to be given explicitly in print, thus we make our results available upon request as a fortran routine.

Bottom/sbottom contributions for vanishing bottom mass
The LET employed in the previous section to compute the top/stop contributions to the box form factors relies on the assumption that the external momenta are negligible with respect to the masses of all particles running in the loops. Obviously, this assumption cannot hold for the contributions involving bottom quarks, nor for those involving quarks of the first two generations. In ref. [60] the bottom/sbottom contributions to single Higgs production were computed with an asymptotic expansion in the heavy supersymmetric masses (  (31). In particular, we find The one-loop parts of the functions F 3 b and D b read, in this approximation, where We obtained the two-loop parts of the functions F 3 b and D b by explicit computation of the relevant two-loop diagrams, setting m b = θ b = 0 from the start and taking the first non-vanishing term of an asymptotic expansion in the heavy superparticle masses (for an outline of this approach, see section 3 of ref. [60]). Under the assumption that the one-loop parts of the form factors are expressed in terms of DR-renormalized sbottom masses at the scale Q, we get with x L,R = m 2 b L,R /m 2 g and the notation (L → R) means a term that is obtained from the previous one with the exchanges x L → x R and d b L → d b R . We find that, when m b = θ b = 0, there are no infrareddivergent parts in the two-loop bottom/sbottom diagrams, therefore our results could also be obtained as the first non-vanishing term of a Taylor expansion of those diagrams in the external momenta. On the other hand, we stress that our results cannot be obtained by setting m t = θ t = 0 in the LET results for the top/stop contributions, because the latter rely on the assumption that the external momenta are much smaller than the quark mass. Finally, the contributions of the first two generations of quarks and squarks can be obtained, by means of trivial substitutions, from eqs. (29)-(31) and from the results presented in this section.

Change of renormalization scheme
The results presented in sections 3.1 and 3.2 were obtained under the assumption that the parameters entering the one-loop part of the form factors are expressed in the DR renormalization scheme. If a different scheme is used, the two-loop part of the form factor receives a shift where δH ij is a function of the shifts of all the parameters in the one-loop part of the form factor that are subject to O(α s ) corrections. 5 In the top/stop sector, the parameters that need shifting are the top mass, the stop masses, the stop mixing angle and the trilinear coupling A t . In particular, the shifts of those parameters to the on-shell (OS) scheme adopted in our numerical discussion can be found in appendix B of ref. [69]. The shifts δH t ij can then be written as where the one-loop parts of the functions F 2 , F 3 , F and F 2 are given in eqs. where If the sbottom masses in the one-loop part of the form factors are expressed in the OS scheme, the

The effect of SUSY contributions to Higgs pair production
In this section we present numerical results for the newly-computed SUSY contributions to the box form factors, and for their effect on the Higgs-production cross section. We focus on the process that is most interesting from the point of view of LHC phenomenology, i.e. the production of a pair of light MSSM scalars hh with mass m h ≈ 125 GeV.

Implementation in HPAIR
For the numerical evaluation of the cross section, we added the contributions of loops involving superparticles to the code HPAIR [68], whose public version includes by default the one-loop top-and bottom-quark contributions with full mass dependence [32] and the QCD corrections to the top-quark contributions in the heavy-top limit [33]. For the LO cross section, we added the one-loop squark contributions to the spin-zero and spin-two form factors, borrowing from ref. [66] the results with full mass dependence. At NLO, we included our results for the two-loop stop and (partial) sbottom contributions in the approximation of vanishing external momenta, derived in section 3. In order to improve on that approximation, the LO cross section factored out of the coefficient function G ab (z) in eq. (2) is computed with full dependence on the top and bottom quark and squark masses. In analogy with the implementation of the top quark loops in HPAIR, the gg-channel contribution to the NLO coefficient function in eqs. (13) and (14) specialized to the production of a hh pair -becomes where the subscript "LET" denotes form factors computed in the limit of vanishing external momenta after setting m b = θ b = 0. The two-loop SUSY contributions enter the last term in the first line of eq. (56), which, if only the top-quark contributions were considered as in ref. [33], would reduce to a simple coefficient c 1 = 11/2. The second line of eq. (56) contains the contributions of diagrams with t/u-channel gluon exchange. Following ref. [33], in those contributions we retain the full momentum dependence in the one-loop form factors that stem from the LO matrix element, but take the limit of vanishing external momenta, see eq. (17), in the double-triangle form factors. We also remark that in the NLO coefficient of eq. (56) all form factors -including those with full momentum dependenceare obtained for m b = θ b = 0. For a precise prediction of the cross section for the production of a pair of MSSM Higgs bosons, it is essential to include the corrections to the trilinear Higgs couplings, which can be as significant as the corresponding corrections to the MSSM Higgs masses and mixing. Indeed, to properly reproduce the decoupling limit in which the lightest scalar h has a SM-like self-coupling, λ SM hhh = 3 m 2 h /m 2 Z , the corrections to the coupling should be computed at the same level of accuracy as the corrections to the mass m h . The trilinear couplings are known at one loop [71][72][73][74][75], but at two loops only the O(α s α t ) corrections have been computed, in the effective-potential approximation, for both the MSSM [76] and the NMSSM [77]. In contrast, in this analysis we compute the MSSM Higgs masses and mixing using the code FeynHiggs [78][79][80][81][82], which includes two-loop corrections beyond the O(α s α t ) effectivepotential ones. Since we are anyway focusing on the effects of the SUSY contributions to the gluonfusion loop, we bypass the calculation of the corrections to the trilinear couplings by relying on a simplifying approach, known as "hMSSM", which was recently proposed in refs. [83][84][85][86]. In this approximation one assumes that the corrections to the elements other than (2, 2) of the Higgs mass matrix are negligible, i.e. ∆M 2 1j ≈ 0 with j = 1, 2. In that case the remaining correction ∆M 2 22 , which includes potentially large logarithmic effects from top/stop loops, can be expressed in terms of the parameters that determine the tree-level mass matrix (i.e. tan β, m Z and the pseudoscalar mass m A ) plus the lightest eigenvalue m h , treated as an input parameter: In this approximation the trilinear couplings relevant to the production of an hh pair become λ Hhh = 2 sin 2α sin (α + β) − cos 2α cos (α + β) Combining eqs. (57) and (58) one can see that in the decoupling limit m A m Z , when α → β − π/2, the coupling λ hhh does indeed tend to its SM limit. As discussed e.g. in refs. [87,88], the approximation of neglecting the corrections ∆M 2 1j might not prove accurate for small m A and rather large µ and tan β. We will therefore avoid those choices of parameters in our numerical example.

A numerical example
The SM parameters entering our computation of the cross section for Higgs pair production are the Z boson mass m Z = 91.1876 GeV, the W boson mass m W = 80.398 GeV, the Fermi constant G F = 1.16637 · 10 −5 GeV −2 and the pole top-quark mass m t = 173.2 GeV. We use the MSTW08 set of parton distribution functions [89][90][91] and the associated LO and NLO values of the strong coupling α s . The hadronic center-of-mass energy is set to √ s = 14 TeV. The factorization and renormalization scales are set to the invariant mass M hh of the Higgs boson pair.
We use the code FeynHiggs [78][79][80][81][82] to compute the masses and mixing angle of the Higgs scalars, taking as input the SM parameters listed above plus α s (m Z ) = 0.119. We consider an MSSM scenario characterized by the following parameters in the OS renormalization scheme: tan β = 10, m A = 500 GeV, µ = −400 GeV, M 3 = 1500 GeV, where M 3 denotes the soft SUSY-breaking gluino mass, we define X t ≡ A t + µ cot β, and mt L , mt R , and mb R denote the soft SUSY-breaking masses of the third-generation squarks. We recall that, in the OS scheme, the soft SUSY-breaking parameters in the squark sector are defined as the parameters entering a tree-level mass matrix that is diagonalized by the OS mixing angle (the latter defined, e.g., in appendix B of ref. [69]) and has the pole squark masses as eigenvalues. In this scheme the parameter mb L differs from its stop counterpart mt L by a finite shift [92,93]. The parameters in eq. (60) -as well as the remaining soft SUSY-breaking parameters, which are not relevant to our discussion -were chosen in such a way that, for M S = 500 GeV, they reproduce the light-stop benchmark scenario proposed in ref. [94] and studied in the context of single-Higgs production in ref. [95]. Our choices of m A and tan β ensure that the lightest Higgs scalar h has essentially SM-like couplings to the top and bottom quarks, and that the contribution of triangle diagrams with s-channel exchange of the heaviest scalar H is significantly suppressed, allowing us to focus on the effects of the SUSY contributions to the box form factor. We then vary the squark mass parameter M S between 500 GeV and 1500 GeV, which results in a lightest stop mass mt 1 ranging between 324 GeV and 1326 GeV, and in a prediction by FeynHiggs for m h ranging between 122.3 GeV and 130.7 GeV.
In figure 1 we plot the box form factor F hh -computed in the vanishing-momentum limit as described in section 3 -as a function of the squark-mass scale M S . The solid lines correspond to the one-loop (dark blue) and two-loop (light blue) part of the form factor, including both the top-quark contribution and the squark contributions. The dashed lines correspond to the one-and two-loop form factors including only the top contributions. The plot shows that the squark contributions can be relevant for small squark masses, and they are significantly larger in the two-loop form factor than in the one-loop form factor. Moreover, the decoupling behavior of the squark contributions for large M S appears to be slower at two loops than at one loop. This can be explained by the occurrence of two-loop terms proportional to m 2 t /M 2 ln(M 2 /m 2 t ) (with M denoting generically a SUSY mass parameter), whereas at one loop all terms decouple at least as fast as m 2 t /M 2 . In figure 2 we plot the cross section for the production of a hh pair as a function of M S , computed as described in section 4.1. The dark-blue lines correspond to the LO cross section, the light-blue lines to the NLO cross section, and again the solid (dashed) lines correspond to form factors including (not including) the SUSY contributions. 6 In addition, the dotted light-blue line corresponds to the NLO cross section computed by omitting the SUSY contributions in the two-loop part of the box form factor. The plot shows that, for the considered choices of MSSM parameters, the squark loops can significantly contribute to the cross section for relatively small M S , although their effect gets quickly suppressed when M S > ∼ 1 TeV. In particular, in the light-stop scenario -corresponding to the left edge of the plot -for our choices of m A and tan β the SUSY contributions increase the NLO cross section for h pair production by more than 30% (in contrast, ref. [95] showed that they reduce the cross section for the production of a single SM-like scalar by about 20%). Finally, the comparison between the solid and dotted light-blue lines shows that the newly-computed two-loop SUSY contributions to the box form factor account for a non-negligible part of the increase in the pair-production cross section.

Discussion
Relying on a low-energy theorem that connects the Higgs-gluon interactions to the derivatives of the gluon self-energy, we obtained analytic results for the contributions to Higgs pair production from oneand two-loop box diagrams involving top quarks and stop squarks in the limit of vanishing external momenta. We also obtained, by direct calculation of the relevant two-loop diagrams, the subset of bottom/sbottom contributions that involve the D-term-induced EW Higgs-squark coupling and survive in the limit of vanishing bottom mass. Combined with the existing results for the triangle diagrams in the same approximations [59,60], our calculation allows for a consistent NLO determination of the SUSY contributions to Higgs pair production in the MSSM. We incorporated our results in a private version of the code HPAIR, and found that the two-loop SUSY contributions to the production of a light-scalar pair can have a non-negligible effect in scenarios with stop masses below the TeV scale.
To conclude, a discussion is in order of the approximation of vanishing external momenta that we employed in our calculation. Our results can be viewed as the first term of an asymptotic expansion of the form factor F φχ, 2 in the heavy masses of all particles running in the loops. Such expansion is in principle valid only for partonic center-of-mass energies up to the lowest threshold encountered in the relevant diagrams, which for the contributions considered in this paper corresponds to √ŝ = 2 m t .
In the SM, the vanishing-momentum approximation is known to work rather well for the top-quark contributions to the production of a single scalar h with m h ≈ 125 GeV, because the region in the partonic phase space with √ŝ > 2 m t gives only a small contribution to the hadronic cross section. In contrast, the same approximation is less reliable for pair production, where it is always √ŝ > 2 m h and the whole region up to √ŝ ∼ 600 GeV gives a significant contribution to the cross section [41]. The factorization of the LO cross section with full momentum dependence is expected to reduce the uncertainty of the NLO result due to the dominance of soft and collinear gluon effects [33]. Nevertheless, a NLO determination of the top-quark contributions to Higgs pair production going beyond the vanishing-momentum -or, equivalently, infinite-top-mass -approximation would be desirable. Of the necessary ingredients, the contribution to F h, 2 ∆ of two-loop triangle diagrams involving top quarks and gluons is known with full top-mass dependence from single-Higgs production; the contribution of oneloop top diagrams to R gg , R qq and R qg is known exactly from ref. [44]; the contribution of two-loop, one-particle-reducible top diagrams to F φχ ∆∆ and G φχ ∆∆ is relatively easy to compute. However, an exact evaluation of the two-loop box diagrams involving top quarks and gluons is currently not available, and represents the bottleneck in the quest for an exact NLO determination of the pair-production cross section. Attempts to go beyond the limit of infinite top mass for the two-loop box diagrams were made in refs. [41,42], where several terms in a heavy-top asymptotic expansion of the cross section, i.e. terms proportional to powers ofŝ/m 2 t or m 2 h /m 2 t , were obtained. However, as shown explicitly for the LO result in refs. [96,97], the inclusion of additional terms in the large-mass expansion does not necessarily improve the evaluation of the cross section. Indeed, by including additional terms one is improving the evaluation of the region with √ŝ < 2 m t at the price of worsening the evaluation of the complementary region with √ŝ > 2 m t , which is approximated by a function that has the wrong behavior asŝ increases. In fact, the appropriate expansion in the region with √ŝ > 2 m t would be a large-momentum expansion as opposed to a large-mass expansion.
In the MSSM, the NLO cross section for the production of a pair of SM-like scalars hh suffers from the same uncertainty as in the SM, stemming from the incomplete knowledge of the two-loop diagrams with top quarks and gluons. For what concerns the SUSY contributions, those from two-loop diagrams involving squarks and gluons or quartic squark couplings should be sufficiently well approximated, in realistic MSSM scenarios, by the results obtained in the vanishing-momentum limit. In contrast, some two-loop diagrams involving top, stop and gluino do have thresholds at √ŝ = 2 m t , thus their contributions are in principle subject to uncertainties comparable to those of the SM contributions. The knowledge of those contributions could however be improved following the same strategy employed in ref. [62] for single scalar production, namely evaluating the top-stop-gluino box diagrams via a largemass expansion in the SUSY masses while treating the top quark as a light particle.
Finally, another feature specific to the MSSM calculation of hh production is the possibility of large resonant contributions from triangle diagrams with s-channel exchange of the heaviest scalar H. In such a scenario, the determination of the NLO cross section could be improved by using for F H, 2 ∆ the quark-gluon contributions with full momentum dependence combined with the heavy-SUSY results of refs. [60,62], while retaining the vanishing-momentum approximation in F h, 2 ∆ to avoid spoiling potential cancellations with the box form factor.
The functions that represent the sub-dominant contributions of diagrams involving D-term induced EW couplings read where δG = 1 6 and

Appendix C Extension to the NMSSM
In this appendix we describe how our results for the box form factor for Higgs pair production in the MSSM can be extended to the case of the NMSSM. Instead of the Higgs mass term µ H 1 H 2 , which in the simplest realization of the NMSSM is forbidden by a Z 3 symmetry, the superpotential contains 7 where S is an additional gauge-singlet superfield. An effective µ term is generated by the singlet VEV as µ = λ S , and the CP-even parts S i of the neutral component of the three Higgs fields -ordered as (H 1 , H 2 , S) -mix into three mass eigenstates which we denote as h a , where R S is an orthogonal matrix. The decompositions of the triangle and box form factors in eqs. (21)-(25) generalize to The extension to the NMSSM of the results of refs. [59,60,62] for the triangle form factors of the MSSM has been presented, in the context of single Higgs production, in ref. [99]. Concerning the box form factors, the terms H 11 , H 12 and H 22 coincide with those obtained for the MSSM in section 3. The top/stop contributions to the remaining terms read where the functions F 2 , F 3 , F and F 2 coincide with those entering the MSSM results, see section 3 and appendix A. In the limit m b = θ b = 0 there are no contributions to H 13 , H 23 and H 33 from bottom/sbottom loops. 7 For consistency with the definition of µ in our MSSM results, here we adopt for the sign of λ the opposite convention with respect to ref. [98] and most public codes for NMSSM calculations. We also note that our normalization of the EW parameter, v ≈ 246 GeV, differs by a factor √ 2 from the one in ref. [98].
Finally, when the parameters entering the top/stop contributions to the one-loop part of the form factors are expressed in a renormalization scheme other than DR, the shifts to the form factors that were not already given in section 3.3 read where the shifts δF 2 , δF 3 , δF and δ F 2 coincide with those defined in appendix B.