Real corrections to Higgs boson pair production at NNLO in the large top quark mass limit

In this paper we consider the next-to-next-to-leading order total cross section of Higgs boson pair production in the large top quark mass limit and compute four expansion terms in 1/mt2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_t^2 $$\end{document}. To this end, we analytically compute the real-virtual and double-real contributions to the total cross section and combine them with the existing virtual contribution. Good convergence is observed below the top quark threshold, which makes our results a valuable input for approximation methods which aim for next-to-next-to-leading order corrections over the whole kinematic range. We present details on various steps of our calculation; in particular, we provide results for three- and four-particle phase-space master integrals and describe in detail the evaluation of the collinear counterterms.


JHEP01(2022)049
momentum, say below p T ≈ 200 GeV, and the fast evaluation of the analytic high-energy expansions could be used for the remaining phase space.
The dependence on the renormalization scheme of the top quark mass has been discussed in ref. [26]. Sizeable uncertainties were observed for large values of the di-Higgs invariant mass. To reduce this uncertainty a NNLO calculation is necessary. A full NNLO calculation is currently out of reach, so it is important to construct approximations valid in different regions of the phase space. In the region where the large-m t expansion -the topic of this paper -is valid the scheme uncertainty is small. It is nevertheless useful to have such corrections at hand since they serve as valuable input for approximation methods such as, e.g., the one described in ref. [18].
At NNLO, only approximations in the large-m t limit are available so far. In the infinitem t limit the cross section has been computed in refs. [27][28][29] and finite 1/m t expansion terms for the virtual corrections are available from [30,31]. In ref. [32] a NNLO approximation has been constructed which incorporates, in addition to the exact NLO corrections and the infinite-m t results at NNLO, the exact expression for the double-real radiation contributions at NNLO.
In the present work we consider the 1/m t corrections at NNLO. The virtual corrections to the form factors have been computed in [31] including terms to order 1/m 8 t . In section 2.6 we describe our procedure to obtain the corresponding contribution to the cross section to this order. The main purpose of this paper is to complement the virtual corrections by the corresponding real-radiation contributions. We compute all contributions up to order 1/m 6 t and the channels which involve quarks to 1/m 8 t . The real-radiation contribution can be divided into one-loop 2 → 4 processes and twoloop 2 → 3 processes; sample Feynman diagrams are shown in figure 1. In the following we denote these contributions as "real-real" and "real-virtual", respectively. We are interested in the total cross section which is conveniently obtained using the optical theorem applied to the forward-scattering amplitude (see figure 2). At NNLO this leads to five-loop diagrams. However, in the large-m t limit one observes a factorization into n-loop tadpole and (5 − n)loop so-called phase-space integrals, where in our case we have n = 2 or n = 3. The tadpole integrals are well studied in the literature, however the phase-space integrals are not; we discuss in detail the various integral families, their reduction to master integrals and provide analytic results.
We want to mention that a part of the real-radiation contribution, where the two Higgs bosons couple to different top quark loops (cf. figure 1(c)) have been computed in ref. [31]. In the forward-scattering picture of figure 2 such contributions have three closed top quark loops and we denote them the "n 3 h contribution". Note that there are further n 3 h contributions which have a closed loop with only gluon couplings (as shown in figure 1(d)). Such terms are not included in ref. [31], but are computed in this paper. The remainder of the paper is organized as follows: in the next section we discuss the individual parts of our calculation. This concerns in particular the setup used for the JHEP01(2022)049 Figure 1. Sample Feynman diagrams contributing to the real radiation. Contributions such as those shown in (c) lead to n 3 h contributions which have already been computed in ref. [31]. The n 3 h contributions of (d) contain a top quark loop without a Higgs coupling and have not been computed in ref. [31]; they are considered here.
(a) (b) (c) (d) Figure 2. Sample Feynman diagrams in the forward-scattering kinematics. Three-and four-particle cuts are shown by blue and green dashed lines, respectively. The n 3 h contributions as shown in (b) have already been considered in [31] but those in (c) have not; they are considered here.
computation of the real-radiation corrections including the asymptotic expansion and the reduction to phase-space master integrals. Furthermore, we discuss the ultraviolet and collinear counterterms to subtract the divergences from initial-state radiation. Section 3 is dedicated to the phase-space master integrals. We provide details on the transformation of the system of differential equations to form and on the computation of the boundary conditions in the soft limit. We discuss our analytic and numerical results in section 4 and summarize our findings in section 5. In the appendix we provide useful additional material such as explicit formulae used for the computation of the collinear counterterms, the integrands of the phase-space master integrals, NNLO virtual corrections to the channel qq → HH and NNLO virtual corrections involving four closed top quark loops. Furthermore, we describe in detail our approach to obtain the leading 1/m t term for double Higgs production from the analytic expressions of the single-Higgs production cross section.

