Reconstruction of top-quark mass effects in Higgs pair production and other gluon-fusion processes

We propose a novel method for the treatment of top-quark mass effects in the production of $H^{(*)}$, $HH$, $HZ$ and $ZZ$ final states in gluon fusion. We show that it is possible to reconstruct the full top-quark mass dependence of the virtual amplitudes from the corresponding large-$m_t$ expansion and the non-analytic part of the amplitude near the top-quark threshold $\hat{s}=4m_t^2$ with a Pad\'e ansatz. The reliability of our method is clearly demonstrated by a comparison with the recent NLO result for Higgs pair production with full top-quark mass dependence.


Introduction
Gluons are ubiquitous at the LHC, and gluon fusion is among the phenomenologically most interesting production mechanisms. Specifically, the production of final states including one or more Higgs bosons is typically dominated by gluon fusion, with a virtual top-quark loop mediating the interaction to the Higgs bosons.
Precise predictions for such processes are indispensable for measuring the properties of the Higgs boson. On the one hand, gluon fusion processes experience large K-factors. 1 Examples include a K-factor of 2.3 for single Higgs and 1.7 for Higgs pair production at next-to-leading order (NLO) [4][5][6][7][8][9][10] which clearly demonstrates the importance of taking higher-order corrections into account. On the other hand, calculating these higher-order corrections is extremely challenging. Gluon fusion is a loop-induced process, and the top-quark mass introduces an additional scale in the loop integrals. While the NLO corrections to single-Higgs production have been known analytically for some time [5][6][7][8], the calculation of NLO corrections to processes with more than one final-state particle is still subject of on-going work. For di-Higgs production, which requires the evaluation of two-loop integrals with four scales, numerical results have only become available recently [9,10].
To make higher-order computations feasible an effective field theory (EFT), where the top quark has been integrated out in the limit of an infinite top-quark mass, m t → ∞, has been used extensively in the literature. In this approximation, results are available at NNNLO for single Higgs production [4,11] and at NNLO for Higgs pair production [12,13], and for other gluon fusion processes, i.e. gg → ZZ, gg → Hj at NNLO [14][15][16][17] and gg → HZ. Beyond the infinite top mass limit, several results have also been obtained in the large-m t expansion (LME) for a number of processes listed here: • gg → H: up to 1/m 6 t at NNLO [18][19][20][21][22], including gg → Hg at NLO • gg → HH: up to 1/m 12 t in [23] and 1/m 8 t in [24] at NLO; up to 1/m 4 t at NNLO [23] • gg → HZ: up to 1/m 8 t [25] at NLO • gg → ZZ: up to 1/m 12 t in [26] and 1/m 8 t in [27] at NLO The expansions can be rescaled with the exact leading order (LO) result dσ rescaled LME NLO /dX = dσ LME NLO /dX dσ LME LO /dX where dσ/dX indicates the differential cross section with respect to some quantity X. For inclusive Higgs production this yields good agreement with the exact NLO result [5][6][7][8]. The comparison with the exact Higgs pair production result has however revealed the shortcomings of the approximation (1) for this process [9,10]. This issue is especially pronounced when distributions are considered. Here, we advocate a different approach, based on conformal mapping and the construction of Padé approximations from expansions in different kinematical regimes of the amplitude. This strategy has first been introduced for heavy-quark current correlators Π (j) (q 2 /(4m 2 q )) [28,29] and applied successfully up to four-loop order [30][31][32]. The approximation can be improved systematically by including more information from the various kinematic limits. In fact, the three-loop approximation is indistinguishable from the results of an exact numeric computation [33]. In [28], it has also been shown for the decay H → γγ that a Padé reconstruction of the top mass effects from the asymptotic expansion in a large top mass yields excellent agreement with the full NLO decay rate. Like for heavy-quark correlators and the H → γγ decay rate, the amplitude for Higgs production in gluon fusion only depends on one ratio of scales and the application of the method is straightforward. However, the amplitudes for the remaining processes listed above depend on 4-5 scales. Padé approximations based on the LME terms alone have been used to reconstruct the interference contribution in gg → ZZ [26]. An attempt to reconstruct the gg → HZ cross section has been made in [25]. 2 In this work, we show how such an approximation can be improved drastically by also taking into account expansions in other kinematic regions, using Higgs pair production as an example.
Measuring di-Higgs production at the LHC allows to directly determine the trilinear Higgs boson self-coupling λ 3 [34][35][36], which serves as a probe of the shape of the Higgs potential and is a crucial test of the mechanism of electroweak symmetry breaking in nature. While the couplings of the Higgs boson to the gauge bosons and third-generation fermions have been firmly established to be Standard Model like within 10-20% [37][38][39], constraining the trilinear self-coupling is highly challenging. With 3000 fb −1 of data the estimated bounds are 0.2 < λ 3 /λ SM 3 < 7.0 (neglecting systematic uncertainties) [40]. Current bounds from Higgs pair production final states limit the trilinear Higgs self-coupling between −8.8 < λ 3 /λ SM 3 < 15.0 [41]. Under the assumption that only the trilinear Higgs self-coupling is modified, bounds can be obtained from single Higgs production through the electroweak corrections [42][43][44][45] or from electroweak precision observables [46,47]. However, the current bounds are still above the limits from perturbativity [48].
Precise theory predictions are crucial in the extraction of λ 3 from the cross section measurements. It is evident already at leading order (LO) that the LME alone is not sufficient. In fact, as shown in Figure 1, the cross section is dominated by energies of about 400 GeV, whereas the LME breaks down at the top pair-production threshold around 2m t ≈ 350 GeV. As we will show, constructing Padé approximations from the LME can ameliorate this problem to some degree, but not solve it completely. The reason for this is that, above the top threshold, the production amplitude receives non-analytic contributions, which cannot be reproduced by the purely rational Padé approximants. Incorporating these non-analytic threshold corrections enhances the quality of the approximation dramatically in the dominant kinematic region and thus leads to a much improved prediction for the total cross section.
The outline of this paper is as follows: In Section 2 we introduce our method for single Higgs production and then show how it can be generalized to the case of Higgs pair production. The computation of the additional input terms from the expansion around the top threshold is described in Section 3. In Section 4 we perform a detailed comparison of both the LO and NLO Padé approximation with the full LO result and the recent NLO results [9,10], respectively. We conclude in Section 5 and offer an outlook over possible applications of our method.

