Transverse-momentum resummation for Higgs boson pair production at the LHC with top-quark mass effects

We consider Higgs boson pair production via gluon fusion in hadronic collisions. We report the calculation of the transverse-momentum ($q_T$) distribution of the Higgs boson pair with top-quark mass ($M_t$) effects fully taken into account. At small values of $q_T$ we resum the logarithmically-enhanced perturbative QCD contributions up to next-to-leading logarithmic (NLL) accuracy. At intermediate and large values of $q_T$ we consistently combine resummation with the ${\cal O}(\alpha_S^3)$ fixed-order results. After integration over $q_T$, we recover the next-to-leading order (NLO) result for the inclusive cross section with full dependence on $M_t$. We present illustrative numerical results at LHC energies, together with an estimate of the corresponding perturbative uncertainties, and we study the impact of the top-quark mass effects.


Introduction
Since the discovery of the Higgs boson (H) [1,2] during Run I of the Large Hadron Collider (LHC), it became a physics goal for Run II and future high-energy collider facilities to complete our understanding of the electroweak symmetry breaking mechanism of the Standard Model (SM). The study with increased precision of the couplings between the Higgs boson and the SM particles shows a determination of the couplings to vector bosons and heavy fermions compatible with the SM values with an overall 15% to 20% uncertainty [3,4,5]. The results obtained in Refs. [3,4,5] concern single Higgs boson production.
In order to probe the Higgs boson self couplings one can consider the process of double Higgs boson (HH) production [6,7,8,9,10,11]. In this case, similarly to single Higgs boson production, the main production mechanism is driven by gluon fusion. At the leading order (LO), the two Higgs bosons can couple to a heavy-quark loop via a box diagram or, via the trilinear Higgs selfcoupling, to an off-shell Higgs boson produced by a triangular heavy-quark loop. For this reason, the observation of Higgs boson pair production gives access to a direct extraction of the Higgs trilinear self-coupling and to the reconstruction of the Higgs boson potential [12,13,14,15,16,17].
Predictions for double Higgs boson production in gluon fusion at LO were obtained in Refs. [18,19,20] including full top-quark mass (M t ) effects. However, since the gluon fusion mechanism for HH production is a loop induced process, the next-to-leading order (NLO) QCD corrections were first obtained in the heavy top quark limit M t → ∞ [21] using the Higgs effective field theory (HEFT). In this approximation the top-quark mass is regarded much larger than any other scale in the process and the top quark is integrated out at the Lagrangian level. This significantly simplifies the calculation of the NLO corrections since the top-quark loops shrink to a pointlike interaction of the Higgs bosons with gluons. More recently, next-to-next-to-leading order (NNLO) predictions for HH production in the HEFT have been completed in Refs. [22,23,24,25]. Threshold resummation up to next-to-next-to-leading logarithmic (NNLL) accuracy in the HEFT has been performed in Refs. [26] and [27] matching the resummed results respectively with NLO and NNLO fixed-order calculations.
However, Higgs boson pairs are produced with an invariant mass (M HH ) which is above the topquark mass threshold where the validity of the HEFT description breaks down. Therefore, various approximations to include finite M t effects beyond LO have been performed in the literature. In the so-called "Born-improved HEFT" approximation a reweighting of the NLO HEFT result is performed using a factor B F T /B HEF T , where B F T and B HEF T denote the LO matrix element squared in the full theory and in the HEFT respectively [21]. In the NLO calculation in Refs. [28,29] the top-quark mass dependence is fully taken into account in the real emission correction, while the virtual amplitude is computed in the heavy top quark limit and reweighted by the Bornimproved factor. HEFT results at NLO and NNLO improved by an expansion in 1/M 2 t have been obtained in Refs. [24,30,31,32].
The NLO calculation including the full top-quark mass effects in both the real and virtual corrections has been performed only recently in Refs. [33,34]. It shows that the total cross section is about 14% smaller than the one obtained within the Born-improved HEFT approximation and that for values of the Higgs boson pair invariant mass beyond M HH ∼ 500 GeV, the top-quark mass effects lead to a reduction of the differential cross section by about 20-30% with respect to the same approximation. Therefore, in order to get reliable predictions for the Higgs boson pair production cross section and corresponding distributions, it is important to include the full top-quark mass dependence.
Among the various kinematical distributions, a particularly significant role is played by the transverse-momentum (q T ) spectrum of the Higgs boson pair. A precise description of this observable is important to improve the statistical significance in the experimental searches and therefore, it is essential to carefully investigate the theoretical uncertainties dominated by the higher-order QCD corrections. In the large q T region (q T ∼ M HH ) fixed-order calculations are theoretically justified. However, in the small q T region (q T ≪ M HH ) the reliability of the fixed-order perturbative expansion is spoiled by the presence of large logarithmic terms of the type α n S log m (M 2 HH /q 2 T ) which make the fixed-order results divergent in the limit q T → 0. In order to obtain reliable predictions at small q T , such large logarithmic contributions have to be systematically resummed to all orders [35]- [42]. At intermediate values of q T the resummed and fixed order results can be consistently matched in order to get a uniform theoretical accuracy for the entire range of transverse momenta.
We have used the formalism introduced in Refs. [38,39] to perform the transverse momentum resummation for Higgs boson pair production up to next-to-leading logarithmic (NLL) accuracy, combining it with the NLO (i.e. O(α 3 S )) result with full top-quark mass dependence. The implementation of our calculation for the q T spectrum was performed starting from the numerical code HqT [39,43].
The paper is organised as follows. In Sect. 2 we briefly review the resummation formalism of Refs. [38,39,42]. In Sect. 3 we present numerical fixed-order and resummed results for the transverse-momentum distribution of Higgs boson pairs and we study the scale dependence of our results in order to estimate the perturbative uncertainty of our predictions. We also comment on the size of the finite top-quark mass effects. In Sect. 4 we present our conclusions.