Partonic channels and generation of amplitudes
We write the perturbative expansion of the partonic cross sections as σ ij→HH+X (s, ρ) = δ ig δ jg σ (0) gg + σ (1) ij + σ (2) ij . . . , (2.1) where the superscripts (0), (1) and (2) (2) ij = σ (2),n 2 h ij + σ (2),n 3 h ij + δ ig δ jg σ (2),n 4 h gg . (2.2) Note that the superscripts n k h counts the number of closed top quark loops which have at least one coupling to a Higgs boson. The contributions σ (1),n 3 h gg and σ (2),n 3 h ij have been computed in ref. [31]. In this paper we concentrate on σ (2),n 2 h ij which can be decomposed as σ (2),n 2 h ij = σ (2),n 2 h ij,virt + σ (2),n 2 h ij,real + σ (2),n 2 h ij,coll , (2.3) where the individual contributions are discussed in the remainder of this section. Most emphasis is put on σ (2),n 2 h ij,real since the complicated calculation for σ (2),n 2 h ij,virt was performed in [31] and for σ (2),n 2 h ij,coll "only" convolutions of lower-order cross sections and splitting functions have to be computed. σ (2),n 4 h gg is discussed in appendix D. We compute the cross section in the large-m t limit. It is thus convenient to introduce the variable such that after the expansion in 1/m t the coefficients of ρ n depend on x = m 2 H /s. The analytic calculation of the phase-space integrals is performed in an expansion around the production threshold of the two Higgs bosons, for which it is convenient to define δ = 1 − 4x , (2.5) such that the threshold corresponds to δ = 0.
To compute σ (2),n 2 h ij,real we first generate the forward-scattering amplitudes for the relevant partonic channels, applying the optical theorem to obtain the real-radiation contributions. These amplitudes involve two external momenta q 1 and q 2 , with s = (q 1 + q 2 ) 2 . One has to consider the gg, gq, qq, qq and qq partonic channels, where q and q denote two different light quark flavours. Note that both quark and anti-quark contributions have to be taken into account, i.e., gq and qq also stand for the gq andqq contributions, respectively. Similarly in qq , all quark-quark, quark-anti-quark and anti-quark-anti-quark contributions have to be considered, where in each case the quark flavours are different. For all channels one obtains the same result after replacing a quark by an anti-quark and vice versa. We also compute channels with one or both external gluons replaced with ghosts (c) and anti-ghosts (c), i.e., in addition to the channels listed above we have gc, gc, cc, cc,cc, cq andcq. This allows us to sum over the gluon polarizations according to where λ runs from 0 to 3.

JHEP01(2022)049
Channel qgraf diagrams gen-filtered diagrams building block diagrams gg 16 We generate the amplitudes using qgraf [40], however this program does not allow one to filter only diagrams which admit a cut through the required final-state particles (here two Higgs bosons and one or more gluons or light fermion pairs), as depicted in figure 2 by the green and blue dashed lines. We thus generate all possible forward-scattering diagrams and post-process them using gen [41] which is able to filter the qgraf output based on our required cuts.
Due to the large qgraf output (for the gg channel, 16.6 million diagrams, 42 GB in total) it proved necessary to separate the diagrams into subsets containing particular numbers of top quark, light quark, cut-and un-cut Higgs lines which could be filtered by gen separately, in parallel. After filtering, the number of diagrams remaining is greatly reduced (for the gg channel, to a total of 160 thousand diagrams, 416 MB in total). In table 1 we show the numbers of diagrams for all channels which we compute.

Real radiation: asymptotic expansion and building blocks
After producing the filtered set of Feynman diagrams, we consider the large-m t limit and apply an asymptotic expansion for m 2 t s, m 2 H . This expresses each diagram as a sum of products of so-called "hard subgraphs" (which have to be expanded in their small quantities: masses and momenta) and "co-subgraphs" which are obtained from the original diagrams by shrinking the subgraph lines to a point. The rules which must be applied to determine all of the relevant subgraphs of a diagram are implemented in the software exp [42,43], which produces FORM [44] code to perform the expansion.
The subgraphs generated by the asymptotic expansion are, in this case, one-and twoloop one-scale vacuum integrals; this class of integrals is well studied, and it is possible to compute tensor integrals which contract with external momenta or line momenta of the cosubgraph. Such routines are implemented in MATAD [45]. The co-subgraphs of the expansion lead to two-and three-loop forward-scattering integrals (or "phase-space" integrals) which depend on s and m 2 H . The two-loop phase-space integrals involve three-particle cuts whereas the three-loop phase-space integrals involve both three-and four-particle cuts and JHEP01(2022)049 Figure 3. The forward-scattering diagrams of figure 2, with black circles representing the expanded hard subgraphs in the large-m t expansion. Three-and four-particle cuts are shown by blue and green dashed lines, respectively. thus contribute to both the real-virtual and real-real contributions. From the diagrams in figure 2 one obtains the diagrams in figure 3 where the black circles represent the expanded subgraphs.
In principle, the computation can now be performed using the method and tools described above and in section 2.1, that is, 1. generate the diagrams with qgraf and filter for valid cuts with gen, 2. use q2e/exp to identify subgraphs and generate FORM code, 3. compute the asymptotic expansion with FORM and MATAD, yielding expressions in terms of the remaining phase-space integrals which must still be computed. In practice however, this approach is very computationally challenging due to the enormous number of five-loop forward-scattering diagrams generated and the complexity of expanding each of those diagrams due to the large number of propagators present. We use this method to compute in full only the leading contribution to the large-m t expansion (terms of order 1/m 0 t ) in order to cross-check an alternative approach described below. For the channels gq, qq and cc, we additionally perform this cross-check to order 1/m 2 t in the large-m t expansion.
Many of the diagrams have identical subgraphs, which are expanded in the same small parameters. This suggests that it is beneficial to pre-compute the subgraphs, expanding them to the required order only once and storing the results. We can then generate two-and three-loop forward-scattering amplitudes with "effective vertices" in place of the subgraphs, of which there are comparatively very few (for the gg channel we generate 4612 diagrams, compared to the 160 thousand diagrams discussed in section 2.1. In table 1 we show the numbers of diagrams for all channels). As an example, in figure 4 we show diagrams which involve the 4-gluon-2-Higgs effective vertex; there are 3600 "full" diagrams, one of which is shown in figure 4(a), and just one effective-vertex diagram shown in figure 4(b). This subset of diagrams has been discussed in ref. [46]. In the course of the calculation, the effective vertices are replaced with the pre-expanded subgraphs. We call this the "building block approach".
In practice, we adopt something of a hybrid approach in order to avoid having to deal with two-loop building blocks which are challenging to compute. This means that we implement the following one-loop building blocks: ggH, gggH, ggggH, ggHH, gggHH, JHEP01(2022)049   ggggHH, shown graphically in figure 5. We compute them to order 1/m 8 t with all external momenta off shell, and dependence on the open external Lorentz and colour indices retained. The resulting expressions range between 84 kB for ggH and 164 MB for ggggHH. This means that they have to be inserted into the phase-space diagrams carefully, to obtain acceptable performance. Initially the coefficients of each term of the 1/m 2 t expansion are inserted only as placeholder symbols, while the rest of the diagram is computed. Once the expression is written in terms of a linear combination of phase-space integrals, the placeholder symbols are replaced by using Term environments in FORM to sort the coefficients of the phase-space integrals and 1/m t expansion coefficients individually, so as not to blow up the intermediate size of the whole expression.
For the building blocks involving up to three gluons, the Lorentz and colour structures factorise and can thus be inserted separately into the computation of the kinematic and colour parts of the diagrams, as is required by the structure of our setup to compute the diagrams. For the four-gluon building blocks this is no longer the case, so we proceed in the following way. First, we simply compute the expansion of the building blocks to arrive at a result in the form where the factors K i contain the kinematic dependence. {p k } is the set of external momenta of the building block. The colour structures are given by the C i factors. We now define the symbol δ i such that