The method
We first discuss the construction of a Padé approximation for the simple case of the virtual amplitude A gg→H ( * ) in Section 2.1 and then generalize the approach to Higgs pair production in Section 2.2.

Padé approximation for gg → H ( * )
The LO diagram for the production of an off-shell Higgs in gluon fusion is shown in Figure 2 (left). The corresponding amplitude can be expressed through a dimensionless form factor F that only depends on the variable z = (ŝ + i0)/(4m 2 t ) whereŝ = (p 1 + p 2 ) 2 = p 2 H , y t = √ 2m t /v is the top Yukawa coupling, T F = 1/2 and The form factor F is normalized such that The leading-order contribution to the form factor is analytic in the entire complex plane with the exception of a branch cut for real z ≥ 1 due to on-shell tt cuts. At NLO, massless cuts like the one shown in the right of Figure 2 introduce a branch cut starting at z = 0. However, the branch cut can be made explicit t t H ( * ) H ( * ) Figure 2: The LO diagram for Higgs production in gluon fusion (left) and an example for a NLO diagram that contains a branch cut starting atŝ = 0 (right).
such that all the F il ,x (with i = 1, 2 and x = C F , C A , (C A , ln)) on the right-hand side are again analytic except for real z ≥ 1. In F 2l ,C A , IR divergences in the amplitude have been subtracted as described in Ref. [24]. We can now apply the conformal transformation [28] to map the entire complex z plane onto the unit disc |ω| ≤ 1 while the branch cut at z ≥ 1 is mapped onto the perimeter. The physical branch Im(z) > 0 corresponds to the upper semicircle, starting at ω(z = 1) = 1 and ending at ω(z → ∞ + i0) = −1.
With this mapping, the F il ,x are analytic functions of ω inside the unit circle. We approximate them using a Padé ansatz with a total of n+m+1 coefficients. They can be fixed by imposing conditions stemming from known expansions of the approximated function. In many cases it is found that diagonal Padé approximants with n = m provide the best description. Indeed, we find that this also holds for our analysis. We therefore discard approximants that are too far away from the diagonal, as detailed below.
The LME for the form factor F has been given up to terms of the order z 4 in [8]. The conformal mapping (6) transforms this into constraints on the derivatives of the Padé approximant at ω = 0. Furthermore the form factor vanishes for z → ∞ as F (z) = O(1/z) sinceŝ ∼ z has been factored out in (2). In a direct approach this would imply the constraint [n/m](ω = −1) = 0. Instead, we construct the Padé approximant for the rescaled form factor where a R is a free parameter. This serves a double purpose. First, it removes the spurious constraint at ω = −1 which implies that the dimensionality of the nonlinear system of equations that determines the coefficients of the Padé approximant is reduced by one. Secondly, the variation of the parameter a R allows us to test the stability of the ansatz and to assign an uncertainty to the reconstruction. A set of Padé approximants with n + m = 4 can be constructed based only on the constraints from the LME up to O(z 4 ). The Padé ansatz (7) has m poles in the ω plane. Here, and in the remainder of this work, we eliminate a subset of Padé approximants based on the positions of these poles. Since the amplitude is analytic inside the unit disc, the canonical selection criterion is to exclude approximants with poles at |ω| ≤ 1 + δ, where δ > 0 should be chosen such that no unphysical resonances, caused by nearby poles, are observed in the amplitude. We find, however, that this criterion proves too restrictive as it excludes almost all approximants. Thus, we relax the selection criterion and exclude approximants with poles in the region corresponding to values of z with 0 ≤ Re(z) ≤ 8 and −1 ≤ Im(z) ≤ 1, thereby excluding poles in the vicinity of the phenomenologically relevant region 0.13 z 5. We have checked the stability of the results under variation of the exclusion region. The result is shown in Figure 3 and compared to the exact expression for the form factor [5][6][7][8]. At LO the agreement is good, whereas at NLO the Padé curves become unstable under variations of a R and n/m and show significant deviations from the exact result for energies near and above the top threshold z 1.
We can gain some insight into this deviation by studying the expansion of the form factor around the top threshold. In particular we are interested in the nonanalytic terms in the expansion in (1 − z) which can be determined with the help of a factorization formula as discussed below in Section 3. Our results take the form F 2l F 2l where we have used the symbol to denote that terms that are analytic in (1 − z) have been dropped on the right-hand side. We observe that threshold logarithms ln(1 − z), which cannot be reproduced by the Padé ansatz, appear at NLO. Having determined the coefficients of the logarithmic terms at the first two orders we can Re Im LME Thr Figure 4: We show the same comparison as in Figure 3 but for Padé approximants based on the LME and the threshold expansion. Only [5/2] however subtract them from the form factor and apply the Padé approximation to the subtracted function. Taking a function f (z) with the threshold expansion as an example we definẽ where s 2,4 are constructed such that their leading non-analytic terms in the threshold expansion are given by respectively. In addition, the subtraction terms must be analytic around z = 0 and at most logarithmically divergent for z → ∞. 3 Apart from these constraints, the exact form of the subtraction functions is arbitrary. Our choice for the functions s 2,4 can be found in Appendix A. The threshold expansion off is free of logarithms up to and including the order (1 − z) 2 . An improved approximation of the original function f is then given by where the Padé approximant [n/m]f is constructed from the expansion terms of the subtracted functionf .
In addition, the non-integer powers of (1 − z) in eqs. (9) -(12) imply constraints on the derivatives of the Padé approximation at ω = 1. By using all the available constraints we can construct approximants with a total of n + m + 1 = 8 coefficients at LO and n + m + 1 = 7 coefficients at NLO. The results are given in Figure 4 and show perfect agreement with the exact LO form factor in the full energy range. At NLO the agreement is excellent up to z ∼ 2.5 where tiny deviations begin to emerge. For very large z, outside the phenomenologically relevant energy range, the approximants have unphysical extrema. We suspect that they could be removed by including information from the small m t expansion (SME) of the form factors in the construction.
An alternative implementation is obtained by performing additional subtractions for the root terms by employing the functions s 1,3,5 in Appendix A, thereby removing all known non-analytic terms in the expansion. This yields the same number of constraints on the Padé approximant. In the following we will only use the subtraction functions s 2,4 , since we find no significant differences between the two approaches.