Transverse-momentum resummation
The resummation formalism used in this paper has been introduced in Refs. [38,39] and can be applied to a generic process where a high-mass system of non strongly-interacting particles is produced in hadronic collisions. In this Section we briefly recall the main points of the formalism, by considering the specific case of the hadroproduction of Higgs boson pairs in gluon fusion. For a detailed discussion we refer to Refs. [38,39,40,41,42].
The transverse-momentum differential cross section for this process can be written as † : (1) where f a/h (x, µ 2 F ) (a = g, q,q) are the parton densities of the colliding hadrons (h 1 and h 2 ), dσ HH a 1 a 2 /dq 2 T are the partonic cross sections, M = M HH is the invariant mass of the Higgs boson pair, s (ŝ = x 1 x 2 s) is the hadronic (partonic) centre-of-mass energy, µ R and µ F are respectively the renormalisation and factorisation scale.
The partonic cross section is decomposed as follows: dσ HH a 1 a 2 = dσ (res.) HH a 1 a 2 , contains all the logarithmically-enhanced contributions at small † In this Section we denote with dσ HH /dq 2 T the double differential cross section M 2 dσ HH /dM 2 dq 2 T .
q T which have to be evaluated to all orders in α S and the 'finite' component, dσ HH a 1 a 2 , is free of such contributions.
The resummation procedure is carried out in the impact-parameter (b) space. The resummed component is obtained by performing the inverse Bessel transformation with respect to the impact parameter: where J 0 (x) is the 0th-order Bessel function. The resummation structure of W HH a 1 a 2 , N can be factorised and organised in exponential form by considering the Mellin N-moments W N of W with respect to z = M 2 /ŝ at fixed M ‡ : where we have defined the logarithmic expansion parameterL = ln ( . is the Euler number). The resummation scale Q [39] parameterises the arbitrariness in the separation (factorisation) between finite and logarithmically-enhanced terms. Variations of Q around the hard scale M can be used to estimate the effect of uncalculated higher-order logarithmic contributions.
The universal (process independent) form factor exp{G N } includes all the large logarithmic terms α n SL m , with 1 ≤ m ≤ 2n that order-by-order in α S are logarithmically divergent as b → ∞. The exponent G N can systematically be expanded in powers of α S ≡ α S (µ 2 R ) as follows: where the termL g (1) collects the leading logarithmic (LL) contributions, the function g N includes the NLL contributions, g N controls the NNLL terms and so forth. The logarithmic variableL is equivalent to L = ln (Q 2 b 2 /b 2 0 ) when Qb ≫ 1 (i.e. small values of q T ), but it leads to a behaviour of the form factor at small values of b such thatL → 0 and exp{G N } → 1 when Qb ≪ 1. The logarithmic expansion with respect toL thus reduces the impact of large and unjustified resummed contributions in the small-b region (i.e. at large values of q T ), and it acts as a perturbative unitarity constraint since it allows us to exactly recover the fixed-order value of the total cross section upon integration over q T .
The hard-collinear function H HH N fully encodes the process dependence of the resummation factor W HH N and it includes all the perturbative terms that behave as constants in the limit b → ∞. It has a customary perturbative expansion: Here, to simplify the notation, flavour indices are understood.
where σ (0) HH is the Born-level partonic cross section for the process gg → HH. The general structure of the hard-collinear function H F N has been obtained in Ref. [42], where it is shown that the process dependent contribution to H F N can be embodied in a single perturbative hard factor which is directly related to the finite part of the virtual amplitude of the corresponding process. The process independent part of the hard-collinear function H F N has been explicitly computed up to NNLO in Refs. [44,45]. For HH production the NLO corrections in the full theory were recently calculated [33,34]. From the values of the virtual amplitude with full topquark mass dependence computed in Refs. [33,34] we extracted numerically the process dependent contribution to the NLO coefficient H We now consider the finite component of the cross section. Since it does not contain large logarithmic terms, it can be computed at fixed order in perturbation theory starting from the standard fixed-order results and subtracting the expansion of the resummed component at the same perturbative order [39].
In summary, the resummation at NLL+NLO accuracy is obtained by including the functions g (1) , g S )) § . We note that the NLL+NLO result includes the full NLO perturbative contribution in the small-q T region and that the NLO result for the total cross section is exactly recovered upon integration over q T of the differential cross section dσ/dq T at NLL+NLO accuracy.
3 Numerical results for HH production at the LHC In this Section we consider Higgs pair production via gluon fusion in pp collisions at the centreof-mass energy of √ s = 14 TeV. We first show the fixed-order results which are valid in the large q T region and then we present our resummed prediction at NLL+NLO focusing on the small q T region. We include the full dependence on the top-quark mass and we comment on the size of the M t effects.
The hadronic cross section is computed using the PDF4LHC15 NLO parton densities [47] with α S evaluated at 2-loop order and α S (m 2 Z ) = 0.118 and we consider N f = 5 flavours of light quarks in the massless approximation. We set the central value of the renormalisation, factorisation and resummation scales at µ R = µ F = Q = M HH /2. The Higgs boson and top-quark masses have been set to M H = 125 GeV and M t = 173 GeV respectively, and the top quark and Higgs boson widths have been set to zero.
Bottom-quark mass effects in the double Higgs boson total cross section contribute well below 1% level and have been neglected in the present study. We thus have a two-scale problem with q T ≪ M, where M is the hard scale of the process M ∼ M HH and the top-quark mass is of the same order of the hard scale. This fact justifies the application of the standard q T resummation formalism to compute the HH q T spectrum with finite top-quark mass effects ¶ . § This matching procedure coincides with that of Refs. [39,46]. We note however that here we are using different labels. The fixed-order label NLO used here directly refers to the perturbative accuracy in the small-q T region and of the total cross section, while the labels LO and NLO used in Refs. [39,46] refer to the perturbative accuracy in the large-q T region. ¶ We note that the case of single Higgs boson production is different, since bottom-quark mass effects are sizeable and their inclusion leads to a three-scale problem [48].

