Three-loop form factors for Higgs boson pair production in the large top mass limit

We consider the virtual corrections to Higgs boson pair production at next-to- next-to-leading order, in the large top quark mass limit. We compute five expansion terms for the box-type form factors and eight expansion terms for the triangle form factor, which serve as useful input for the construction of approximations. We present analytic results for the form factors in the soft-virtual approximation. From a technical point of view the calculation is quite challenging since huge intermediate expressions are produced. We describe our methods and optimizations to overcome these difficulties, which might be useful for other calculations.


Introduction
The simultaneous production of two Higgs bosons is a promising process to obtain information about the self-coupling of the Higgs boson and thus the structure of the scalar potential.Although it is experimentally very challenging it is expected that this process can be observed after the high-luminosity upgrade of the LHC.
On the theoretical side there has been quite some effort to obtain precise predictions for differential and total cross sections for Higgs boson pair production.In analogy to single Higgs production the numerically most important contribution is provided by gluon fusion, followed by vector boson fusion, associated production with top quarks and the Higgs-strahlung process (see, e.g., Ref. [1]).
Exact leading order (LO) results for gg → HH have been available for more than thirty years [2,3].Next-to-leading order (NLO) corrections have been computed numerically much more recently, and are available from two independent groups [4][5][6].Note that the numerical evaluations are quite expensive.For this reason it is important to have approximations at hand, which are valid in certain regions of the phase space.Among them are large top quark mass expansions [7][8][9] which are available up to order 1/m 12  t [8].Furthermore, in Ref. [10] an expansion around small transverse momentum has been performed and results in the high-energy region are available from [11,12].They have been combined in Ref. [13] with the exact calculation from [4,5] to provide a precise grid for the NLO virtual corrections [14].In Ref. [15] exact results for the real radiation contribution have been combined with the effective-theory virtual corrections.Interesting approximations for gg → HH at NLO have been constructed in Ref. [16] where expansion terms from various regions have been combined with the help of a conformal mapping and Padé approximation.The same method has been been applied in Ref. [17] (using the triangle form factor results of this paper) to the Higgs-gluon form factor, an important ingredient of single-Higgs boson production, in order to reconstruct the full quark mass dependence. 1t NNLO exact results are currently out of range, which makes it even more important to obtain approximations, if possible from various kinematic regions.Within the effective theory, where the top quark mass is assumed to be infinitely heavy, NNLO corrections have been computed in Refs.[8,19,20].Power-suppressed terms have been obtained in Ref. [21], where the soft-virtual approximation was constructed.Real corrections which originate from three closed top quark loops have been computed in Ref. [22].In Ref. [23] approximate NNLO expressions are constructed on the basis of the exact NLO results [5] and further NNLO building blocks which are also available for finite top quark mass.Other NNLO contributions, such as the three-loop virtual corrections, are taken in the infinite top quark mass limit.The results of this paper provide additional 1/m 2 t corrections to the three-loop gg → HH amplitude which could improve the approximations of Ref. [23].
The resummation of threshold-enhanced logarithms to next-to-next-to-leading logarithmic (NNLL) accuracy has been performed in Refs.[24,25] and differential distributions up to NNLO for various observables were computed in Ref. [26] in the heavy-top limit.More recently, finite top quark mass effects have also been included [27].
At N 3 LO first results are available in the limit of infinitely heavy top quarks.In Ref. [28] massless two-loop box contributions have been computed and four-loop corrections to the effective coupling of two Higgs bosons and two, three or four gluons became available from [29,30].
In this work we consider NNLO virtual corrections to gg → HH and compute the three relevant form factors for a large top quark mass.We evaluate five expansion terms for the box-type form factors and eight expansion terms for the triangle form factor, i.e., up to order 1/m 8 t and 1/m 14 t , respectively.The results for the two box-type form factors are new.The results for the triangle form factor have been obtained in Ref. [31,32] up to order 1/m 8 t , the higher-order expansion terms presented here are new.In a previous work [21] expansion terms up to 1/m 4 t were computed for the (differential) cross section, but not for the form factors.Our results constitute important input for the construction of approximations.For example, it is possible to extend the consideration of Ref. [16] to NNLO.Furthermore, as already mentioned above, it might be possible to improve the approximations of Ref. [23].
The remainder of the paper is organized as follows: in the next section we introduce our notation and define the form factors.We provide technical details in Section 3, and mention several optimizations which were crucial to be able to perform the calculations.Ultraviolet renormalization and infrared subtraction are discussed in Section 4 and both analytical and numerical results are shown in Section 5. We conclude in Section 6.

