Soft gluon resummation for Higgs boson pair production including finite Mt effects

We perform the all orders resummation of threshold enhanced contributions for the Higgs boson pair production cross section via gluon fusion, including finite top quark mass ($M_t$) effects. We present results for the total cross section and Higgs pair invariant mass ($M_{hh}$) distribution. We obtain results at next-to-leading logarithmic accuracy (NLL) which retain the full $M_t$ dependence, and are matched to the full next-to-leading order (NLO) prediction. Our NLL+NLO results represent the most advanced prediction with full $M_t$ dependence for this process, and produce an increase of about 4% in the total cross section with respect to the NLO result for LHC energies, and for a central scale $\mu_0 = M_{hh}/2$. We also consistently combine the full NLL with the next-to-next-to-leading logarithmically (NNLL) accurate resummation computed in the Born-improved large-$M_t$ limit, and match it to the next-to-next-to-leading order approximation of Ref. \cite{Grazzini:2018bsd}, so called NNLO_FTa. We find that the resummation effects are very small at NNLL for $\mu_0 = M_{hh}/2$, in particular below 1% at 13 TeV, indicating that the perturbative expansion is under control. In all cases the resummation effects are found to be substantially larger for the central scale $\mu_0 = M_{hh}$, resulting in a more stable cross section with respect to scale variations than the fixed order calculation.


Introduction
The study of the properties of the Higgs boson discovered by the ATLAS and CMS collaborations is one of the main goals of the present and future runs of the LHC. Among the different measurements that can help to distinguish between the Standard Model (SM) and new physics scenarios, the measurement of the Higgs self coupling is one of particular interest, as in the SM it is determined by the scalar potential, responsible for the electroweak symmetry breaking mechanism.
The production of Higgs boson pairs provides a direct way of measuring the Higgs trilinear coupling, and the high-luminosity upgrade of the LHC is expected to provide constraints on its value by measuring the double Higgs production cross section [2,3]. In the SM, the main production mechanism is the fusion of gluons via a heavy quark (mainly top quark) loop, and the corresponding cross section has been computed at leading order (LO) in Refs. [4][5][6]. The QCD corrections for this process have been computed first in the heavy top-quark mass (M t ) limit (HTL), both at next-to-leading order [7] (NLO) and next-to-next-to-leading order [8][9][10][11] (NNLO), and more recently the NLO corrections with full M t dependence also became available [12,13], later also supplemented by transverse momentum resummation [14] and parton shower effects [15,16]. The size of the QCD corrections was found to be large -about a 70% increase in the total cross section at NLO for LHC energies-, and also the difference with respect to the HTL was found to be significant, the latter being around 15% larger than the full NLO result at 14 TeV.
Very recently, an improved and fully differential NNLO prediction -labeled NNLO FTa for full-theory approximation, see also Refs. [17,18]-was presented in Ref. [1], which in particular features the full loop-induced double-real corrections. This result predicts an additional increase in the total cross section with respect to the full NLO calculation of about 12% at the LHC, and a residual uncertainty due to missing finite-M t effects estimated to be about 2.5%.
Besides the previously described fixed-order calculations, the all-orders resummation of soft gluon emissions has also been performed -again within the HTL-at next-to-next-to-leading logarithmic accuracy (NNLL) in Refs. [19,20]. The resummed contributions, which account for the dominant effect of the missing higher-orders in the perturbative expansion in the threshold limit, are found to further stabilize the cross section leading to smaller theoretical uncertainties.
In this work we perform the resummation of the threshold enhanced contributions including finite M t effects. In particular, up to next-to-leading logarithmic accuracy (NLL) we retain the full M t dependence, therefore obtaining NLL+NLO results that represent the most advanced prediction computed in the full theory. Finally, by performing matching to the NNLO FTa cross section, we achieve the state of the art results for Higgs pair production by reaching NNLL accuracy within the best available approximation for the M t effects .
This work is organized as follows: in section 2 we collect all the analytical expressions needed to perform threshold resummation up to NNLL, then in section 3 we present our numerical predictions for the LHC and future colliders, and in section 4 we summarize the results.