blue solid), HEFT (black dotted) and Bornimproved reweighted HEFT (red dashed). The band is obtained by varying µ R and µ F as described in the text. Right panel: resummed prediction at NLL+NLO accuracy in the full theory. The resummed result (blue solid) is compared to the corresponding fixed-order result (red dashed) and to the finite component (black dotted). The lower panels show the ratios with respect to the central value of the full theory result.
We start the presentation of our numerical results by considering the calculation at fixed order. As explained in the previous Section we performed the calculation of the double Higgs boson q T spectrum at first order in QCD (i.e. O(α 3 S )). The relevant partonic subprocesses are gg → HHg, qg → HHq,qg → HHq and qq → HHg and we have generated all the relevant one-loop amplitudes using GoSam [49,50] retaining the full top-quark mass dependence. The corresponding matrix elements in the HEFT have been computed analytically. The phase space integration was performed using the CUBA library [51].
In Fig. 1 (left panel) we present the double Higgs boson q T spectrum for an invariant mass in the range 300 < M HH < 500 GeV at the LHC ( √ s = 14 TeV). We show the first order prediction in the full theory which includes the exact top-quark mass dependence (blue solid line) and in the HEFT (black dotted). In addition, we present also the approximation obtained by reweighting the HEFT result by the Born-level matrix elements for HH production with the full top-quark mass dependence (red dashed). The reweighting is performed at the matrix element level using the initial-initial antenna phase space mapping [52] to generate a Higgs boson pair Born-like configuration from the real-emission kinematics.
We show the scale dependence band (blue solid) of the full theory result which is obtained by varying independently the renormalisation and factorisation scales by a factor 2 around their central value, with the constraint 1/2 ≤ µ R /µ F ≤ 2. The scale dependence band is about ±35% at small q T , and it slightly increases up to ±40% at q T ∼ 400 − 500 GeV. We observe that the band is rather large and flat. This is not unexpected since the bulk of the scale dependence is due to the µ R dependence which is driven by the overall factor α S (µ R ) 3 .
By comparing the effective theory and the full theory results we observe that the pure HEFT calculation gives a poor approximation of the full theory spectrum for the entire range of q T (the two predictions accidentally cross each other for q T ∼ 250 GeV). The Born-improved reweighted HEFT result gives a good approximation of the exact calculation in the small q T region (q T ∼ < 50 GeV) where however both results diverge logarithmically. This is expected since in the small q T limit the phase space is restricted to soft and collinear emissions and in this limit the fixed-order cross section factorises into the Born contribution and process independent logarithmic terms. At larger q T (q T ∼ > 50 GeV), the agreement between the Born-improved HEFT and the full theory result rapidly deteriorates. The top-quark mass effects in the loop diagrams produce deviations in the fixed order q T spectrum of about 20 − 25% at q T ∼ 100 GeV and 80 − 100% at q T ∼ 175 GeV. The q T spectrum in the Born-improved HEFT approximation is much harder than in the full theory and generates an unphysical tail at large q T .
We now turn to present the resummed results. In Fig. 1 (right panel) we compare the NLL+NLO spectrum (blue solid line) at the default scales (µ F = µ R = Q = M HH /2) with the fixed-order result (red dashed). The finite component is also shown for comparison (black dotted). In the inset plot of the figure it is shown the region from intermediate to large values of q T . We observe that while the fixed-order calculation diverges at q T → 0, the resummation leads to a well-behaved distribution: it vanishes as q T → 0, has a kinematical peak at q T ∼ 18 GeV and tends to the corresponding fixed-order result for q T ∼ M HH . The finite component vanishes as q T → 0 and gives a contribution to the NLL+NLO result that is around 4% in the peak region and it increases to about 15% at q T ∼ 125 GeV and 30% at q T ∼ 200 GeV. We notice that in a wide region of intermediate values of q T the difference between the NLL+NLO and the fixed-order result is quite large (around 40 − 50% for 80 ∼ < q T ∼ < 250 GeV), thus indicating that the effect of the logarithmic terms included in the resummation is important even outside the small-q T region. The contribution of the finite component sizeably increases at large values of q T (q T ∼ M HH ) and the resummed spectrum approaches the fixed order prediction. We have checked the numerical accuracy of our calculation by computing the integral over q T of the NLL+NLO resummed spectrum. The result is in agreement with the value of the NLO total cross section calculated in Refs. [33,34] at the percent level, thus proving that the uncertainty associated to the numerical extraction of the H HH (1) N coefficient is completely under control.
We now discuss the scale dependence of the NLL+NLO result. As previously discussed the resummation formalism we are using is strictly valid for a two-scale problem with the resummation scale of the order of the top-quark mass. For this reason we explicitly avoided values of the resummation scale parametrically too large with respect to the top-quark mass M t . In Fig. 2 (left panel) we show the resummation scale dependence band (red dashed lines) obtained by varying Q in the region M HH /4 ≤ Q ≤ M HH /2 at fixed values of µ R and µ F (µ R = µ F = M HH /2). The resummation scale dependence is about ±3% at the peak, decreases to about ±1.5% at q T ∼ 30 GeV and increases again to about ±12% at q T ∼ 200 GeV.
Additionally we show in Fig. 2 (left panel) the resummed prediction for the resummation scale choice Q = M top (magenta dot-dashed line) and we quantitatively estimated the effect of this choice. We observe no significant differences with respect to the choice Q = M HH /2. The percentual difference with respect to our default value (Q = M HH /2) is around 1% at the peak, it decreases to few permille at q T ∼ 30 GeV, it increases to 2% at q T ∼ 100 GeV and it remains ∼ < 4% for 100 ∼ < q T ∼ < 300 GeV. This quantitative effect is widely covered by the resummation scale uncertainty band and can be explained by the fact that the HH cross section is peaked at an invariant mass of the order of M HH ≃ 400 GeV (for which M HH /2 ≈ M top ) [33,34]. In Fig. 2 (left panel) we also considered the renormalisation and factorisation scale dependence band (blue solid) obtained by varying independently µ R and µ F by a factor 2 around their central value (with the constraint 1/2 ≤ µ R /µ F ≤ 2) at fixed value of the resummation scale (Q = M HH /2). The µ R and µ F scale dependence band is about ±10% at the peak and it increases to about ±12% at q T ∼ 30 GeV, to about ±17% at q T ∼ 100 GeV and to about ±20% at q T ∼ > 250 GeV. We observe that the size of the µ R and µ F band is larger than the Q band for a wide region of q T (q T ∼ < 250 GeV). By comparing the µ R and µ F scale dependence bands of fixed order and resummed calculations, we observe that the resummed scale dependence band is not flat and its size is smaller than the fixed-order one at small and intermediate values of q T (q T ∼ < 250 GeV). This behaviour is not unexpected since the NLL+NLO resummed result, contrary to the fixed-order case, includes the full NLO correction in the small-q T region and satisfies the NLO unitarity constraint described at the end of Section 2 (see the discussion after Eq. (4)) and these NLO effects are spread on a region from small to intermediate values of q T . Nevertheless we point out that the µ R and µ F scale dependence is only a part of the perturbative uncertainty of the resummed prediction which includes also the resummation scale dependence.
We add a comment on the relation between the scale variation bands and the normalisation of the resummed q T spectra which is given by the corresponding NLO total cross section. On the one hand, the total cross section does not depend on the resummation scale. For this reason, the corresponding uncertainty band is independent on normalisation effects. On the other hand, the µ R and µ F scale variation band depends on normalisation effects and can be substantially reduced if we consider the normalised q T spectrum, 1/σ × dσ/dq T (i.e. if we are interested only on the shape of the q T distribution and not on its normalisation). The µ R and µ F scale dependence band for the normalised q T spectrum is shown in the lower left panel in Fig. 2 (black dotted lines), and we observe that it becomes smaller than the resummation scale uncertainty band for the entire q T range.
We conclude this Section with an assessment of the size of the finite top-quark mass effects which are included in our calculation. In order to study the impact of the M t effects we computed the resummed spectrum also in the HEFT approximation and in the Born-improved HEFT (the latter approximation was obtained by reweighting the HEFT result with the Born-improved factor at a differential level) . In Fig. 2 (right panel) we compare the NLL+NLO prediction in the full theory with exact M t dependence (blue solid line), with the pure HEFT (black dotted) and with the Born-improved HEFT (red dashed) results. We observe, similarly to the fixed-order case, that the Born-improved HEFT gives a good approximation (within 5% accuracy) of the full theory result for q T ∼ < 70 GeV * * . At higher values of q T we observe that the finite top mass effects are large and have a strong q T dependence. The effect is about 12% at q T ∼ 100 GeV, about 60% at q T ∼ 200 GeV and larger than 200% for q T ∼ > 250 GeV, showing that the inclusion of the full top-quark mass dependence is essential to obtain a reliable description of the double Higgs boson q T spectrum over a wide region of q T .