7)
JHEP01(2022)049 Figure 6. Examples for diagrams for which we perform a direct expansion of the two-loop tadpole subgraph, so as to avoid having to compute two-loop building blocks.
which allows us to re-write eq. (2.6) as We can thus compute the Lorentz and colour parts of diagrams involving this building block separately, and finally multiply them according to eq. (2.7). If a diagram contains two insertions of this building block (such as the diagram shown in figure 4) we introduce two independent sets of the δ symbols.
For diagrams for which the top quark dependence can not be completely constructed from these building blocks (i.e. diagrams which would require two-loop building blocks), we can generate diagrams with at least one (one-loop) building block vertex and treat the remaining top quark propagators explicitly with exp. This means that at most, we have to directly expand the subgraphs of four-loop forward-scattering diagrams. Two examples of such diagrams are shown in figure 6. For some diagrams of this type (such as the second example) exp also generates a one-loop subgraph which we discard, since this case is already included in the building block computation.

Real radiation: partial fraction decomposition and IBP reduction
After integrating out the top quark loops we end up with two-and three-loop phase-space integrals. It is useful to distinguish three different sets of integral families according to the number of loops and the cuts through two Higgs bosons and massless particles (see also figure 7 below): • Two-loop families with one or two three-particle cuts. These appear in the real corrections at NLO and real-virtual corrections involving three top quark loops at NNLO.
• Three-loop families with a three-particle cut. These appear in the remaining realvirtual corrections at NNLO.
• Three-loop families with one or two four-particle cuts. These produce the real-real corrections at NNLO.
In total 16 two-loop families, 28 three-loop families with a three-particle cut and 64 four-particle cut families are required to map all diagrams. We generate the four-point amplitudes with three independent external momenta, q 1 , q 2 and q 3 . After using the forward-scattering kinematics, i.e. q 3 = −q 2 , some of the propagators become linearly dependent; a partial fraction decomposition is needed to obtain proper input for the IBP reduction. For this purpose we have developed the program LIMIT [47] which automates this step of our calculation. This program is based on ideas presented in ref. [48].

JHEP01(2022)049
In the following we describe the workflow of LIMIT. First, LIMIT identifies linear relations among the propagators of a given family, derives relations to perform a partial fraction decomposition, and generates FORM routines to apply these relations to the integrals. Each initial, linearly dependent, family may yield multiple linearly independent families.
Next, families with multiple cuts are split into multiple families with one cut each, since each integral with two cuts corresponds to the sum of two integrals with one cut. Splitting the families in this way is necessary because LiteRed [49,50] can only handle families with a single cut. Now we use LiteRed to identify vanishing sectors of the families and to derive symmetry relations for the non-vanishing sectors. Based on this, LIMIT then generates FORM code to nullify vanishing integrals and apply the symmetry relations, reducing the number of terms present in each amplitude.
Finally, LIMIT identifies equivalent families based on their graph polynomials and groups them into classes. For each class, LiteRed is employed to derive relations for translating integrals of each family to integrals of a single representative family. As not all linearly independent families have the same number of propagators, LIMIT then tries to embed classes of families with fewer lines into classes of families with more lines. The rules for mapping integrals into this minimal number of representative families are translated into FORM code and applied to the amplitudes.
By applying LIMIT to the three sets of integral families mentioned above, we obtain 4, 5 and 14 integral families, respectively. Their graphical representation is shown in figure 7. For more details we refer to [47]. Using LiteRed we relate equivalent sectors of families in the three different classes and perform an IBP reduction to master integrals. In total O (500 000) integrals need to be reduced, which takes LiteRed less than two days. The two-loop reduction leads to 16 master integrals with a three-particle cut. They have already been studied in ref. [31] where the n 3 h contributions were computed where both Higgs bosons couple to different top-quark loops. We furthermore obtain 17 three-loop three-particle cut and 57 four-particle cut master integrals. They are shown in figure 8, as well as figures 9 and 10, respectively.
Details on the computation of the master integrals are provided in section 3. After inserting the results for the master integrals we obtain an expansion for each Feynman diagram in 1/m t and . For the real-radiation contribution we observe poles of order 1/ 4 .