Threshold resummation
We consider the hadronic production of Higgs boson pairs via gluon fusion. The hadronic cross section for a collider center-of-mass energy s H , differential in the Higgs pair system invariant mass M hh , can be expressed in the following way where τ = M 2 hh /s H , µ R and µ F are the renormalization and factorization scales respectively, andσ 0 represents the Born level partonic cross section. The parton densities of the colliding hadrons are denoted by f a/h (x, µ 2 F ) with the subscripts a, b labeling the type of massless partons (a, b = g, q f ,q f , with N f = 5 different flavours of light quarks). The hard coefficient function G ab can be computed in perturbation theory, expanding it in terms of powers of the (MS renormalized) QCD coupling α S (µ 2 R ) as: We introduce now the notation needed to perform the soft gluon resummation in Mellin space [21,22]. We start by considering the Mellin transform of the hadronic cross section, which takes the following factorized form Here we have introduced the N -moments of the hard coefficient function and parton distributions, specifically Once all the ingredients in N -space are known, we can obtain the physical cross section via Mellin inversion, where the constant C M P defining the integration contour in the N -plane is on the right of all the possible singularities of the integrand [23].
We will perform the all-order summation of the threshold enhanced contributions, which corresponds to the limit z → 1 or equivalently N → ∞ in Mellin space, and appear as α n S ln m N terms with 1 ≤ m ≤ 2n. We will therefore consider (for the resummed contributions) only the gluon-initiated configuration, given that it is the only partonic channel that is not suppressed in this limit. The soft-gluon contributions in the large-N limit can be organized in the following all-order resummation formula for the partonic coefficient function in Mellin space, All the large logarithmic corrections are exponentiated in the Sudakov factor ∆ N , only depending on the dynamics of soft gluon emissions from the initial state partons. It can be expanded as The term ln N g (1) resums all the LL contributions α n S ln n+1 N , g (2) collects the NLL terms α n S ln n N , α S g (3) contains the NNLL terms α n+1 S ln n N , and so forth. The perturbative coefficients g (n) needed to perform NNLL resummation are known and only depend on the type of incoming partons, and their explicit expression can be found, for instance, in Refs. [24,25].
All the contributions that are constant in the large-N limit are contained in the function C gg (α S ). They originate in non-logarithmic soft contributions and hard virtual corrections, and can be expanded in powers of the strong coupling: In particular, in order to perform N i LL resummation we need up to the C (i) gg coefficient. At the same time, this coefficient can be obtained from the N i LO fixed order computation; even more, given that the soft gluon contributions in C (i) gg are universal, the only process dependence enters via the virtual corrections. The explicit (universal) relation between C (i) gg and the loop corrections has been derived up to i = 2 in Ref. [26], and later at one order higher in Ref. [27], and reads (for µ R = µ F = M hh ) where ζ n represents the Riemann zeta function, γ E is the Euler number and β 0 = (11C A − 2N f )/12π. The infrared-regulated one and two-loop correctionsσ (1) fin andσ (2) fin can be obtained from the corresponding matrix elements after applying the corresponding subtraction operator. The explicit formulas can be found in Ref. [26]. For the particular case of Higgs boson pair production, their explicit expression valid in the HTL can be found in Ref. [20], while for the NLL resummation with full M t dependence we can obtain numerical results forσ (1) fin , and therefore C (1) gg , using the publicly available grid interpolation of the two-loop NLO virtual corrections [15].
Finally, in order to fully profit from the knowledge of the fixed order calculation, we implement the corresponding matching. As usual, we expand the resummed N i LL cross section to O(α i s ) * , add the full N i LO cross section, and subtract the expanded result of the resummed one to avoid a double counting of logarithmic fixed order effects, as

Numerical results
In this section we present the numerical predictions for the LHC and future hadron colliders. We use the values M h = 125 GeV and M t = 173 GeV for the Higgs boson and top quark masses, and do not consider bottom quark contributions. We use the PDF4LHC15 sets [28][29][30][31][32][33] for the parton densities and strong coupling, evaluated at each corresponding perturbative order. The fixed order cross sections are obtained from the implementation of Ref. [1], which is based on the publicly available computational framework Matrix [34].
In the first place, we present in section 3.1 the NLL+NLO predictions. It is worth to point out that, even if more advanced predictions have been obtained for this process (specifically the so-called NNLO FTa defined in Ref. [1]), these results represent the most advanced prediction computed in the full theory, i.e. with full M t dependence. * Relative to the LO α 2 S power, which is always understood. Based on the knowledge of the threshold enhanced contributions at NLL with full M t dependence, and in particular on the O(α 2 S ) of its expansion, we can also provide an improved fixed order (approximated) NNLO prediction. This is presented in section 3.2. Finally, we combine the full NLL calculation with the NNLL contributions computed in the heavy top limit. This is presented in section 3.3.