Padé approximation for gg → HH
The amplitude for the process gg → HH can be parametrized by two dimensionless form factors F 1,2 with Given that there are four independent scales the dimensionless form factors depend on three ratios This implies that their analytic structure is much more complicated than it was the case for F . For instance, there are branch cuts in the complext andû planes above the thresholdst ≥ 4m 2 t andû ≥ 4m 2 t . These are, however, not kinematically accessible for external momenta that are both real and on shell. Furthermore, for z ≥ 1/r H ≥ 4 there is also a discontinuity from cuts corresponding to the processes gg → ttH and H → tt which are, however, not accessible for the physical Higgs and top masses. In the limit of small quark masses, z → ∞, where this type of cut is present, the recent analytical computation of the NLO virtual amplitudes for Higgs plus jet production [49,50] has revealed a rather complicated structure of logarithms in the soft and (in particular) the collinear limit which is presently not fully understood.
Here, we take a practitioners approach and note that when r H and r p T are kept fixed we can separate massless cuts as in (5) and again end up with functions that are analytic in z apart from a branch cut for real z > 1. Therefore it is possible to approximate the top-quark mass dependence of the form factors at a given phasespace point, i.e. for fixed m 2 H ,ŝ and p 2 T , by constructing a Padé approximant that describes the dependence on the variable z.
We find that the inclusion of the top threshold terms, as described for the triangle form factor (5) in Section 2.1, is of even greater importance for the construction of Padé approximants for the form factors (19) than for F . The computation of these terms is described in the following Section 3 and our results are given in Appendix C. Readers who are mostly interested in the phenomenological aspects may prefer to proceed to Section 4. There, we assess the reliability of our approach for Higgs pair production by performing a detailed comparison with the exact NLO results.