Ultraviolet counterterms
The NNLO real-radiation amplitude receives a one-loop counterterm contribution induced by the NLO corrections, which are required to order . To be more precise, we have to renormalize α s , m t and the (external) gluon field. The corresponding renormalization constants are given by where µ is the renormalization scale and α s = α s (µ). We first perform the renormalization in the full theory (i.e., we have n f = 6 = n l + n h = 5 + 1). Afterwards, we perform a decoupling from α (6) s to α (5) s . Note that we renormalize the top quark mass first in the MS scheme; the transition to the on-shell scheme in the final result is straightforward.

Collinear counterterms
The NNLO collinear counterterms are obtained from the convolution of the LO and NLO cross sections σ (0) gg and σ (1) ij and the gluon or quark splitting functions. In general such a convolution has the form In eq. (2.10) the superscript "bare" refers to the renormalization of the parton distribution functions, i.e. the subtraction of the 1/ poles in connection to radiation from incoming partons. It does not refer to the renormalization of the ultraviolet divergences; here we assume that all ultraviolet counterterms are already taken into account. The kernels Γ ij in eq. (2.10) are obtained from the splitting functions as follows where P (0) ij and P (1) ij are one-and two-loop splitting functions which can be found in refs. [51][52][53][54][55]. For further details we refer to refs. [56,57] where the computation of the collinear counterterms is discussed for NNLO and N 3 LO single Higgs boson production. The combinatorics for single and double Higgs production is identical and thus the formalism outlined in detail in [56] can be applied in the context of this calculation. Additionally, this reference contains explicit formulae which solve eq. (2.10) for the renormalized cross section σ ij within perturbation theory. The resulting collinear counterterms are computed from convolution integrals of the following form, where the ellipses represent further contributions to the collinear NNLO counterterm σ (2) gg,coll . Note that in eq. (2.13) we distinguish the renormalization scale µ from the factorization scale µ f . We compute the collinear counterterms in the n l flavour theory, i.e., in eq. (2.12) the one-loop coefficient of the QCD beta function is given by Figure 11. Sample forward-scattering Feynman diagrams for the virtual corrections at LO, NLO and NNLO.
where n l is the number of massless quarks. We compute all contributions in terms of SU(N c ) colour factors C A = N c and C F = (N 2 c − 1)/(2N c ) and the trace normalization T f = 1/2. Inspection of eq. (2.12) shows that at NNLO the collinear counterterm starts at O(1/ 2 ) whereas the real and virtual corrections contain poles up to O(1/ 4 ). In appendix A we provide useful formulae for the analytic computation of the convolution integrals, such as the one shown in eq. (2.13). These are obtained by expanding in δ; we have computed all contributions up to order ρ 3 and δ 30 .

Virtual corrections
Virtual corrections only exist for the gg and qq channels. In the latter case they contribute for the first time at NNLO and are suppressed by 1/m 2 t at the level of the amplitude (and so by 1/m 4 t at the level of the cross section). We discuss them in appendix C. In this section we concentrate on the channel gg → HH. Sample Feynman diagrams for the forward-scattering kinematics are shown in figure 11. At NNLO we distinguish contributions which have two, three and four closed top quark loops and we denote the corresponding contributions by n 2 h , n 3 h and n 4 h respectively. Note that all n 3 h virtual corrections are already contained in ref. [31], even those with a top quark loop which has no coupling to the Higgs bosons (such as the diagram in figure 11(g)). The n 2 h terms develop poles up to 1/ 4 which cancel against those of the real-radiation corrections. The n 4 h contribution appears only in the virtual corrections and is therefore finite. The corresponding analytic results are presented in appendix D.
For the virtual corrections we do not use the optical theorem but compute the contribution to the cross section from the form factors obtained in ref. [58] in an expansion up to order 1/m 8 t . The triangle form factors are even available up to order 1/m 14 t . The starting point is the amplitude for the process gg → HH which has two independent tensor structures. Expressed in terms of the "box" and "triangle" form factors it is given by where G F is the Fermi constant, a and b are adjoint colour indices and the two Lorentz structures are given by Analytic results for the form factors F tri , F box1 , F box2 are available in the ancillary files of ref. [58]. The differential cross section is constructed in terms of the amplitude of eq. (2.15) according to where the factors 1/2, 1/N 2 A = 1/8 2 and 1/(2 − 2 ) 2 are due to the identical particles in the final state, and the gluon colour and spin averages, respectively. The factor f 2PS comes from the d-dimensional two-particle phase space and is given by We perform the integration to obtain the cross section after expanding in δ. If we were to integrate over t, the integration boundaries would be δ dependent. We therefore make the following change of variables: In terms of the new integration variable y the integration boundaries become [0, 1] and the δ dependence is isolated to the integrand. This allows us to expand in δ at the level of the integrand, and thus integrate each term of the expansion individually. In the supplementary material to this paper we provide the results in such a way that the (divergent) virtual corrections can be extracted. We perform the construction described above for the ultraviolet renormalized form factors as published in ref. [58]. The renormalization of the gluon wave function, top quark mass and strong coupling in the LO expression induces 1/ poles at NLO and 1/ 2 poles at NNLO. The insertion of the one-loop counterterms in the virtual NLO expression (which has an expansion starting at 1/ 2 ) induces 1/ 3 poles at NNLO.