Setup
The amplitude for the process g(q 1 )g(q 2 ) → H(q 3 )H(q 4 ) is conveniently decomposed into three form factors.In the following we outline their precise definition.We start with the amplitude which is given by where a and b are adjoint colour indices and the two Lorentz structures are given by with M 1 and M 2 can be projected from M µν using the projectors where ǫ = (4 − d)/2 is the standard dimensional regularization parameter.
The Feynman diagrams involving the triple-Higgs boson coupling only contribute to A µν 1 , which is the only structure relevant for single-Higgs production, therefore it is convenient to decompose M 1 and M 2 into "triangle" and "box" form factors with the prefactor where T = 1/2 and n h = 1 have been introduced for convenience.We furthermore define the expansion of the form factors in α s as with X ∈ {tri, box1, box2}.Note that F (i) X corresponds to the (i + 1)-loop result.In our final expressions the strong coupling constant is defined with five active quark flavours, which is an appropriate choice since we consider the top quark mass to be large.
In the course of the calculation it is convenient to introduce the Mandelstam variables with T n h It is furthermore convenient to express the final result in terms of the transverse momentum of one of the Higgs bosons which is given in terms of the Mandelstam variables by (equivalent to Eq. ( 3))

Calculation details
We generate the Feynman amplitudes with the help of qgraf [33] and obtain 11,197 and 5703 diagrams at one, two and three loops.Note that both one-particle irreducible (1PI) and one-particle reducible (1PR) contributions have to be considered.Sample diagrams are shown in Fig. 1 together with the corresponding colour factors expressed in terms of the Casimir invariants of SU(N c ): ). Furthermore we have T = 1/2 and use the labels n l and n h for closed massless and massive fermion loops respectively.For numerical evaluation we set n l = 5 and n h = 1.In the following subsections we provide several technical details of the calculation of the form factors.