The amplitude near threshold
In this section the computation of the non-analytic terms in the threshold expansion of the form factors defined in Section 2 is described. Factorization formulae for the inclusive production cross section of heavy-particle pairs near threshold have been developed in [51][52][53][54][55] and applied to a number of processes [56][57][58][59][60][61][62][63]. The approach is based on the factorization of forward-scattering amplitudes which are related to the inclusive cross section by the optical theorem. We have extended the factorization formula to the gg → H ( * ) , HH, HZ, ZZ amplitudes. Only the basic aspects are sketched here and the reader is referred to the original literature [51][52][53][54][55] for a detailed derivation and discussion.

Structure of the amplitude near threshold
Near the threshold, z → 1, the top quarks can only be on shell if they are nonrelativistic. This implies a large hierarchy between the top mass m t , its typical momentum m t √ 1 − z and its kinetic energy m t (1 − z) which set the hard, soft and ultrasoft scale, respectively. Therefore, an effective field theory (EFT) can be constructed by integrating out the hard and soft scale. Then, the only dynamical modes left are non-relativistic top quarks, collinear and ultrasoft gluons and the external fields. The EFT describes the interactions of the remaining modes and is based on potential non-relativistic QCD (PNRQCD) [64][65][66][67][68][69] and Soft Collinear Effective Theory (SCET) [70][71][72][73][74][75]. The amplitudes for gg → F with final states F = H ( * ) , HH, HZ, ZZ are given by the master formula (cf. [51,52]) where the matrix elements have to be evaluated in the EFT. In analogy with [51,52] we call the contributions in the first and second line of (20) line the 'resonant' and 'non-resonant' amplitude, respectively. This structure is shown in Figure 5 in diagrammatic form.
The 'resonant' part in the first line of (20) contains the contributions that involve a non-relativistic top quark pair, i.e. a top pair that is close to being on resonance. This entails that only a soft spatial momentum can be exchanged between the initial and final state. Since the incoming gluons contain hard momentum components they must be connected by a hard subgraph. The same holds for the two final state particles. Integrating out these hard subgraphs yields local production operators that annihilate the incoming gluons and create a non-relativistic top pair and local annihilation operators O that annihilate the top pair and create the final-state particles. Here A ⊥ c is the collinear gluon field given in [54], the non-relativistic two-component spinor fields ψ and χ annihilate a top quark and produce an anti-top quark respectively, Γ (k) contains a combination of Pauli matrices, SU (3) c generators and potentially covariant derivatives and φ † F represents a combination of fields that produces the final state. Both types of operators have associated hard-matching coefficients that absorb the higherorder corrections from hard modes. The propagation of the non-relativistic top pair is subject to a non-local color Coulomb interaction that manifests as α s / √ 1 − z corrections in the amplitude. These so-called Coulomb singularities can be resummed to all orders within PNRQCD. The 'resonant' contribution contains non-analytic √ 1 − z and ln(1 − z) terms that correspond to on-shell cuts of the non-relativistic top pair.
Contributions where a hard momentum component is exchanged between the initial and the final state are contained in the 'non-resonant' part in the second line of (20). In the EFT they are represented by the matrix element of the local operator that annihilates the incoming state and creates the final state. Since the top quarks cannot be on shell near threshold when they carry hard momentum, there are no discontinuities associated with tt cuts. Therefore, this contribution admits the form of a Taylor expansion in (1 − z) once massless cuts have been separated as described in Section 2.1. The computation of this contribution is very involved since already the leading term in the Taylor expansion has the complexity of the full amplitude evaluated directly at the threshold z = 1. However, we expect the Padé approximation to predict this unknown analytic part of the amplitude very accurately, even when using only the LME as input. Indeed, as we showed explicitly in Section 2.1, adding the knowledge of just the non-analytic terms near threshold is already sufficient to reconstruct the full top-quark mass dependence with high accuracy. Therefore we can safely ignore the non-resonant contribution and only focus on the much simpler factorizable part.