Transformation to form
In this section, we discuss the computation of the 17 three-loop three-particle cut master integrals depicted in figure 8, as well as the 57 three-loop four-particle cut master integrals JHEP01(2022)049 depicted in figures 9 and 10. We follow the method of differential equations [59,60] and use LiteRed to take the derivative of each master integral w.r.t. x and reduce the resulting integrals back to master integrals. Thus, we obtain two systems of differential equations: one for the three-particle cut integrals and one for the four-particle cut integrals. These systems have the form L is the vector of master integrals, and the matrix M is block-triangular with entries rational in both x and . In the following we take two approaches to solve eq. (3.1): in section 3.2 we seek a canonical basis of master integrals [61] and in section 3.3 we compute them in terms of an asymptotic series in δ.

Finding a canonical basis
A canonical basis C is a particularly useful basis of integrals in which the differential equations take the form where the x i are constants and the matricesM i do not depend on or x. A system of the form of eq. (3.2) can be integrated order-by-order in in terms of Goncharov Polylogarithms (GPLs) [62] over the letters x i . To find a rational transformation from L to C we use the program Libra [63] and adapt Lee's algorithm [64] to the problem at hand. As the matrix M (x, ) contains rational functions with numerator and denominator degrees as high as 20 in our case, a direct application of Lee's algorithm would lead to unmanageable intermediate expression swell. We thus deviate from the standard algorithm by first bringing all diagonal blocks to Fuchsian form (that is, they have only simple poles in x), followed by bringing the off-diagonal blocks to Fuchsian form. We then normalize the eigenvalues of all residues to be proportional to or take the form a ± 1/2, where a is an integer.
Consequently, we arrive at an intermediate basis I, related to L by a rational transformation, in which the differential equations are in Fuchsian form and do not contain spurious poles: For the three-particle cut master integrals the poles x i are given by 0, 1/4 and 1, corresponding to the kinematic limits s → ∞, s → 4m 2 H and s → m 2 H , respectively. In the case of the four-particle cut master integrals there are also poles at −1/4 and −1, corresponding to the kinematic limits s → −4m 2 H and s → −m 2 H respectively. 1

JHEP01(2022)049
Rational transformations render the eigenvalues of all poles proportional to , apart from ±1/4. Therefore to find a canonical basis, we need to rationalize the square roots √ 1 − 4x and √ 1 + 4x. For the three-particle cut master integrals only the first square root plays a role. For this system, we perform the variable change mapping the physical interval x ∈ [1/4, 0) to z ∈ [1, 0). In the case of the four-particle cut master integrals we need to rationalize both square roots simultaneously, requiring the more complicated variable change . With all roots rationalized, we now proceed by normalizing the corresponding eigenvalues to be proportional to integers, bringing the off-diagonal blocks into Fuchsian form once more and factoring out to arrive at a canonical basis in z and k for the three-and four-particle cut systems, respectively. These resulting system can be integrated in terms of GPLs and their boundary conditions can be fixed in the limit x → 1/4.
Once we have the leading-order term in the δ expansion for the master integrals (see section 3.3), the exact expressions for the master integrals can be computed in two ways. The first is to transform the leading-δ terms of all master integrals into the canonical basis. The differential system can then be integrated order-by-order in and the integration constants are fixed by comparing the resulting expressions to the boundary conditions in the limit δ → 0.
The second option is to first compute the path-ordered exponential U of the canonical differential equation matrix M in terms of which we can write the master integrals of the canonical basis as where I 0 does not depend on x. Libra can decompose I 0 as where K depends only on and c is the minimal number of necessary boundary conditions. The entries of c are given by the leading-order terms in the δ expansion of a certain subset of original master integrals, determined by Libra to be sufficient to fix all boundary conditions. The main difference between these two approaches is the apparent number of necessary boundary conditions. For example using the first approach, we need the integrals I

JHEP01(2022)049
In our calculation we have used both methods for the three-loop three-particle cut master integrals. For the three-loop four-particle cut master integrals we have only used the first method since our non-trivial letters make the construction of the matrix U difficult.
In the following section we describe how we compute the boundary conditions. We will compute the leading-order term in the δ expansion for all master integrals; the integrals which are not necessary to integrate the differential system serve as a cross check.

Computation of the boundary conditions
In the previous subsection we have presented the construction of the exact results for the phase-space master integrals. The purpose of this subsection is to provide boundary conditions needed to compute the exact expressions. In the next subsection we briefly describe a formalism which allows us to compute the δ expansion of all master integrals to a high enough order in δ to be sufficient for practical applications. In fact, in the combination of the real-radiation contribution with the virtual corrections and the collinear counterterm it is necessary to use the δ-expanded results in order check the cancellation of the poles analytically.
We start with some definitions which are common to the three-and four-particle phase-space integrals. For the phase-space measures of the Higgs bosons we write and the phase-space measures of massless particles are given by (3.10) When necessary, we use a parametrization for the momenta such that the scalar products take the form where q i and p j are momenta of the initial and final state, respectively. In the following, we use the symbol "≈" to denote that the leading-order term in δ expansion agrees on the left-and right-hand sides of the equation, but that the sub-leading terms may differ.
For the (d − 1)-dimensional angular integration, we use the following property: (3.12)