NLL+NLO with full M t dependence
The results for the total cross section are shown in Table 1 for different center-of-mass energies. We use as the central scale µ 0 = M hh /2, though we also present results for µ 0 = M hh . Scale uncertainties are obtained via the usual 7-point variation.
We can observe that the size of the threshold effects goes down for larger collider energies, as expected from the fact that more energy is available and therefore soft gluon contributions become less dominant. As it was also observed in the heavy M t limit, we can appreciate that the size of the threshold corrections is much larger for µ 0 = M hh , ranging from 21.3% at 7 TeV to 10.2% at 100 TeV. The corresponding values for µ 0 = M hh /2 are 6.0% and 1.7%, respectively. For LHC energies, the soft gluon resummation effects are of the order of 4% for the central scale µ 0 = M hh /2. We can observe a reduction in the scale uncertainties (except for the 100 TeV predictions, where fixed-order and resummed results are comparable), this reduction being stronger for smaller center-of-mass energies. In fact, the NLL relative scale uncertainties remain practically unchanged when varying the collider energy, being always about ±10%.
In Table 2 we present the ratio of the central values for the predictions corresponding to µ 0 = M hh /2 and µ 0 = M hh , both for the fixed-order and resummed results. We can observe that the variation is substantially smaller in the resummed case, pointing towards a clear √ s  improvement in the stability of the cross section when taking into account the all-orders soft gluon effects.
We also present NLL predictions (with µ 0 = M hh /2) for the Higgs pair invariant mass M hh , at 7 TeV, 13 TeV (Figure 1), 27 TeV and 100 TeV (Figure 2). The lower plots show the ratio to the NLO result. We can see that the effect of the resummed contributions becomes larger as the invariant mass of the system increases, which again is expected due to the fact that less energy is available for extra emission. The increase in the Sudakov factor is however partially compensated by a suppression at large M hh in the NLO virtual corrections entering in C (1) gg , leading to a rather mild increase in the tail. Also here we can clearly observe that the resummation effects decrease with the collider energy.
It is interesting to compare our results with the ones obtained in the heavy-M t limit [20]. In order to do so, we present in Figure 3 the ratio between the NLL and NLO predictions as a function of M hh for different collider energies, both in the full theory and in the HTL. We can observe that there are clear differences in the shape, with the results with full M t dependence growing faster for lower invariant masses but showing a relative suppression with respect to the large-M t results in the tail. Still, this difference in the M hh spectrum between the two predictions is of the order of ±1%, and it is moderate compared to the overall effect of the resummed contributions. This indicates certain stability in the M t dependence of the threshold effects, and therefore the lack of full M t dependence at NNLL should lead to a rather small residual uncertainty due to missing finite-M t effects

Improved NNLO FTa
As it was mentioned in the previous section, the NLL+NLO results represent the most advanced prediction available for double Higgs production in the full theory. However, higher order corrections are still sizeable and therefore they need to be included in order to obtain accurate results, even if they are known only in an approximated way. The best fixed order prediction  available in the literature is the so-called NNLO FTa [1], which is obtained by working in the heavy M t limit but improved via a reweighting technique in order to account for finite M t effects. In particular, the NNLO FTa includes the full double-real loop induced squared matrix elements.
Before presenting combined NNLL+NNLO FTa predictions in the following section, it is worth to discuss possible improvements to the approximated NNLO result of Ref. [1] based on the knowledge of the full NLL+NLO result. Expanding the NLL+NLO results to O(α 2 S ) -where an overall α 2 S from the Born cross section is understood-, we can obtain the exact threshold enhanced contributions proportional to α 2 S ln 2 N † . Even if it features the full doublereal corrections, these contributions are obtained only within the (Born-improved) heavy M t limit in the NNLO FTa , because of the approximation performed in the real-virtual piece of the calculation. Therefore, we can define an improved NNLO FTa (denoted as NNLO FTa−i ) in the following way ‡ σ NNLO In Table 3 we show the comparison between the NNLO FTa and NNLO FTa−i predictions for the total cross section. We can observe that the difference is very small, being always below 0.5%. Even if this does not represent a proof of the accuracy of the NNLO FTa , the smallness of this effect points in this direction, and the difference is largely included within the estimated † Contributions proportional to α 2 S ln 3 N and α 2 S ln 4 N are already obtained in an exact way at LL, and are also reproduced with full M t dependence by the NNLO FTa . ‡ Besides having the full M t dependence in the α 2 S ln 2 N term, the NNLO FTa−i differs from the NNLO FTa result also in the term proportional to α 2 S ln N , though in this case the full M t dependence is only in those contributions generated by the NLL resummation.  M t uncertainty reported in Ref. [1].
In Figure 4 we present the Higgs boson pair invariant mass distribution for both NNLO approximations, for a collider energy of 13 TeV. We can observe that the difference between them is again very small in the whole invariant mass range, slowly growing with M hh but always within the scale uncertainties. This behavior is not surprising since the NNLO FTa is expected to be less accurate for large values of M hh , and also because the difference between NNLO FTa and NNLO FTa−i is only in threshold enhanced terms, which become more relevant for larger invariant masses. We can also observe that the scale uncertainties are larger for the NNLO FTa−i in the tail, being the central value corresponding to µ 0 = M hh /2 in the middle of the uncertainty band, while for the NNLO FTa it is located close to the upper limit. This fact reflects in the slightly larger scale uncertainties for the NNLO FTa−i total cross section that can also be observed in Table 3.
In summary, both for the total cross section and the invariant mass distribution we find that the differences between the NNLO FTa and NNLO FTa−i predictions are well within the estimated uncertainties inherent to these approximations.