Computation of the non-analytic terms
In this section we describe the computation of the 'resonant' part of the amplitude (20). We adopt here the non-relativistic power counting where α s ∼ √ 1 − z and denote the k'th order in this counting by nrN k LO to distinguish it from the fixed-order expansion in the strong coupling constant. At nrLO, the matrix element is given by a non-relativistic Green function which resums the 1/ √ 1 − z enhanced effects from the ladder-exchange of Coulomb gluons as indicated in Figure 6. Hence, at any loop order, the leading non-analytic term in the threshold expansion of the amplitude can be determined by expanding the nrLO result to the respective order in α s . Up to nrNNLO, terms of the relative order = . . . must be included, where A 0 (z = 1) is the LO amplitude evaluated at the top threshold, l = 0, 1, . . . denotes the angular momentum of the top pair and the global factor √ 1 − z accounts for the suppression of the phase-space near threshold. Fig. 7 illustrates the relation between different orders in standard relativistic perturbation theory and in the non-relativistic effective theory. For example, the following terms on the right-hand side of Eq. (24) contribute to the fixed-order expansion up to NLO: • The nrLO terms with relative factors • The nrNLO terms with relative factors • The nrNNLO terms with relative factors For the processes gg → H ( * ) and gg → HH there is no contribution from Swave tt states due to parity and C-parity conservation. 4 The leading 'resonant' contribution therefore contains the P-wave Green function [76] which is suppressed by (1 − z) near threshold. We want to determine the 'resonant' amplitude up to nrNLO in the scaling (24), which contains the next-to-leading non-analytic terms in the threshold expansion at any loop order. In addition we compute the first two terms in the fixed-order expansion of the nrNNLO result in α s , i.e. those of relative orders (1 − z) 5/2 and α s (1 − z) 2 . They correspond to the next-to-next-to leading threshold terms for the one and two loop amplitude which we study in Section 4.
The matrix elements in (20) receive corrections from the higher-order non-local potentials and the dynamical modes contained in the EFT. The EFT contains no interactions of collinear modes with non-relativistic modes or between collinear modes of different directions. They cannot be present because the combination of the involved momenta yields hard modes which have been integrated out. Therefore the LO NLO NNLO nrLO nrNLO nrNNLO Figure 7: Relation between relativistic (LO, NLO, NNLO) and non-relativistic (nrLO, nrNLO, nrNNLO) power counting up to next-to-next-to-leading order. The axes show the powers of α s and √ 1 − z in the various coefficients represented by the markers. Note that the normalization is chosen such that α 0 s corresponds to LO.
only collinear corrections at nrNLO are from the left diagram in Figure 8. The corresponding loop integral is scaleless and therefore vanishes in dimensional regularization.
Ultrasoft gluons couple to the collinear and non-relativistic sector as well as to the P-wave production and annihilation operators. The exchange of ultrasoft gluons between the collinear states shown in the diagram on the right of Figure 8 yields only scaleless integrals. The interactions in the EFT must be multipole expanded. At leading order in the multipole expansion ultrasoft gluons couple to the net color charge of the tt state since the large wavelength λ ∼ 1/(m t (1 − z)) gluons cannot resolve the spatial separation a B ∼ 1/(m t √ 1 − z) of the top pair. The first non- vanishing term in the multipole expansion for color singlet states is therefore the chromoelectric term ψ † x · E ψ which is suppressed by α s . Similarly the ultrasoft gluon term in the covariant derivative in the P-wave operators is suppressed by α with respect to the derivative term. A single insertion of either of these subleading terms vanishes by rotational invariance [55]. Thus, contributions from the subleading interactions require at least two insertions and first appear at nrNNNLO.
The effects of higher-order potentials enter as corrections to the non-relativistic Green function. The nrNNNLO S-wave and nrNLO P-wave Green functions have been computed for tt production in e + e − collisions near threshold [76,77]. We determine the α 0,1 s terms in the nrNNLO P-wave Green function in Appendix B. Up to the considered order the resonant amplitudes hence take the simple factorized form The Wilson coefficients C  As mentioned before, it is sufficient to know the nrNNLO terms proportional to α 0 s and α 1 s in order to construct approximations to two-loop (NLO) fixed-order amplitudes (cf. Fig 7). The remaining nrNNLO terms of the relative order α 2 s (1 − z) 3/2 will be important for the reconstruction of the three-loop amplitude. Since its determination requires the calculation of the two-loop matching coefficients C tt→F are finite after field and mass renormalization. The one-loop coefficients C (k) gg→tt , however, require additional IR subtractions since the virtual amplitude by itself is not IR safe. Our results for the threshold expansion of the form factors are given in (9)- (12) and Appendix C together with the details of the IR subtractions. Together with the nrNLO expression for the P-wave Green function [76] these results are sufficient to determine the leading and next-to-leading non-analytic terms in the threshold expansion of the form factors at any order in α s .
Another interesting, yet more involved, application of our formalism is Higgs plus jet production. Here, we shortly comment on that, but leave a more careful assessment to future work. The amplitudes gg → Hg, gq → Hq and qq → Hg obey the same structure of (20) near the top threshold but the corresponding 'resonant' matrix elements are more complicated since the final state now contains a colorcharged particle. Ultrasoft gluons can then be exchanged between the initial state, the final state and the intermediate top pair which is in a color octet state and no longer decouples. In [53][54][55] it was demonstrated for arbitrary color structures that the 'resonant' matrix elements in forward-scattering amplitudes factorize into the convolution of a non-relativistic Green function, therein called the potential function, and an ultrasoft function, therein called the soft function. At leading power this follows from field transformations that decouple the collinear and non-relativistic fields from the ultrasoft fields. The extension to higher orders requires a careful assessment of the subleading interactions and was performed to NNLL in [53][54][55]. Following these derivations we identified no aspect that would obstruct the extension to Higgs plus jet production and therefore conjecture that an analogous factorization formula holds for the corresponding amplitudes.