Three-particle phase space
All 17 three-particle cut master integrals can be written as where Q i are one-loop integrals, which can easily be read off from figure 8. For convenience we provide explicit expressions in appendix E. Note that the Q (3) i do not depend on p 3 and p 4 , since the external momenta of the loop diagrams are q 1 , q 2 , p 5 and the combination p 3 + p 4 . The latter can be eliminated by using momentum conservation. By performing the integration over p 4 and taking the leading order in δ, we obtain There are six one-loop integrals appearing in the 17 master integrals. Three of them (L 1 , L 2 , L 3 ) are massless one-loop two-point functions and are known for a general dimension d: . (3.16) Note that in the soft limit there are only five independent one-loop integrals since L 2 ≈ L 1 . The remaining three integrals (L 4 , L 5 , L 6 ) are either triangle or box diagrams, and their soft limit can be computed using expansion by regions [65][66][67]. In particular, we use the implementation in the Mathematica package asy2.m [67]. Applying this method is straightforward but the details of the computation are rather technical; we do not discuss them here but simply refer to ref. [67]. Rather, in this paper, we discuss an alternative way to obtain the same results by making use of the Mellin-Barnes representation, which serves as a cross-check of the calculation based on expansion by regions. In the following, it is assumed that the integration contours of the Mellin-Barnes integrals are chosen properly, (see, e.g., ref. [68]) and we abbreviate the product of Γ-functions as

JHEP01(2022)049
We begin by expressing the integrals in the Mellin-Barnes representation: Note that the Mellin-Barnes representation of a given Feynman integral is not unique; here we show the representations used in our calculation. We choose them in such a way that the exponent of κ 5 is z 1 . From eq. (3.15), we see that the leading-order term in δ is obtained from the leading-order term in κ 5 . We solve z 1 integrals above by closing the integration contour in the right half plane and taking the residues of the poles of Γ-functions. For example, in L 5 only the poles of Γ(z 2 − z 1 ) contribute and their locations are z 1 = z 2 − n where n = 0, 1, 2, . . .. Thus the leading order in κ 5 corresponds to z 1 = z 2 pole. Keeping this in mind, we first solve the z 2 integral also by closing the integration contour in its right half plane. This time, the poles from Γ(−z 2 ) and Γ(−z 2 − − 1) contribute. Among these, the pole at z 2 = − − 1 gives the smallest value of z 2 , so we consider the residue at this pole. The resulting expression is now a one-dimensional integral over z 1 , and its leading-order term in κ 5 is obtained by taking the residue at z 1 (= z 2 ) = − − 1. In this way, we extract the leading-order term of L 5 in κ 5 . The calculation of L 4 and L 6 proceeds in a similar way. For L 4 , the two poles at z 2 = 0 and z 2 = − contribute to the leading order and we have to take both of them into account. For L 6 , the leading-order term is obtained from the pole at z 1 = z 2 − − 1 instead of z 1 = z 2 . The results are given by , .
We follow ref. [31] to perform the three-particle phase-space integration and obtain the following results for the leading-δ contribution of the master integrals: where for brevity for some integrals we have truncated the expansion at 1/ 2 and do not display the finite term, and we set µ 2 = s. Note that the µ dependence can easily be reconstructed by multiplying (µ 2 /s) 3 . The complete set of orders needed for our calculation and also additional terms in the δ expansion can be found in the supplementary material of this paper [69]. In order to fix the boundary conditions of the differential equations discussed in section 3.2, only the results for I 1 , . . . , I 6 and I 8 , . . . , I 11 are needed. Therefore in practice, only L 1 , L 2 , L 3 , L 5 and L 6 need to be computed. All other master integrals in eq. (3.20) and also the higher-order terms of the δ expansion serve as cross-check.