NNLL resummation
We present now the NNLL predictions. In order to account for the NLL contributions with full M t dependence, we add the difference between the full theory and HTL predictions at NLL. Specifically, defining we have that our NNLL+NNLO FTa−i cross section is given by For the sake of brevity, we will denote this result NNLL FTa−i . Note that the NNLL result is matched to the NNLO FTa−i prediction instead of NNLO FTa , though as it was seen in the previous section the difference between the two is very small. In Table 3 we present the NNLL FTa−i predictions for the total cross section, for µ 0 = M hh /2. We can observe that the resummed contributions result in a small increase with respect to the NNLO FTa−i result, ranging from 1.3% at 7 TeV to 0.1% at 100 TeV, and being around 0.8% at the LHC. Again, the effect is much larger for the central scale µ 0 = M hh , where for instance the increase in the total cross section at 13 TeV is above 8%.
From Table 3 we can also compare the NNLL predictions with the NNLO FTa results of Ref. [1]. We can observe that the increase due to the resummed contributions is partially compensated with the existing decrease from the NNLO FTa to the NNLO FTa−i predictions, accidentally making the difference between the NNLO FTa and NNLL FTa−i results even smaller. The largest difference between these two predictions is in the scale uncertainties, which are comparable in size but turn out to be more symmetric for the NNLL FTa−i result.
In Table 4 we compare the fixed order NNLO FTa−i and resummed NNLL FTa−i predictions for the scale choices µ 0 = M hh /2 and µ 0 = M hh . In accordance with what was observed at NLO and NLL, we can see that the fixed order results present a larger variation in the central value when changing the renormalization and factorization scales, while the resummed results show a better stability. Again, this effect is less strong when we increase the collider energy.
Finally, in Figures 5 and 6 we present the Higgs pair invariant mass distribution at different collider energies. We can see again that the threshold effects increase with M hh by comparing the NNLO FTa−i and NNLL FTa−i curves. We observe that also at a differential level that the difference between the NNLO FTa and NNLL FTa−i predictions is very small, being below or  around 1% in the mass range under study. The difference in the scale uncertainty bands between these two predictions can also be appreciated, specially in the tail.
In conclusion, the difference between the resummed NNLL FTa−i prediction and the NNLO FTa result turns out to be small for µ 0 = M hh /2 compared to the size of the theoretical uncertainties, except only for the effect in the shape of the scale uncertainty bands. The small impact of the all orders soft gluon resummation is an indication of the good control over the perturbative expansion.

Summary
In this work we have computed the threshold resummation for Higgs boson pair production at hadron colliders via gluon fusion, including finite M t effects. We presented results both at NLL and NNLL accuracy, consistently matched to the corresponding fixed order cross sections.
Our NLL+NLO predictions retain the full M t dependence, and represent the most advanced prediction for this process computed in the full theory, i.e. not relying on the large-M t limit. We found that at 13 TeV the NLL+NLO cross section is larger than the NLO result by about 4.1% for the central scale µ 0 = M hh /2, while this effect goes up to 16.7% for µ 0 = M hh . The size of the resummed contributions decreases with the energy, going down to 2.8% and 1.7% at 27 and 100 TeV respectively, again for µ 0 = M hh /2. We observed clear differences in the shape of the corrections as a function of M hh with respect to the large-M t result, but moderate compared to the overall size of the threshold effects.
Using the knowledge of the full NLL contributions, we have defined an improved NNLO approximation, NNLO FTa−i . We found that the difference with respect to the NNLO FTa of Ref. [1] is very small, always below 0.5% for all the collider energies under consideration and well within the estimated M t uncertainties of the approximation, pointing towards the reliability  of the NNLO FTa result.
Finally, we have also consistently combined our full NLL predictions with the NNLL resummation computed in the large-M t limit, and matched it to the NNLO FTa−i result, thus providing a prediction for the Higgs boson pair production cross section with the most advanced ingredients available to date. We found that the effect of the resummed contributions is small at this order, being about 0.8% at the LHC and smaller for larger collider energies. The effect is again larger for µ 0 = M hh , being around 8.1% at 13 TeV. The small size of the threshold resummation effects at NNLL, specially for µ 0 = M hh /2, is an indication of the fact that the perturbative expansion is under good control, and that no sizeable higher order effects are expected beyond the order reached within this calculation.