Comparison with the exact result
As a proof of method, we compare our results at LO and NLO with the results in full top mass dependence for Higgs pair production. While at LO, the Higgs pair production cross section is known in full mass dependence since the late 80's [78][79][80], the computation of the NLO QCD corrections is quite involved, due to the many scales of the problem. The first work on the NLO corrections was based on the heavy top mass limit [81] reweighted with the matrix elements squared of the full LO results (HEFT). The real corrections in full top mass dependence have been computed in [82,83], while the virtual corrections have been kept in HEFT. The computation of the virtual corrections in full top mass dependence became available only recently in [9,10].

Numerical setup
For the numerical evaluation we choose a centre-of-mass energy of √ s = 14 TeV.
The Higgs boson mass has been set equal to m H = 125 GeV and the top quark mass to m t = 173 GeV. We do not account for bottom quark loops as they contribute with less than 1% at LO. We have adopted the PDF set NNPDF3.0 [84]. The strong coupling constant is set to α s (M Z ) = 0.118 at LO and NLO. The renormalization scale has been set to M HH /2, where M HH denotes the invariant mass of the Higgs boson pair, as suggested by the NNLL soft gluon resummation performed in [85,86].
We construct our Padé approximants at LO (NLO) as described in Section 2 by solving numerically the 8 (7) equations from the LME [24] and threshold expansion, given in Section 2.1 and Appendix C, by means of the FORTRAN routine MINPACK [87]. 5 For every phase space point we construct a total of 100 Padé approximants [n/m], where a R takes a random value between [0.1,10], n, m ∈ [1,6] at LO and n, m ∈ [1,5] at NLO, and take the mean value. From that we obtain an error estimate on every form factor by taking the standard deviation. For the computation of the cross section or the virtual corrections we add up the errors stemming from the different form factors quadratically. Padé approximants with poles in Re(z) ∈ [0, 8] and Im(z) ∈ [−1, 1] were excluded, since functions with poles close-by in the complex plane could have an unwanted resonant behaviour. The running time per phase space point for the construction of 100 Padé approximants at NLO is usually below 6 s.

Comparison at LO
In Table 1 we give the results for the LO cross section in different approximations. The first row, [n/m] w/o THR, symbolizes the cross section obtained with Padé approximants constructed without input from the threshold expansion, where n, m ∈ [1, 3] and approximants with poles as described above have been excluded. The result we obtain when including the threshold information and using the specifications described in Section 4.1 is denoted by [n/m]. With [n/n ± 1, 3] we symbolize the results we find when only the Padé approximants [5/2], [4/3], [3/4] and [2/5] are used. 6 Finally, we give the full LO cross section (obtained with HPAIR [88]) in the fourth row of Table 1. As can be inferred from the table, the Padé approximants provide a very good approximation for the full cross section, in particular if only the most diagonal and next-to-diagonal Padé approximants are constructed. The threshold expansion proves to be essential for a good approximation. As expected, the standard deviation computed from the construction of 100 [n/m] Padé approximants with random a R and different n, m becomes smaller if we construct only the most diagonal and next-to-diagonal Padé approximants.
In Fig. 9 we show the invariant Higgs mass distribution for the full result (dark blue), the [n/n ± 1, 3] Padé approximants (pink) and the Padé approximants without the threshold expansion (light blue). While the [n/n ± 1, 3] full Padé approximants fit the shape of the invariant mass distribution in full mass dependence almost per-  fectly, the approximation where the threshold expansion is not included (hence the approximation is only built from the LME) fits the shape only for small invariant mass. The error on the construction of the approximation including the threshold expansion is rather small whereas if the approximation is constructed only from the LME, the error becomes much larger in particular above the threshold. We thus conclude that at LO our approximation of the mass effects by Padé approximants works well as long as the conditions obtained from the threshold expansion are included. Using only nearly diagonal Padé approximants leads to a result with smaller error with values closer to the true result.