Four-particle phase space
The 57 four-particle phase-space master integrals can be parametrized as (3.21) where the soft limit all Q (4) i are given in appendix E. We group the Higgs momenta (p 3 , p 4 ) into p HH ≡ p 3 + p 4 and the parton momenta (p 5 , p 6 ) into p gg ≡ p 5 + p 6 by inserting where m 2 HH = p 2 HH and m 2 gg = p 2 gg . The momentum assignments of Q (4) i are chosen such that they do not depend on p 3 , p 4 or p HH in the soft limit, which is possible because of the fact that p 3 and p 4 appear only in the combination p HH , which can be eliminated using momentum conservation. The dependence of p 5 , p 6 and p gg is chosen such that Q (4) i factorizes in these variables, i.e., when a propagator contains p 5 + p 6 we replace it with p gg and in the remaining cases we keep p 5 and p 6 .
With the help of this parametrisation, it is straightforward to perform the phase-space integration over p 3 , p 4 , p 5 and p 6 using the well-known trick of boosting to the centre-of-mass frames of p HH and p gg . The resulting integrals can be evaluated in an expansion in δ. The leading-δ terms read As in eq. (3.20) the expansion of some integrals has been truncated before the finite term can be displayed, and we set µ 2 = s. We provide higher order and additional terms in the δ expansion in the ancillary files [69]. Of the 57 boundary conditions given in eq. (3.24) we only require the information from the following integrals, and a dedicated calculation is only necessary for one integral of each set: {I

Deep δ expansion of the master integrals
In this subsection, we briefly mention a method to obtain for the master integrals higher order terms in the δ expansion with the help of differential equations. This method is used also in refs. [15,31,70,71]. After changing from x to δ the differential equations in eq. (3.1) can be written as We are interested in L as a series expansion in δ and which has the form L = n 1 ,n 2 ,n 3 c n 1 ,n 2 ,n 3 δ n 1 (log δ) n 2 n 3 , (3.26) where the coefficients c n 1 ,n 2 ,n 3 are constants containing ζ 2 , ζ 3 , ζ 4 and log 2. The matrix M (δ, ) in eq. (3.25) consists of rational functions of δ and and it can be expanded in both variables and truncated at a certain order. By substituting the ansatz in eq. (3.26) into eq. (3.25) we obtain a set of recurrence relations for the coefficients c n 1 ,n 2 ,n 3 where the initial conditions are obtained by the leading terms calculated in the previous subsection. Thus, by solving the recurrence relations, it is possible to obtain the deep δ expansion of the master integrals efficiently.

Comparison of exact and δ-expanded phase-space master integrals
Following the approach described in the previous subsections we were able to obtain analytic results for all three-and four-particle cut master integrals, which are exact in s and m 2 H . However, the expressions are quite large (in some cases of the order of several MB) and contain GPLs which have to be evaluated at some fixed value of their argument. We did not make any effort to obtain minimal expressions, since for our application it is sufficient The panel in the top row correspond to three-particle cut master integrals and the corresponding pictures can be found in figure 8. The pictures for the four-particle cut master integrals in the bottom row can found in figure 9. The considered order is indicated on the y-axis labels.
to consider an expansion in the limit δ → 0. This can be obtained in a straightforward way by either expanding the exact results or by using the boundary conditions of section 3.3 in combination with the differential equations. In figure 12, for four typical examples of coefficients of master integrals, we compare the δ-expanded and exact results. We normalize the curves to the exact expressions and plot the expansions as functions of δ. In each panel three curves are shown, which correspond to approximations including 10, 30 and 50 terms of the δ expansion. All curves tend to 1 for δ → 0, so the plots focus on the region 0.8 < δ < 1 where the expansions start to diverge from the exact results.
The 0 term of the three-particle cut integral I (11) (the first panel of figure 12) is an example for which as few as 10 δ-expansion terms provide an excellent approximation. Even for δ = 1 the deviation from the exact result is at the level of 0.01%. Including additional expansion terms leads to further improvement. For the other three examples shown, we observe that the δ 10 curves show deviations from the exact result at the percent level, even for δ values as small as 0.8. The approximations including expansion terms to order δ 30 and δ 50 show percent-level agreement with the exact results for δ values up to around 0.9 and 0.95, respectively.
In our final expression for the cross sections, we use the expansions up to δ 30 ; from a practical point of view this is sufficient, since our large-m t approximation is only valid below the top threshold at s = 4m 2 t . This corresponds to δ = 1 − m 2 H /m 2 t ≈ 0.48, where the deviation of the expansions up to δ 30 from the exact results is below 0.0001% for all master integrals. On the other hand, the expansions up to δ 10 show deviations of several percent.

JHEP01(2022)049 4 Results
We are now in a position to combine the virtual, real-radiation and collinear counterterm corrections according to eq. (2.3). It is interesting to check their contributions to the cancellation of the poles. In the gg channel the highest-order poles, i.e. the 1/ 4 and 1/ 3 terms, receive contributions only from the virtual and real-radiation corrections. Starting from 1/ 2 the collinear counterterm also contributes. In fact, there are two sources which can be identified from eq. (2.12): either the term proportional to 1/ 2 is convoluted with the LO cross section, or there are two convolutions of the one-loop splitting functions in the O(α s ) term, again with the LO cross section. The convolutions involving the two-loop splitting function and the convolutions of the one-loop splitting function with NLO partonic cross sections develop only 1/ poles.
For the partonic channels with quarks or anti-quarks in the initial state there are no virtual corrections. 2 As a consequence the renormalized real-radiation contribution develops at most 1/ 2 poles which cancel against those of the collinear counterterm contributions.
In the following we present analytic results for the renormalized partonic NNLO cross sections for all five channels. We restrict the expressions to the two leading terms in the δ expansion and we present only results up to order 1/m 2 t (or, ρ 1 ). To obtain compact expressions we set µ 2 = m 2 H . The top quark mass is renormalized in the on-shell scheme. In the supplementary material to this paper [69] we provide expansions for general renormalization scale µ up to δ 30

JHEP01(2022)049
Note that the subset of the n 3 h real-radiation terms with a closed top quark loop which has no coupling to Higgs bosons only contribute to higher-order δ terms not shown in eq. (4.1). The n 4 h terms (see eq. (D.1) in appendix D) are not included in this expression. The result for the gq channel is given by where σ xy = ∆ xy + σ qq .
In the ancillary files [69] we provide the renormalized cross sections with tags for the virtual, real-radiation and collinear counterterm contributions, which allows one to extract the individual contributions if desired. This means that these expressions contain poles in , which cancel upon removing the tags.
In figure 13 we plot the partonic cross section for all five channels as a function of √ s, with the renormalization scale µ 2 = m 2 H . Furthermore we use m t = 173.21 GeV, m H = 125 GeV and α (5) s (m H ) = 0.1127. We combine all corrections, including the n 3 h terms from ref. [31], the n 4 h terms from appendix D and the qq → HH virtual corrections from appendix C. In the lower part of each plot we show the ratio to the highest ρ term, which is ρ 3 for the gg and ρ 4 for all other channels.
The overall picture is similar to the LO, NLO and NNLO n 3 h contributions [13,31]; a reasonable convergence pattern is observed for √ s 320 GeV, which rapidly deteriorates for values of √ s above the top quark threshold. For the qg, qq and qq channels one observes a better convergence behaviour than for the gg channel. Similar to NLO, the qq channel shows the worst behaviour. Note that there is a strong hierarchy in the numerical size of the various contributions, ranging from O(10 −2 ) fb for the gg channel to O(10 −5 ) fb for the qq and qq channels. It is interesting to note that below threshold the mass corrections in all cases significantly increase the cross section, as can be seen in the ratio plots.
By construction all curves have to approach zero in a continuous way for √ s → 2m H . This is clearly visible from the plots for all channels with quarks in the initial state, however Ratio to Ratio to Ratio to for the gg channel the plot gives the impression that this is not the case for the ρ 2 and ρ 3 curves. Inspecting the analytic result reveals that there are terms containing ρ 2 √ δ log 4 (δ) and ρ 3 √ δ log 4 (δ) which are responsible for a rapid rise in the cross section close to the threshold. Figure 14 shows the LO, NLO and NNLO contribution from the gg channel to the cross section as a function of √ s. At LO we use the exact expressions and at NLO and NNLO expansions up to ρ 4 and ρ 3 , respectively. The NLO corrections increase the cross section by about a factor two and the NNLO contributions add approximantly another 20%.
Although the radius of convergence of the large-m t expansion is quite limited here, it contains useful information for the construction of approximated NNLO results. For example, in ref. [18] NLO virtual corrections to double Higgs boson production are considered. The analytic structure of the form factors close to the top quark threshold is combined with results obtained in the large-m t limit. Stable approximations are only obtained after including power-suppressed 1/m t terms. This approach has been applied in refs. [72,73] to the Higgs-gluon form factor and also there it was found that it is important to include higher-order 1/m t terms to obtain a stable result in the high-energy region.

Conclusions
The main achievement of this paper is the computation of the NNLO real-radiation contribution to the total cross section of the process gg → HH in an expansion for large top quark mass. This improves the description beyond the infinite top quark mass limit and provides important information for NNLO approximation procedures.
We provide a detailed description of the various steps of our calculation and in particular discuss the computation of the three-and four-particle phase-space master integrals. We provide various intermediate results which might be useful for other calculations. In

JHEP01(2022)049
the ancillary files of this paper we provide separate results for the virtual corrections, real-radiation contribution and the collinear counterterm.
The inclusion of the power-suppressed terms significantly increases the cross section below the top quark threshold. For example, at NNLO, for a partonic center-of-mass energy √ s ≈ 300 GeV the cross section increases almost by a factor of two after including the 1/m 4 t and 1/m 6 t terms.
where Q(z) can be rational functions of z, a plus distribution or a delta function. It may also contain logarithms or dilogarithms. In the latter case it is convenient to rewrite the arguments using the identities In this way, we obtain Q(z) in terms of Li 2 (1 − z) and Li 2 (1 − z 2 ) and the expansion around δ = 0, which corresponds to z = 1, is straightforward.

A.1 σ (0)
gg Let us first consider the contributions involving the LO cross section which has a simple series expansion in δ where for our application we have N max = 2. Note that the coefficients c (n,u),(0) gg are available from an explicit calculation of σ (0) gg (x). In the following we discuss the various cases for Q(z) which appear.
For the case of Q(z) = z r , where in our case only r ≥ −1 is needed, we have and for Q(z) = 1/(1 + z),  4) which are used to replace z by ν. Afterwards we insert eq. (A.8), expand in δ and integrate each term individually. This approach is also used for the results which we present below. Due to the more complicated functions Q(z) some of the formulae are more involved. Although the formulae we present are exact, we apply an upper limit to the δ expansion when we use them in our calculations.