Conclusions
We have considered Higgs boson pairs produced in gluon fusion in hadronic collisions and we performed the calculation of the transverse-momentum (q T ) distribution of the double Higgs boson system taking into account finite top-quark mass (M t ) effects.
At small values of q T we have resummed the logarithmically-enhanced perturbative QCD contributions using the formalism introduced in Refs. [38,39]. We have presented the results of the resummed calculation at next-to-leading logarithmic accuracy (NLL), and we have combined them with the fixed-order computation at O(α 3 S ). Our calculation includes the complete next-toleading order (NLO) contributions at small q T and exactly reproduces the NLO total cross section with the full top-quark mass dependence upon integration over q T .
We have presented illustrative numerical results in pp collisions at √ s = 14 TeV, performing a study of the scale dependence of our predictions to estimate the corresponding perturbative uncertainty. Comparing the NLL+NLO and fixed-order results, we have shown that the higherorder terms contained in the resummed calculation are essential to obtain reliable predictions at small q T and give an important contribution ( ∼ > 40 − 50%) to the fixed-order result, in a wide region of intermediate values of q T (q T ∼ < 250 GeV). Finally, by comparing our results with the Born-improved Higgs effective field theory (HEFT) approximation in the M t → ∞ limit, we have quantified the size of the finite M t effects which turn out to be large ( ∼ > 60%) for q T ∼ > 200 GeV and very large ( ∼ > 200%) for q T ∼ > 250 GeV.
Our results show that both q T resummation and finite top-quark mass effects are necessary to obtain reliable predictions for the double Higgs boson q T spectrum over the full transverse We stress that the full theory result contains the complete finite M t effects both in the resummed part, through the Born-level partonic cross section σ  5)), and in the finite component. * * We note that the Born-improved HEFT approximation works particularly well (within 1% accuracy) for q T values around the peak. This agreement is not general and it depends on the particular Higgs boson pair invariant mass window. Considering Higgs boson pairs with an invariant mass in the range 350 < M HH < 400 GeV the agreement in the peak region is about 7%. momentum range.