Asymptotic expansion
The programs q2e and exp [34][35][36] have been designed to work hand-in-hand when applying a (possibly nested) asymptotic expansion involving a large external momentum or a large internal mass to an amplitude generated by qgraf [33].The output of exp is FORM [37] code2 for each sub-diagram which has to be considered according to the rules of asymptotic expansion (see, e.g., Ref. [38]).
In this case we apply the rules of asymptotic expansion for the limit where q 2 1 = q 2 2 = 0 are the incoming gluon momenta and q 2 3 = m 2 H . Equation ( 11) implies that the Feynman amplitudes are expanded in powers of possibly multiplied by logarithms of these ratios.
The main purpose of Eq. ( 11) is the reduction of the number of scales in the loop integrals.Furthermore, the three-loop integrals are factorized into products of lower-loop integrals.
In the box diagrams we initially have the scales s, t, m 2 H and m 2 t and in the triangle diagrams s and m 2 t .After asymptotic expansion we find the following products of integrals Type of integrals for scales hard subgraph co-subgraph 3-loop vacuum - Note that integrals with more than one scale only have to be considered at one-loop order; the corresponding integral families are well-studied in the literature [39][40][41].In the above table "massless" refers to the propagator masses only.Dependence on the Higgs boson mass is retained.In the one-loop massless box case, degenerate cases also occur for which one of the scales is absent.
In cases in which one has to deal with products of integrals we organize the output of exp in such a way that we perform the vacuum integrals first, since it is simpler to compute vacuum tensor integrals than tensor integrals for families with external momenta.In fact, at one and two loops vacuum tensor integrals with arbitrary rank can be treated. 3or three-loop vacuum integrals we implement projectors which are discussed in detail in Figure 2: Sample three-loop diagrams (left) and the corresponding co-subgraphs (right) which result from the asymptotic expansion according to Eq. ( 11).The blobs represent effective vertices from the hard subgraphs.They correspond to vacuum integrals.Subsection 3.2.For the remaining massless integrations no tensor integrals have to be solved.
The vacuum integrals are performed with the FORM package MATAD [43] and for the massless integral families we use FIRE [44] to obtain integral tables which express all integrals appearing in the amplitudes in terms of master integrals (see Fig. 7 of Ref. [21] for graphical representations).Analytic expressions for the latter can be found in Refs.[21,[39][40][41].
Let us illustrate the procedure described above using two typical Feynman diagrams shown in Fig. 2. We show the three-loop diagrams which have to be expanded in all external momenta, and the corresponding lower-loop co-subgraphs which appear after applying the scale hierarchy of Eq. ( 11).The blobs represent effective vertices from the expanded hard subgraphs which we do not show explicitly.
Note that due to the rules of asymptotic expansion the hard subgraphs have to expanded in all small quantities, which in this case are the external momenta q i but also the loop momenta of the co-subgraphs.This results in a multi-dimensional Taylor expansion which we want to compute up to 5 th order (i.e. up to order 1/(m 2 t ) 4 ) for the box form factors and up to 8 th order (1/(m 2 t ) 7 ) for the triangle form factor.At this point the intermediate expressions can become huge and special measures and optimizations have to be applied in order to obtain the results with the computing resources available.These methods are described in the following subsection.

Projection
A major bottleneck in the computation of [21] is the calculation of three-loop tensor vacuum integrals.After expansion in 1/m 2 t the intermediate expressions become rather large, which cause these routines to perform very poorly.In order to avoid this issue in this work we project the sub-diagrams which contain a three-loop vacuum integral onto a suitable ansatz, and compute individual terms of this ansatz separately.The intermediate expressions for each term become much smaller, and we no longer have to compute tensor integrals.The diagrams contributing to the triangle form factor have a simpler structure, and thus use a simplified version of the method discussed below.For this reason we are able to compute an additional three terms in the expansion, compared to the depth of the box-type form factors.We elaborate on this at the end of the subsection.
Each diagram can be written in the following way, (see also [45], here we extend the ansatz to account for the additional external momentum), where L max depends on the depth of the 1/m 2 t expansion being considered.Since we consider the process g(q 1 )g(q 2 ) → H(q 3 )H(q 4 ) we have that q 2 1 = q 2 2 = 0; we can therefore set i = j = 0 in the ansatz here.Associated with the six possible scalar products between the momenta are six derivative operators with which one can construct projection operators P i,j,k,l,m,n to project particular coefficients C i,j,k,l,m,n of the ansatz A from the amplitude, i.e.
We compute all necessary derivative operators applied to the diagrams after the expansion in 1/m 2 t , before we perform the three-loop vacuum integral procedures.Each derivative operator (that is, each i,j,k,l,m,n required) is applied as a separate task and all results are combined at the end.This ensures that intermediate expressions remain a manageable size, and that no derivative operator is computed more than once.
For reasonable performance it is crucial to not repeat the 1/m 2 t expansion of the diagrams for each of the above tasks, since it is a very computationally expensive procedure.The expansion is performed just once; the intermediate result is then split into parts containing particular numbers of each external momentum and stored.The projection tasks can load exactly the part which will yield a non-zero result after taking the derivatives with respect to the external momenta.
The structure of the computation is summarized below.For some aspects we provide a more detailed description in Section 3.3.
1. 1/m 2 t expansion: (a) Sum all diagrams with the same colour factor to make "super-diagrams".Many terms are common to multiple diagrams, so summing them reduces the total size of the intermediate expressions.At three loops there are 5703 Feynman diagrams which form nine super-diagrams with the colour factors (considering only three-loop vacuum sub-diagrams), The super-diagrams with colour factors proportional to d abc d abc , which arise from Feynman diagrams with two closed fermion loops with three gluon couplings each, are found to vanish after expansion in 1/m 2 t in Step 1.(d) (see below), which is why this colour structure is not listed in Fig. 1.Note that of the eight three-loop colour structures listed in Fig. 1 only (T n h ) 3 has no 1PI three-loop vacuum contribution.
This produces 5 × 9 = 45 projected super-diagrams, to be expanded in 1/m 2 t .Apply Feynman rules and perform Dirac algebra.The coefficients of these Lorentz structures (in the round brackets of Eq. ( 4)) will be included when everything is combined in Step 2. (b).
(c) Use graph symmetries to reduce the number of terms and size of expressions.
(d) Perform the 1/m 2 t expansions.These are heavy computations, for which we use computing nodes with relatively large amounts of memory and processing cores (at least 96GB memory and 12 cores).It is crucial to not duplicate any work here; we make extensive use of the FORM statements Collect (to reduce the number of terms to be processed) and ArgToExtraSymbol (to temporarily reduce the size of the expressions).After expansion, graph symmetries are again used to reduce the number of terms and size of the expressions.The five most difficult projected super-diagrams are those with colour factor C 2 A T n h .To expand to 1/m 8 t these each require a wall time of around 10 days.The total size of the (gzip compressed) stored expressions for the expansions of the 45 projected super-diagrams is 324GB.

Projection:
(a) For each of the necessary operators (see Eq. ( 17)), load the relevant part of the expanded super-diagram (for the example of Eq. ( 17), the part containing terms with one q 1 , one q 2 and two q 3 ).All other parts would yield zero after differentiation, so do not need to be loaded.The differentiation must be performed inside FORM CFunction arguments to avoid an enormous blow-up of intermediate expression sizes.These tasks are much easier, computationally, than those of the expansion steps.They are run requiring only 8GB of memory and 4 processing cores.To obtain the 1/m {0,2,4,6,8} t terms of the expansion there are {15, 38, 88, 174, 324} derivatives to compute for each of the 45 projected super-diagrams, yielding {675, 1,710, 3,960, 7,830, 14,580} tasks to be run respectively.These tasks required a total time of approximately 1,600 days to complete; running tasks concurrently on our cluster this corresponds to a wall time of about 1 month.
(b) The results of these operators applied to the diagrams allow one to construct the result in the form of the ansatz Eq. ( 13).Combining all terms, along with the coefficients of the Lorentz structures of Eq. ( 4), yields the final result for the form factors M 1 and M 2 .
As mentioned above, some simplifications are possible when computing the triangle form factor.It comes only with the Lorentz structures g µν and q 1,ν q 2,µ (see Eq. ( 2)), thus in step 1.(b) fewer projected super-diagrams are produced since we can ignore the additional three structures required by the box-type form factors.The ansatz of Eq. ( 13) can also be simplified; only the index l needs to be non-zero.Thus, fewer derivative operators need to be computed in step 2. (a): for the 1/m {0,2,4,6,8,10,12,14} t terms of the expansion we must apply just {2, 2, 3, 3, 4, 4, 5, 5} derivative operators.

Calculation optimizations
In this Section we outline a few methods by which we were able to optimize the computation, in addition to the projection procedure described above.We note that without such optimizations, computing the expansion to a depth 1/m 8 t (and likely 1/m 6 t ) for the box-type and 1/m 14 t for the triangle form factors would not have been possible with the computing resources available to us.

Graph Symmetries
In Step 1. (c) and 1.(d) we use graph symmetries to reduce the size of the intermediate expressions.We map the vacuum integrals to a minimal set by using rotation and reflection symmetries, implemented by re-labelling the line momenta of equivalent graphs such that they coincide.Some examples of this procedure are shown in Table 1.
After expansion in 1/m 2 t many integrals appear with higher-power ("dotted") propagators.One can move the dots around the graph, using the same symmetry relations as described above, to obtain a smaller set of integrals.

ArgToExtraSymbol
In step 1. (d), the 1/m 2 t expansions are performed.At this point, the FORM representation of the terms in the expressions looks something like + Den(l1,mt) * Den(l1+q1,mt) * ... * Den(l2-q3,mt) * ( many terms ) where the Den functions represent the propagators to be expanded; they are of the form 1/(m 2 t − (l 1 + q 1 ) 2 ), for example.The "many terms" inside the brackets are coefficients which do not take part in the expansion.Since there can be many thousand such coefficients, it is crucial to keep them bracketed away during the multi-module expansion routine, to keep the number of terms small and avoid expanding the same product of Den functions many times.One typically achieves this with a construction such as Bracket Den; .sortCFunction f; Collect f; which moves the bracketed terms inside the argument of f.While this does indeed keep the number of terms small, it does not (greatly) reduce the size of the expression.If the expression is large enough to require disk-based sorting at the end of each module of Top-level Topology Graph 1 Graph 2 Relabelling  the expansion routine, one still has a severe performance bottleneck.We resolve this by additionally making use of the statement ArgToExtraSymbol f; after Collect f;, which replaces the (large) arguments of the fs with unique symbols, whose definitions are stored by FORM.More memory is required to store these definitions, but the resulting reduction in size of the expression involved in disk-based sorting provides a large speed-up of the expansion routine.After expansion is complete the original coefficients may be recovered with the FromPolynomial statement.
Let us remark that the use of ArgToExtraSymbol is also essential to make possible and speed up the calculation of the subdiagrams where two-loop vacuum tensor integrals are needed.

Compression
In step 1.(d) it was stated that the intermediate results of the 1/m 2 t expansions are compressed with gzip and stored, for use in step 2. (a).Since these compressed results occupy 324GB, they cannot easily be stored uncompressed on the storage available to us.In step 2. (a) the tasks can easily retrieve the relevant compressed intermediate result from network storage and decompress it onto the local storage of the node on which they are running, by making use of FORM's #system preprocessor command: As well as reducing the capacity required for the storage of the intermediate results, this also provides a large performance improvement by reducing the I/O load of the network and storage hardware when hundreds of tasks are running concurrently.

Renormalization and infrared subtraction 4.1 Ultraviolet divergences
The renormalization of the ultraviolet (UV) divergences is straightforward: • The top quark mass (m t ) renormalization is needed up to two loops.We first renormalize m t in the MS scheme, and then transform m t to the on-shell scheme.Note that higher order ǫ terms are needed in the corresponding one-loop expression since the NLO (two-loop) amplitude develops 1/ǫ 2 poles, even after all UV counterterms are taken into account.Since the LO (one-loop) amplitude is finite the twoloop term in the MS-on-shell conversion formula is only needed up to O(ǫ 0 ).
• The gluon wave function renormalization is also needed up to two loops.Note that, since we work in dimensional regularization, where scaleless integrals are set to zero, only diagrams with virtual top quarks contribute.These two-point functions have to be computed for q 2 = 0 which corresponds to on-shell gluons.Note that the gluon wave function renormalization agrees with the decoupling constant of the gluon field needed to relate five-and six-flavour QCD [46].
• The strong coupling constant renormalization up to two loops is performed in full six-flavour theory.
• Finally the decoupling relation for α s is needed up to two loops in order to express α s in terms of α s .Similar to the MS-on-shell mass relation also here the one-loop expression is needed up to order ǫ 2 .
The final result is expressed in terms of the top quark pole mass, and the five-flavour strong coupling, α s .It still contains poles up to order 1/ǫ 4 which are of infrared nature.They will be treated in the next subsection.

Subtraction of infrared divergences
For the subtraction of the infrared (IR) poles we follow Ref. [47], see also Refs.[8,48].For convenience we provide explicit expressions for the subtraction terms.We apply the IR subtraction at the amplitude level since we want the obtain finite expressions for the form factors.
After UV renormalization we have the following colour factors at one-, two-and three-loop order: In the following discussion we omit the overall factor T n h which is contained in the quantity X 0 , see Eq. ( 6).Note that the structures T 2 n h n l , T 3 n h n 2 l and T 3 n 2 h n l are not present in the two-and three-loop diagrams (cf.Fig. 1) but only arise from UV counterterms and IR subtraction (see below).
After UV renormalization, at two-loop order the colour factors {C A , T n l } come with 1/ǫ poles and C A also has a 1/ǫ 2 pole.At three-loop order, highest-order pole appearing with each colour factor is summarized in the following table, We have checked that all these poles cancel after applying the following IR subtraction procedure: finite form factors, F fin X , at NLO and NNLO are obtained via X − X − X − where, as in Eq. ( 7), X ∈ {tri, box1, box2}.The quantities on the r.h.s. of Eq. ( 21) are UV-renormalized.I (2) g on the r.h.s. of Eq. ( 21) are given by [47,48] with 5 Results In the following we discuss the results for the finite form factors constructed according to the prescription of the previous section.Note that the one-loop form factors have no dependence on the renormalization scale µ.At two and three loops the coefficients of the log(µ) terms depend on the choice of the IR subtraction terms.In our case it is convenient to cast the results for the two-and three-loop form factors in the following form where (µ 2 = −s) with β 0 as defined in Eq. ( 24) and The one-and two-loop results are expanded up to order 1/m 14 t , the three-loop expressions up to 1/m 8 t (box) and 1/m 14 t (triangle).For illustration we show the analytic result for the leading term (m 0 t ) of the three-loop for factors.The corresponding one-and two-loop results can be found in Ref. [12] and the triangle form factor up to 1/m 12 t with numerical values for the colour factors can be found in Ref. [17].Our results read where T = 1/2 has been chosen and the overall factor n h is contained in Eq. ( 6).Furthermore, we have introduced We refrain from showing explicit results for F (2) box2 which has a similar structure to F (2) box1 .Note that for most colour structures F (2) box2 starts at order 1/m 2 t except for the three colour structures C F n h , C A n h and n l n h which arise from (1PR and 1PI) diagrams with two closed top quark loops.The analytic results for the one-and two-loop box and triangle form factors expanded up to 1/m 12 t and 1/m 14 t respectively, the three-loop box form factors F (2) box1 and F (2) box2 expanded up to 1/m 8 t , and the three-loop triangle form factor F (2) tri expanded up to up to 1/m 14 t can be found in the supplementary file of this paper [49].Note that at two-loop order the 1PI (colour structures C A T n h and C F T n h ) and 1PR ((T n h ) 2 ) contributions are separately finite.At three-loop order this is not the case and the whole contribution has to be considered in order to arrive at finite form factors, see also discussion in Refs.[8,30].
Let us now briefly discuss the numerical impact of our calculation.For the numerical evaluation we use m t = 173 GeV and m H = 125 GeV and for the transverse momentum we introduce the parameter with r p T = 0.01 as default value.Furthermore we choose for the parameters introduced for closed fermion loops n l = 5 and n h = 1.
In Fig. 3 we show real parts of the one-and two-loop results for the three form factors as a function of √ s.We include terms up to order 1/m 14 t for the triangle and order 1/m 12 t for the box form factors. Lines with longer dashes include more expansion terms.Below the corrections show a very similar behaviour.Since the two-loop terms have proven to provide useful and important input into the Padé procedure [16], which can be used to obtain approximations valid in the whole √ s region, we expect that the three-loop terms are of similar importance.
For some applications it is advantageous to rescale the higher order corrections by the exact leading order contributions using where F (n),exp X and F (0),exp X are expanded up to the same order in 1/m t .We refrain from showing the corresponding results but simply want to mention that the differences between the F (n),exp X and the rescaled expression (31) become smaller with increasing order in 1/m t .In fact the curves which correspond to the deepest expansions are very close to each other.

Conclusions
We compute three-loop corrections to the process gg → HH in the large-m t limit and provide results for five expansion terms (up to order 1/m 8 t ) for the two box-type form factors and for eight expansion terms (up to order 1/m 14 t ) for the triangle form factor.
As compared to previous work [21] we have computed two 4 more expansion terms, which required significant reorganization and optimization of the calculations since huge expressions are obtained at various intermediate stages.We discuss these modifications in Section 3. Furthermore in Ref. [21] only partonic cross sections, rather than individual form factors, are available.
The analytic results for the form factors, which are provided in a computer-readable form in supplementary material [49], are useful input for the construction of approximations for gg → HH at NNLO, both for total cross sections and differential distributions.This concerns both the construction of Padé approximants along the lines of [16] (indeed these new results for the triangle form factor have already been used in [17]) but also approximation procedures which have been employed in Ref. [23].

Figure 1 :
Figure 1: Sample Feynman diagrams contributing to gg → HH.For simplicity we show diagrams with a triple-Higgs boson coupling only at one-loop order.A sample colour factor is shown below each diagram.However, note that in general a diagram contributes to more than one colour structure.Solid, dashed and curly lines denote quarks, Higgs bosons and gluons respectively.

Figure 4 :
Figure 4: Real parts of three-loop form factors as a function of √ s for r p T = 0.1.

Table 1 :
Graphs 1 and 2 are derived from the Top-level Topology, with different lines missing.This yields different, but equivalent, graphs.Line momenta are relabelled to make this equivalence manifest; we show the replacements required to map Graph 1 onto Graph 2. The arrows denote the direction of momentum flow.