JHEP01(2022)049
Now we consider Q(z) = (log(1−z)) r 1−z + . Using the definition of the plus distribution we obtain In order to compute the finite integral on the r.h.s., we consider Q(z) = (1 − z) −1+η which leads to . (A.12) Note that this formula is exact in η. Using σ After combining both expressions we obtain where Here ψ(x) is the polygamma function and ψ (x) ≡ ψ (1) its derivative. The results for eq. (A.11) are obtained by expanding the l.h.s. of eq. (A.14) in η according to (A.16)

JHEP01(2022)049
Note that the r.h.s of eq. (A.14) is already expanded in η. The comparison of the η 0 and η 1 terms provides results for eq. (A.11) for r = 0 and r = 1, which we need for our calculations.
Next we treat cases such as Q(z) = log(z) log(1 − z) 1 − z , Q(z) = z 2 log (1 − z), . . . , (A. 17) i.e., Q(z) contains factors of z and 1 − z with positive or negative exponents, and/or logarithms with z or 1 − z as argument. Here it is convenient first to consider To obtain the result for Q(z) = z 2 log(1 − z) we set r 1 = 2 and r 2 = η. Afterwards we expand eq. (A. 19) in η and take the coefficient of the linear term, on both sides. In a similar way one can treat Q(z) = log(z) by setting r 1 = η and r 2 = 0.
In case Q(z) (from eq. (A.17)) contains a log(z) term a non-integer value for r 1 must be chosen.
Cases such as Q(z) = log(z) log(1 + z) 1 + z , z 2 log(1 + z), . . . (A.20) can be treated in a similar way. Here a convenient auxiliary function is given by where the desired result is again obtained by a suitable choice of r 1 and r 2 supplemented by proper expansions. The generic integral is given by Note that the expression on the r.h.s. is significantly more complex than eq. (A. 19). This is due to the fact that the factor (1 + z) r 2 introduces a new type of denominator whereas (1 − z) ∼ δ and thus no new structure is introduced.

JHEP01(2022)049
Having dealt with logarithms, we now turn to functions Q(z) which involve dilogarithms. For Q(z) = z r Li 2 (1 − z) we have , (A. 25) and for Q(z) = z r Li 2 (1 − z 2 ) the integral is given by In figure 16 we show the size of σ qq→HH compared to σ qq . The virtual corrections form a significant part of the total cross section, particularly near the production threshold. Therefore they should be included for a proper NNLO description of the qq channel.

JHEP01(2022)049
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.