Comparison at NLO
Finally, we compare our results to the computation of the NLO corrections in full top mass dependence of Ref. [9,10]. In the framework of Ref. [89] a grid and an interpolation function with numerical values for the virtual corrections of Ref. [9,10] have been provided. In order to fit the conventions of Ref. [89] we define the finite part of the virtual corrections as with and F 1 defined in eq. (36). For F 2l,[n/m] x we use the matrix elements constructed with the Padé approximant [n/m]f . All other matrix elements are used in full top mass dependence. The form factors F 2∆ x stem from the double triangle contribution to the virtual corrections and can be expressed in terms of one-loop integrals. They are given in Ref. [24] in full top mass dependence. In the heavy top mass limit they become The contribution of the double triangle diagrams to the virtual corrections is only of the order of a few per cent [90].
In Table 2 we compare values for the full computation of the virtual corrections obtained from the grid of Ref. [89], the HEFT results rescaled with the full Born  ). The errors given in the table are, in case of the Padé-approximated results, due to the construction of the different approximants and due to the rescaling with a R . For the full results the error stems from internal binning in the grid. As can be inferred from the table, the Padé construction approximates the full result quite well. It provides a much better approximation than the HEFT results with a generally reliable error estimate. While up to M HH = 450 GeV the Padé method provides an excellent approximation on the level of 2%, for larger invariant masses and p T the results worsen gradually. As already anticipated from the LO results, constructing only diagonal and next-to diagonal Padé approximants improves both the error and the values of the virtual corrections with respect to the full result. Indeed we even find that only constructing diagonal Padé approximants gives results even closer to the full result. Since this does not allow for a reliable error estimate any more (the error would then solely stem from the variation of a R ) we do not discuss this here any further.
In Fig. 10 we show for p T = 100 GeV the virtual corrections V f in for varying M HH for the Padé approximations [n/n ± 0, 2], the Padé approximants constructed only from the LME, the full result and the reweighted HEFT results. Again, we can see that contrary to the HEFT results the Padé approximation can reproduce the correct scaling with the invariant mass of the full result. The quality of the approximation is improved significantly with the inclusion of the threshold expansion. The error of the Padé approximation increases with the invariant mass. Note that the full result has, apart from the previous error from the internal binning, also an error due to the interpolation procedure. We do not quantify this error but in comparison to the HEFT grid provided with Ref. [89] we conclude that while in the range up to M HH 570 GeV this error is negligible, it will be a few % for larger M HH . The comparison with the numerical results of [89] demonstrates that our prescription for the uncertainty related to the construction of Padé approximants also provides a reasonable error estimate at NLO.
In conclusion, we see that for the NLO corrections the Padé approximation reproduces the correct scaling behaviour for small and moderate invariant mass and p T . Since the cross section peaks around M HH ≈ 400 GeV and p T ≈ 150 GeV this will lead to a reliable approximation and reliable error estimate also for the full cross section. It can be expected that both the error and the difference with respect to the full result improves once more input is used (i.e. higher orders in the threshold expansion, higher orders in the LME, possibly input from a small mass expansion).

Conclusions and outlook
We have reconstructed the top-quark mass dependence of the one and two loop virtual amplitudes for Higgs pair production in gluon fusion with Padé approximants based on the LME of the amplitude [24] and new analytic results near the top thresholdŝ = 4m 2 t . We observe perfect agreement of the one-loop results with the exact expressions once the additional conditions from the threshold terms are imposed. Significant deviations are observed when only the LME is used to construct Padé approximants, but we still find agreement within the uncertainty estimate of our reconstruction, which is based on variation of the rescaling parameter a R and the use of different [n/m] approximants. At the two-loop level the full result can be reproduced in the entire phenomenologically relevant range within typical uncertainties ranging from below ±3% in the region M HH ≤ 450 GeV up to about ±20% for M HH = 700 GeV. Thus, our method allows for a determination of the total cross section including top-quark mass effects at NLO where the uncertainty due to the reconstruction is negligible compared to the scale uncertainty which is of the size of ±13% [9,10]. This represents considerable progress compared to the rescaled HEFT and LME approximations where a reliable uncertainty estimate is not possible. Our method can also be systematically improved by including higher orders in the LME or threshold expansions. We expect even better behaviour if one also considers the leading term in the small-mass expansion z → ∞ which corresponds to the bottom-quark contribution expanded for small m b . An approach for computations in this limit has recently been introduced [49,50,91]. Furthermore our results strongly suggest that the combination of the Padé approximants of the NLO virtual corrections with the exact evaluation of the real corrections [82,83] can reproduce differential distributions to high accuracy.
There is a large number of possible applications for our method. To further increase the precision for Higgs pair production one needs to consider NNLO QCD corrections. The rescaled HEFT approximation for the NNLO corrections increases the cross section by 18% [10] which exceeds the estimate from scale variation at NLO. A NNLO computation which retains the full top-quark mass effects is clearly out of reach of the current technology. On the other hand, the LME has already been computed up to 1/m 4 t in [23] and we have determined the two first non-analytic terms in the threshold expansion. This presently available input only allows for the construction of Padé approximants with n + m = 3 where we do not expect stable behaviour, but a calculation of two or three more expansion parameters would allow the evaluation of NNLO corrections in the soft-virtual approximation of [23,92]. Additionally, one can study the NLO electroweak corrections involving topquark loops. Of particular interest are the contributions involving additional Higgs bosons which alter the dependence of the cross section on the values of the Higgs self couplings.
It is straightforward to apply our method to gg → HZ and the top-quark mediated gg → ZZ amplitude and at higher orders in perturbation theory. In all these cases, results in the LME have been obtained at two loops [25][26][27] and for gg → H ( * ) even at three loops [18][19][20][21][22]. The determination of the threshold terms only requires the computation of the respective one-loop matching coefficients in (25). Another phenomenologically very interesting case is Higgs plus jet production. The construction of Padé approximants is also possible here but the computation of the threshold expansion is more involved as we outlined in Section 3.2. Beyond LME results, also the leading term in the small-mass expansion is know for the relevant two-loop amplitudes [49,50]. Hence, the effects of this additional input on the reconstruction of top-quark mass effects can be studied in this case.
where we have used the symbol to denote that terms analytical in (1 − z) have been dropped on the right-hand side, and Π (1),v is the well-known two-loop correction to the vacuum polarization [93] in the convention of [31]. The functions s i in (29) are constant as z → 0 and only diverge logarithmically as z → ∞.

B Expansion of the P-wave Green function
The P-wave Green function has been computed up to nrNLO in [76]. In addition we have determined the terms of order α 0 s and α 1 s in the nrNNLO correction. Those are given by the insertion of the 'kinetic potential' [69] and the 1/m 2 potential [69] where q = p − p , the term proportional to V 1/m 2 vanishes for the P-wave due to asymmetry under the integration over the spatial momentum components and V p = 1 + O(α s ). Our result for the P-wave Green function expanded in α s and (1 − z) reads where and again indicates that terms analytic in (1 − z) have been dropped.

C Results for the gg → HH form factors near threshold
We give the results for the threshold expansion of the gg → HH form factors up to three-loop order. The expansion of the form factors F 1 and F 2 in the strong coupling constant takes the form At the two-loop level the contributions F 2 i that involve two top-quark loops are known exactly [24] and have therefore been separated in (35). They will not be considered further because their threshold expansion does not contain any non-analytic terms. The form factor F 1 is further decomposed into a 'triangle' and 'box' contribution as indicated in Figure 11. The contributions of the 'triangle' diagrams to the form factor F 2 vanish. As discussed in Section 2 we make massless cuts explicit such that all the F 's on the right-hand side of (37) and (38) are analytic in z except for a branch cut for real z ≥ 1.
'Triangle' topology 'Box' topology Like the previous works [5,6,24] we use the following MS scheme convention in our calculation. The renormalized form factors still contain IR divergences which cancel with contributions involving unresolved real radiation that are not considered here. We use subtractions of a minimal type as Refs. [5,6,24] The full form of the subtraction term at NNLO is known [23,94] and includes a contribution proportional to F 1l i which has been omitted here because it only affects the three-loop results beyond the considered order in the threshold expansion.
Our results for the triangle form factor at one and two loops are given in (9)(10)(11)(12). At the one-loop order we determine the remaining form factors up to nrNNLO in the threshold expansion In spite of some notational differences, our convention is identical to [24]. There is an exact cancellation between the β 0 / contribution in (40) and the charge and gluon field renormalization terms. This has been exploited in [24] where both effects are not written explicitly. Furthermore, the sign in the factor (−ŝ − i0) − has been ignored in [24] because the induced imaginary part is not relevant within the LME at the considered order in the strong coupling constant.