Higgs boson contribution to the leading two-loop Yukawa corrections to $gg\to HH$

We analytically compute two-loop Yukawa corrections to Higgs boson pair production in the high-energy limit. Such corrections are generated by an exchange of a Higgs boson between the virtual top quark lines. We propose two approaches to obtain expansions of the massive two-loop box integrals and show that precise results are obtained for transverse momenta of the Higgs bosons above about 150 GeV. We discuss in detail the computation of all 140 master integrals and present analytic results.


Introduction
Higgs boson pair production is a promising process which can provide experimental information about the Higgs boson self coupling (see, e.g., Ref. [1]). It is thus important to provide precise theoretical predictions of this process. The dominant contribution to Higgs boson pair production comes from gluon fusion, mediated by a top quark loop. There are a number of works in the literature in which QCD corrections to gg → HH have been considered. The NLO QCD corrections are known exactly [2][3][4], however, the numerical approach is quite computationally demanding. In practice it is therefore advantageous to construct approximations based on several expansions, valid in different regions of phase space [5][6][7][8][9][10][11][12][13]. A subsequent combination of the numerical approach with these expansions leads to fast and precise results which cover the whole phase space [14,15]. At NNLO [16][17][18][19][20][21] and N 3 LO [22][23][24][25][26] only the large-m t expansion has been considered. The to date most precise predictions have been obtained in Ref. [27] where a NNLO approximation has been constructed, based partly on exact and partly on large-m t results.
Electroweak corrections are expected to be of the order of a few percent and thus they should be included in the theoretical description. In the Standard Model there are several couplings (gauge, Yukawa, Higgs boson self coupling) which are of different nature and can be treated separately. In this paper we take a first step towards the electroweak corrections and compute top quark Yukawa corrections originating from Higgs boson exchange in the top quark loop. More precisely, we consider diagrams like the one shown in Fig. 1. For this subclass only planar diagrams contribute and thus only planar integral families have to be considered.
Note that in the R ξ gauge there are also other Yukawa corrections from the exchange of neutral and charged Goldstone bosons. They are not considered in this paper. Rather we concentrate on corrections with a virtual Higgs boson.
In the case of QCD corrections the top quark is the only massive particle in the loop. As additional scales, one has the Mandelstam variables s and t and the Higgs boson mass from the final-state particles. Electroweak corrections introduce additional masses in the propagators of the loop integrals, which increases the complexity significantly.
There are further classes of diagrams with a Higgs boson exchange. In contrast to the diagram in Fig. 1 they either involve Higgs boson self couplings (see Fig. 2(a)-(d)) or are one-particle reducible (see Fig. 2(e)-(h)). The results for the triangle diagrams in Fig. 2(a) can be obtained from the integral families discussed in this paper. Note that the diagram classes (b), (c) and (d) also involve non-planar contributions. Diagrams (e)-(h) are one-particle reducible and factorize into a product of one-loop integrals.
The master integrals which are computed in this paper are sufficient to compute the contributions from Figs. 2(a), (e), (f), (g) and (h). However, in this paper we concentrate on the two-loop box contribution of Fig. 1 and pursue the following goals: • Develop a method to obtain high-energy approximations of two-loop four-point in- tegrals where two different masses are present inside the loops.
• Provide details of the analytic computation of the master integrals which appear in the subclass of diagrams considered in this paper.
• Provide explicit analytic results for the master integrals in the high-energy limit.
The remainder of the paper is organized as follows: in the next section we introduce our notation and in Section 3 we outline the expansions which we apply to the Feynman diagrams. In Section 4 details of the computation of the amplitudes in terms of master integrals are provided. In Section 5 we provide a detailed description of the computation of the master integrals and numerical results of the form factors are are given in Section 6. We conclude in Section 7. In the appendix we present results for three-dimensional Mellin-Barnes integrals which enter our result.

Notation
The Mandelstam variables for the amplitude g(q 1 )g(q 2 ) → H(q 3 )H(q 4 ), with all momenta (q i ) defined to be incoming, are given by with It is convenient to introduce the scattering angle θ and the transverse momentum of the Higgs bosons in the center-of-mass frame, which are given by Due to Lorentz and gauge invariance it is possible to define two scalar matrix elements M 1 and M 2 as M ab = ε 1,µ ε 2,ν M µν,ab = ε 1,µ ε 2,ν δ ab (M 1 A µν where a and b are adjoint colour indices and the two Lorentz structures are given by A µν 2 = g µν + 1 p 2 T q 12 (q 33 q ν 1 q µ 2 − 2q 23 q ν 1 q µ 3 − 2q 13 q ν 3 q µ 2 + 2q 12 q µ 3 q ν 3 ) , with q ij = q i · q j . The Feynman diagrams involving the Higgs boson self coupling only contribute to A µν 1 and thus, it is convenient to decompose M 1 and M 2 into "triangle" and "box" form factors with where T = 1/2, µ is the renormalization scale and G F is the Fermi constant.
We define the perturbative expansion of the form factors as where α t is given by α is the fine structure constant and s 2 W ≡ sin 2 θ W is the square of the sine of the weak mixing angle. Throughout this paper the strong coupling constant is defined with six active quark flavours. Note that the form factors are defined such that the one-loop colour factor T is contained in the prefactor X 0 .
In this paper we only consider the contribution of the diagram class shown in Fig. 1 to F box1 and F box2 .

Asymptotic expansion
For the computation of the two-loop integrals we follow two approaches, which we describe in the following. For this purpose it is convenient to distinguish the mass of the finalstate Higgs bosons (m ext H ) from that of the Higgs boson which propagates in the loops (m int H ). This means that for the process gg → HH we have the following dimensionful quantities: the Mandelstam variables s and t, and the masses m t , m int H and m ext H . In our two approaches we assume the following hierarchies: In approach (A) we treat the inequality m 2 t (m int H ) 2 at the level of the integrand by applying the hard-mass expansion procedure as implemented in the program exp [28,29]. For each Feynman diagram this leads to two subgraphs: the two-loop diagram itself and the one-loop diagram which contains all top quark lines. In the latter case the co-subgraph consists only of the Higgs boson propagator.
The two-loop subgraph is Taylor-expanded in m int H whereas the one-loop subgraph is expanded in the loop momentum of the co-subgraph, which is a one-loop vacuum integral with mass scale m int H . In addition, each subgraph is then expanded in m ext H , which is performed at the level of scalar integrals with the help of LiteRed [30,31].
At this point one has to deal with one-and two-loop four-point integrals which only depend on the variables s, t and m t . These integrals belong to the same set of topologies used in the calculation of the QCD corrections presented in Refs. [8,9]; we are able to re-use those results here.
where δ = 1 − (m int H ) 2 /m 2 t , and expand D h (p) in the limit δ → 0 at the level of the integrand. The expansion in m ext H is then performed in the same way as for approach (A), described above. The remaining integrals are two-loop four-point integrals with massless legs, where all internal propagators have the mass m t ; this is a different set of integral topologies to those of the QCD corrections and approach (A).
In the final result, it is advantageous to introduce δ = 1 − m int H /m t . By making the replacement we obtain an expansion in δ which often has better convergence properties than the expansion in δ (see also discussion at the end of Section 2 in Ref. [32]).

gg → HH amplitude and form factors
In this section we provide some details regarding how the two expansion approaches discussed in Section 3 are implemented. We generate the amplitude with qgraf [33] and process the output with q2e and exp [28,29] in order to generate FORM [34] code for the amplitudes. This yields 6 one-loop diagrams and 60 two-loop diagrams.
As mentioned above, in approach (A) exp identifies a one-and a two-loop sub-graph for each of the two-loop diagrams. The corresponding four-point integrals are expanded in m ext H using LiteRed [30,31] and then integration-by-parts (IBP) reduced to a set of master integrals using FIRE [35]. These master integrals, which depend on s, t and m t , are well-studied in the literature and the results of Refs. [8,9] can be re-used here.
In approach (B), exp does not perform any expansion but simply maps each diagram to a predefined integral family with massive final-state Higgs bosons and an internal Higgs boson propagator with mass m int H . These integrals are expanded in δ at the level of the integrand by FORM, and the resulting scalar integrals are expanded in m ext H by LiteRed and IBP reduced using FIRE. The number of master integrals is minimized using the FIRE command FindRules, which equates identical integrals which belong to different integral families; this procedure yields a basis of 167 master integrals. We also apply FindRules to the entire list of unreduced integrals, as discussed in Ref. [9]. Applying the IBP reduction tables to the equalities found here yields an additional 27 non-trivial relations between master integrals, thus we finally obtain a basis of 140 two-loop master integrals. We additionally perform the IBP reduction of a set of test integrals using Kira [36,37]; here we also find a basis of 140 two-loop master integrals after minimizing between the different families.
These master integrals are four-point integrals with massless external legs, and all propagators have the mass m t . Up to permutations of the external momenta, they belong to one of two integral families, shown in Fig. 3. The computation of these master integrals in the limit s, t m 2 t is described in Section 5. The amplitudes for the two form factors are linear combinations of the master integrals, and we expand their coefficients to order (m ext H ) 4 and (δ ) 3 . This expansion depth requires the IBP reduction of around 350,000 scalar integrals. We also pre-expand the coefficients in m t and , and the final expansions of the form factors are obtained after inserting the m t -and -expanded master integrals.
The freedom in the choice of basis for the master integrals can lead to some undesirable properties; the first is that the denominators of the coefficients of the master integrals in the reduction rules do not factorize in the dimensional regulator and the kinematic invariants and masses, s, t and m t . The second is that the coefficients contain poles in , which imply that the master integrals need to be computed to higher orders in , to produce the finite contribution of the amplitude. The first point complicates the reduction and subsequent expansions of the amplitude, leading to poor computational performance. The second leads to unnecessarily difficult master integral computations involving functions and constants of higher transcendental weight which, ultimately, will cancel in the physical amplitude. We use an improved version of the program ImproveMasters.m [38] to find, for each family, a basis of master integrals for which both of these issues are avoided.

Fully massive two-loop box master integrals
The main purpose of this section is to provide details on the computation of the master integrals, which is based on differential equations [39][40][41][42]. The technically most challenging part is the computation of the boundary conditions which is described in Subsection 5.2.

Differential equations
The master integrals have a non-trivial dependence on two scaleless parameterst = t/s andm 2 t = m 2 t /s. We use LiteRed in combination with FIRE to derive linear systems of coupled differential equations with respect to each of the variables. In principle one can try to solve these sets of differential equations analytically, however the results are not expressible in terms of iterated integrals but rather involve more complicated structures like elliptic integrals. To obtain, nevertheless, precise and easy-to-evaluate results we follow the ideas of Ref. [8,9] and evaluate the master integrals analytically in the highenergy expansion, i.e., for m 2 t s, t.
To construct the asymptotic expansion we insert a power-log ansatz for each master integral into the differential equation in the variablem 2 t and re-expand in andm t . Since there are no spurious poles in either in the physical amplitudes or the differential equations we can choose i max = 0 for all master integrals. We have produced the expansion up to j max = 120 for each of the master integrals, however the amplitudes contain spurious negative powers in m t in the coefficients of the master integrals and additionally, one factor of m 2 t is moved to the prefactor α t . Thus the final expansion depth of the form factors for approaches (A) and (B) arem 112 t andm 114 t , respectively. After inserting the ansatz given in Eq. (12) for each master integral into the differential equation we can compare the coefficients of ,m t and log(m 2 t ) between the right-and left-hand side of the differential equation to obtain a system of linear equations for the expansion coefficients c (n) ijk . We use Kira together with FireFly [43,44] to solve this system of equations in terms of a minimal set of boundary conditions, making sure to favour coefficients which belong to simpler master integrals and low expansion depth in the reduction. The main challenge is then to compute the remaining undetermined set of boundary conditions, which still depend on the second kinematic variablet. We note that the set of boundary conditions required is independent of the value of j max , thus the final expansion depth of the master integrals is limited only by the ability of Kira and FireFly to solve the large system of equations generated by high values of j max . Deeper expansions than we have presented here are certainly possible, if required.
For the calculation of thet dependence of the boundary conditions in the limitm t → 0 we use the methods developed in Refs. [8,9,45]. In particular, we use the method of expansion-by-regions [46,47] to obtain integral representations for the required boundary coefficients. These are subsequently solved with the help of Mellin-Barnes integrals, either by analytically solving summations over residues or by high-precision numerical evaluations together with the PSLQ algorithm [48]. In the following section we will describe in detail how to obtain the integral representations of the asymptotic expansion in general and show details of the calculation of a few examples explicitly.
We have performed numerical cross checks for all 140 master integrals with the help of FIESTA [49]. Using Euclidean kinematics, where both s and t are negative, we typically obtain six digits of agreement for the real-valued master integrals belonging to the integral families in Fig. 3. For these checks we use m t = 173 GeV and set s = t/2 with √ s ≥ 1200 GeV. In the physical region, where s > 0 and t ≤ −s/2 < 0, we find agreement for all 140 master integrals within the numerical uncertainty of FIESTA which provides between two and six significant digits. The lowest precision is obtained for the sevenline master integrals with dots, which are numerically very challenging. We have also performed consistency checks by inserting our analytic high-energy expansions into the system of t-differential equations and found that they are satisfied, order-by-order in m t .

Boundary conditions: Mellin-Barnes approach
In this subsection, we demonstrate the Mellin-Barnes (MB) approach for the calculation of the boundary conditions for them t -differential equations for the master integrals. We only consider the subset of master integrals for which the Euclidean region is defined by S, T > 0 and U < 0, where S = −s, T = −t and U = −u. The remaining master integrals can then be found by crossing relations. The analytic continuation to the physical region is done at the end of the calculation.

Basics of Mellin-Barnes representations and template integrals
We start with a short review of the basics of MB representations and the usage of so-called "template integrals" in the asymptotic m t expansion. 1 For a two-loop master integral with n lines we employ the following α representation, where U and F are Symanzik polynomials, δ i are additional regulators associated with the denominators D i , and the integration measure is chosen as For later convenience, we further adopt the notation for the α-parameter measure as We realize the asymptotic expansion in the high-energy region with the help of version 2.1 of the program asy [50]. Using as input asy provides the possible scalings of the α parameters in the asymptotic expansion; it provides a list of replacements α i → χ n i α i which describe the different regions contributing to the asymptotic expansion.
In the "hard region" we have m 2 t ∼ χ, while all α parameters scale as 1. Therefore, it corresponds to a simple Taylor expansion in m t which can be realized via The integrals on the r.h.s. can be reduced to known massless master integrals (see, e.g., Ref. [51,52]) using IBP methods.
For the "soft regions", i.e. the regions in which at least one of the α parameters scales ∼ χ, we can expand the α representation of Eq. (13) according to the region's α-parameter scaling as 2 U (r) and F (r) are the Symanzik polynomials where χ has been introduced by applying the scaling of region r. Note that contrary to the hard region, which always starts at O(m 0 t ), the soft regions can have different leading powers.
Taking the derivatives w.r.t. χ in Eq. (18) essentially produces the content of the curly brackets multiplied by polynomials in α i , dimensionful quantities, the dimension d, and negative powers of U (r) . This allows us to define "shift operators"Ŝ k r which reproduce the k th derivative in the region r without computing the derivative explicitly. Schematically, these shift operators can be written aŝ where σ runs over the various combinations of at most k th order monomials constructed from v j ∈ {m 2 t , d, S, T, U }, ρ σ ≥ 0 is an integer, and we have introduced the notation: The χ-expansion of a region r can now be interpreted, not in terms of derivatives, but as the shifting of the indices of the single template integral of the region, T r . This template integral represents the leading integral in the region's χ-expansion and is given by We provide Mathematica expressions for all template integrals in the ancillary files [53]. The action of one possible term of the shift operators on the template integrals is given by: where β i ≥ 0 and ρ ≥ 0 are integers and P β i In this way, the higher-order χ-expansion terms for the master integrals without numerators 3 can be obtained from a single template integral per region. The full expansion of a master integral in the soft regions can therefore be written as The MB representation of the template integrals can be obtained by means of direct integration over the α parameters and the application of Mellin-Barnes representations, where the integration path has to be chosen in such a way as to separate the poles of the Γ(· · · + z) and Γ(· · · − z) factors. Note that the individual template integrals contain spurious poles in the regulators δ i , which cancel in the sum of all soft regions.

Mellin-Barnes representations for master integrals with numerators
In the following we introduce a parametric method to directly obtain the MB representations for the boundary conditions of master integrals with numerators. Another method would be to reduce master integrals with numerators to a basis of master integrals with only dots via IBP reductions. However in such a basis, deeper expansions in and m t are often required due to the presence of spurious poles in the IBP relations. 4 The numerators in the α representation can be introduced on the same footing as propagator denominators [31,55] via The α representation of the n-line master integral with m additional numerators can then be obtained as where in our case we have m = 1, 2. In the second line,Ũ andF are Symanzik polynomials in terms of (n + m) α parameters, while in the last line U and F are those of Eq. (13) in terms of only n α parameters. The functionÔ m comes from the derivatives in the second line; it has a similar form as the shift operators of Eq. (19). Note that in Eq. (26) no expansion in χ has been performed.
At this stage, having derived the n-dimensional α representation, we are ready to apply all the techniques developed for the n-line master integrals to Eq. (26). By performing the asymptotic expansions as described in Eq. (16), the resulting hard-region integral I (hard) n,m can be solved in the same way as Eq. (17), and the integrals in the soft regions can be expressed as where the action of the expanded shift operatorŜ k+m r follows the same rule as in Eq. (22), and the template integrals T r are the same n-line integrals defined in Eq. (21). We emphasize that the shifts from operatorsŜ m r yield the leading-order terms in the asymptotic m t expansions of these master integrals with numerators.
Eq. (27) provides an algorithmic way to obtain the MB representations of master integrals with arbitrary numerators. Compared with using IBP reduction to change to a basis of master integrals without numerators, our method has the advantage of avoiding spurious higher-order poles in and m t . Hence, one obtains a much more compact expression in terms of MB integrals, and the cancellation of δ i -poles among different regions can be obtained more easily. 5

Solving Mellin-Barnes integrals
In order to solve the MB representations derived in Eqs. (23) and (27), the first step is to fix the integration contour and perform analytic continuation and series expansions in the δ i and regulators accordingly. 6 This step can be performed with the help of the package MB.m [56]. We now obtain a large number of multi-dimensional MB representations for complicated integrals, which requires a systematic approach for their calculation.
In general our method aims to find infinite sums of residues of the MB integrals, that are suitable for summation procedures. Such residue sum representations are passed to EvaluateMultiSums.m [57] and HarmonicSums.m [58] which internally use Sigma.m [59] for the analytic summation. This step is non-trivial, especially for multi-dimensional MB integrals, and it involves various supplementary techniques such as adding auxiliary scales, the T -expansion of MB integrals and ansatz fitting procedures, as well as numerical evaluation and the PSLQ algorithm. We will describe these methods by providing three examples in the following subsections.
In the following, we adopt the abbreviations and denote the Harmonic PolyLogarithms (HPLs) as H(m 1 , . . . , m n , x) (see Ref. [60] for their definition). We use log(x) and H(0, x) interchangeably.

Example 1: three-line integral
We start by considering the three-line massive sunrise integral with massless external lines. The diagram is shown in Fig. 4, where solid and dotted lines denote massive and massless propagators, respectively. The Symanzik polynomials are given by 5 Note that the complexity of Eq. (27) at O(χ k ) is similar to the O(χ m+k ) expansions in Eq. (23). 6 For details on the analytic continuation of multiple regulators, we refer to [45]. and involve only one scale, m 2 t ; hence there is no need to perform asymptotic expansion. We obtain a one-dimensional MB integral representation: Since no expansion in m t has to be performed the regulators δ i are not required and have been dropped. We first fix the integration contour and the value of such that the left-and right-poles of the Gamma functions are separated by a straight line. In this case Re(z 1 ) = −1/7 and = 1 satisfy this condition. We then perform the analytic continuation → 0 such that we can expand the integrand in . These manipulations can be performed using MB.m [56], which yields where the remaining MB integral In order to solve the integral I (MB) we can close the integration contour to the right and then sum the residues. We obtain: where S i (n) denote harmonic sums, i.e., S i (n) = n k=1 sign(i) k /k |i| . As can be seen from the sum representation, the solution will be given by inverse binomial sums at infinity which have, for example, been studied in Ref. [61,62]. However, their associated constants are not as well studied as those associated to the usual harmonic sums and thus a simplification of the final result is difficult; for this reason we proceed with a different method. The first step is to introduce a parameter into the sum and define This allows us to find a generating function of the sum in Eq. (34) with the help of the command ComputeGeneratingFunction implemented in HarmonicSums.m. The result is expressed in terms of iterated integrals over the letters Afterwards we rationalize the square-root valued letters with the command SpecialGLToH and take the limit ξ → 1 to reconstruct I (MB) in Eq. (31). The result is given by We see that the solution can be written in terms of iterated integrals with cyclotomic letters [63]. They can be further reduced to known constants that are represented by multiple polylogarithms evaluated at the sixth roots of unity [64] which yields, for the -finite part of the massive sunrise diagram, the result Here ψ (1) 1 3 is the PolyGamma function that is related to the Clausen function by When reconstructing analytic expressions from numerical evaluations using the PSLQ algorithm we therefore have to to use the following basis of constants as well as all possible products up to transcendental weight 4:
The scaling (0, 0, 0, 0, 0) corresponds to the hard region, in which only m 2 t ∼ χ and all α parameters scale as α i ∼ χ 0 . In the remaining eight regions a subset of the α parameters scale as α j ∼ χ.
Hard region: For the hard region, we proceed in the same way as Eq. (17). The leading term at O(m 0 t ) can be obtained by setting m t = 0, which corresponds to one of the known massless master integrals given in Refs. [65,66]. For the sub-leading term at O(m 2 t ), we first perform a Taylor expansion at the integrand level, and then perform an IBP reduction with LiteRed [30,31] to reduce again to the set of known massless master integrals to obtain the final result.
Soft regions: For the soft regions, we apply the eight scalings from Eq. (40) to the Symanzik polynomials in Eq. (39), and expand the α representation to the sub-leading order in χ as described in Eq. (18). For region (1) we find, for example, with the expanded Symanzik polynomials Note that U 1 is the coefficient of χ 0 and F 1 is the coefficient of χ 1 . The eight template integrals, which correspond to the leading contributions, can be extracted according to Eq. (21). They are represented by, at most, one-dimensional MB integrals.
The template integral for region 1 is given by which is obtained from Eq. (21) through straightforward integration. The expansion in Eq. (41) can also be reinterpreted in terms of shift operators acting on this template integral: Results: Solving the MB integrals in the soft regions and combining them with the hard region, we obtain the solution of the five-line integral of Fig. 5: which is free from δ i -and -poles. As a final example we consider the seven-line double box integral (see Fig. 6) with two additional numerators, which needs to be evaluated to O( 0 , m 0 t ) for the boundary conditions. 7
Asymptotic expansions: With the representation in terms of seven α parameters for this "7+2"-line integral in hand, we can again apply the asymptotic expansions for the scaling of Eq.
For the hard region, we proceed in the standard way, i.e. we take the massless limit and perform IBP reductions to the known massless master integrals. For the 13 soft regions, we expand the α representation in Eq. (48) according to Eq. (27), The expanded shift operatorŜ 2 r is the leading term of the operatorÔ 2 in Eq. (49) where the r th region scales according to Eq. (50). The 13 template integrals can be identified by U r and F r according to Eq. (21). By performing parametric integrations and Mellin transformations, we obtain up to three-dimensional MB representations for the template integrals. By applying the shift operators in Eq. (22) to Eq. (51), we obtain the MB representations of the soft regions. 8 The next step is to perform an analytic continuation w.r.t. the eight regulators δ 1 , . . . , δ 7 and . We fix the integration contours at {Re(z 1 ) = −1/7, Re(z 2 ) = −1/11, Re(z 3 ) = −1/17} as straight lines. Then we perform the continuation with the MB.m package and expand the expression to order O(δ 0 i ) and O( 0 ). This yields a large number of one-, two-and three-dimensional MB integrals; 2003, 515 and 14 respectively. In the following paragraphs, we will demonstrate our method to solve multi-dimensional MB integrals, focussing in particular on non-trivial examples which have a non-zero contribution from the contour-closing arc at infinity which must be taken into account.
Arc and residue sums: Here we start with a simple but non-trivial example which appears in our calculations, which demonstrates the importance of the arc contribution. The example is a one-dimensional scaleless MB integral with the integrand 8 An explicit example of applying the shift operators is shown in Eq. (44).
where the integration contour is fixed at Re(z 2 ) = −1/11. Cauchy's residue theorem states that where the (−) sign comes from the fact that we close the contour clockwise. One usually assumes that the arc contribution vanishes. However, this is not the case for Eq. (52). Closing the integration contour to the right and summing the residues we obtain On the other hand, regularizing the integrand by multiplying with ξ z 2 and summing the residues we obtain The same result can be found by precise numerical integration and employing the PSLQ algorithm. The difference between the two results in Eqs. (54) and (55) is the missing contribution from the arc in Eq. (53): Therefore, in order to systematically take the arc contribution into account, we always rely on numerical integration of the MB integrals accompanied by the PSLQ algorithm to cross-check results obtained from the residue summations for scaleless MB integrals. However, the problem becomes more complicated when a non-vanishing arc contribution like Eq. (56) is nested in two-dimensional MB integrals involving the kinematic invariants T /S. In the following we will introduce a method which can deal with such situations.
Nested arc contribution: For two-dimensional MB integrals, we always first try to reduce their dimensionality using Barnes' lemmas as implemented in barnesroutines.m [67] and other simplification tricks. For the remaining two-dimensional MB integrals involving kinematic invariants and a nested arc contribution, we need a more careful analysis. Let us now consider two-dimensional MB integrals of the form whereΓ denotes the product of Gamma functions with common integration variables. In our case we have two types of z 1 residues from the Gamma functions, which are given by    type 1: z 1 = 0, 1, 2, . . .
From the type 1 residues with integer z 1 we obtain whereF (k 1 , z 2 ) denotes the resulting residue function. From the type 2 residues in Eq. (58), we have We can then take the nested z 2 residues in Eqs. (59) and (60), which introduces a second infinite sum over k 2 , and then perform the residue summations over both k 1 and k 2 with the help of Sigma.m and EvaluateMultiSums.m. However, this two-dimensional (k 1 , k 2 ) residue summation will miss the arc contributions in the first type, given in Eq. (59), from scaleless one-dimensional MB integrals in z 2 . The residue summation for the second type, given in Eq. (60), is correct, since the kinematic scale choice 0 < T /S < 1 will suppress the asymptotic behaviour of the integrands and ensure that the arc contributions in Eq. (60) are vanishing. Instead of introducing another regulator into the two-dimensional MB integrals, which would increase the computational complexity significantly, we use precise numerical integration together with the PSLQ algorithm in order to find the correct results at fixed values of k 2 . Clearly we can not compute the infinite sum in this way, so we introduce the method of T -expansion and ansatz fitting procedures to obtain the correct result for Eq. (57).
Ansatz fitting and T -expansions: The basic idea of this method is to start with an ansatz for the sum of MB integrals of the type given in Eq. (57) which contains rational functions and HPLs up to weight 4, and perform a series expansion in T to a finite power n. Then we expand Eqs. (59) and (60) up to O(T n ) by taking residues, and compute the remaining one-dimensional MB integrals. The result can then be fitted to the series expansion of the ansatz; the fitting procedure consists of solving a system of linear equations to determine the unknown coefficients of the ansatz.
An ansatz which includes weight 4 functions is rather large, requiring a series expansion to a high power n to completely fix its coefficients. In practice, our experience shows that the arc does not contribute to the higher-transcendental-weight contributions, allowing us to limit the size of the ansatz and thus the required depth of the series expansions.
In the following, we demonstrate this idea with an explicit example that is present in our calculation. We have a two-dimensional MB expression I and perform the residue summation as described above. This leads to where and x = T /S. I with the nine free parameters c 1 , . . . , c 9 .
Using numerical integration and the PSLQ algorithm we can construct a series expansion of I which is given by Note that here the arc contributions are included correctly. By performing a series expansion of Eq. (63) and comparing to Eq. (64) we obtain an over Results: After solving all MB integrals and adding the result from the hard region, we derive the final solution of this "7+2"-line master integral where the constants C T and C S originate from three-dimensional MB integrals which are discussed in Appendix A.

Crossing and analytic continuation
As stated above, we only calculate the boundary conditions for the subset of master integrals for which the Euclidean region is defined for S, T > 0 and U < 0. The boundary conditions for all other master integrals can be obtained by applying one of the five crossing relations: While the rational dependence can be easily obtained via these replacements, the HPLs need analytic continuation.
Due to our choice of the Euclidean region we start with HPLs of the argument x = T /S, which are real in this region. To analytically continue to the physical region, we have to arrive at the argument x = −T /S = T /s = −x. The transformation of HPLs to the negative argument is implemented in HarmonicSums and HPL. However, we have to take care to use the correct sign for the analytic continuation. We have s = s+i ε, so x = x+i ε and therefore have to use the '+' sign for the analytic continuation which leads to Using HarmonicSums or HPL we can transform the argument of all occurring HPLs to the physical region. For example, we have The analytic continuation of the HPLs after the application of the different crossings can be obtained in a similar manner, but require more involved transformations. For example, after the crossing T → U we end up with HPLs of the argument y = −(1 + T /S + i ε). We can map these HPLs back to argument x by first applying the transformation y → −y = y and afterwards y → 1 − y = x . The sign for the analytic continuation has to be chosen as '−' for the first and '+' for the second transformation. This results, for example, in As a final example, let us look at the crossing S ↔ T . Here, we find HPLs of argument w = S/T − i ε. We can map these HPLs back to argument x by first applying the transformation w → 1/w = x and afterwards continue as for the first example. We find The analytic continuation for the other crossings can be derived analogously. In total we can express all 140 master integrals through the following set of HPLs: While the expression in terms of HPLs is more convenient for analytic manipulations, the expressions in terms of polylogarithms (Li n (x)) and Nielsen polylogarithms (S n,m (x)) might be more convenient for numerical evaluations, since many standard math libraries already contain implementations.
In the supplementary material to this paper [53] we provide the analytic results for all 140 master integrals.

Form factors for gg → HH
The contribution to the form factors of gg → HH from diagrams of Fig. 1 is infrared finite and has only ultraviolet divergences. They are removed by renormalizing the top quark mass and Yukawa coupling in the leading order contributions. The counterterms are well known in the on-shell scheme, see, e.g., Ref. [68]. In this work it is sufficient to perform the renormalization in the MS scheme. The corresponding mass counterterm is given by (see, e.g., Eq. (31) of Ref. [69]) where α is the fine structure constant, s W ≡ sin θ W is the sine of the weak mixing angle and N C = 3. The second term inside the round brackets originates from the tadpole contribution 9 and is only provided for completeness; it is not used in this paper. Note that one factor m 2 t is collected in α t (see Eq. (9)) such that the expansion up to (m 2 t ) 56 and (m 2 t ) 57 are available for the Padé method. We follow Ref. [71] and construct the so-called "pole distance re-weighted" Padé approximants and the corresponding uncertainties (see Section 4 of [71] for a detailed discussion), in which Padé approximants [n/m] are included which satisfy N low ≤ n + m ≤ N high and N low ≤ n + m − |n − m| . massive. This is one of our master integrals, which we have expanded up to (m 2 t ) 60 . For this integral it is possible to obtain precise numerical results using FIESTA. In Fig. 7(b) we compare, for p T = 120 GeV, the real and imaginary parts of the Padé method to the numerical results. For the Padé method we use {N low , N high } = {50, 60}, the same choice as we make for the form factors. For values of √ s ≈ 500 GeV and higher the Padé uncertainties are very small and we find perfect agreement between the Padé and FIESTA results. For lower √ s the Padé uncertainties in the real part grow. It is nevertheless interesting to see that the central values are close to the numerical results. On the contrary, for the imaginary part the Padé uncertainties remain small but there is a clear deviation from the exact result. This can be explained as follows: The integral we consider admits two-and three-particle cuts. For the latter √ s = 3m t = 519 GeV which is about the starting point for the deviations; the Padé method is not expected to be able to approximate the exact function below the cut, which we clearly see in the imaginary part in Fig. 7(b).
In Fig. 7(a) we show the analogous result for the seven-line master integral of approach (A) where middle line is massless. This integral only has cuts through two massive lines (and possibly also a massless line) and indeed, we observe good agreement of the Padé and FIESTA results, even close to the top quark pair threshold at 2m t = 346 GeV.
Let us now move to the form factors F box1 and F box2 and discuss the quality of the expansions in m H and δ. For this purpose we fix p T and plot various different depths. We normalize all curves to the highest-available depth of approach (B), which includes m 4 H and δ 3 .
In Fig. 8(a) the result is shown for the real part of F box1 for p T = 500 GeV. The colours correspond to approach (B) and the results from approach (A) are shown in gray and black. The y axis spans a range below 1% and all approximations which include at least m 2 H terms in approach (A) and m 2 H and δ 1 terms in approach (B) are visible in the plot and thus show a deviation well below the percent level.
In Fig. 8(b) we show results for p T = 200 GeV and √ s values between 480 GeV and 580 GeV. For larger values of √ s the form factor crosses zero and the ratios inflate. Beyond the zero crossing the ratios are a similar size to those of Fig. 8(a). The result from approach (A) show a deviation of about 10% in case the Higgs mass is neglected. It reduced to below 5% after including the m 2 H terms and is of order 1% after including also the quartic terms. The situation is similar for approach (B): Once quadratic terms in m H and δ are include the deviation from 1 is below 5%. Including more expansion terms in m H and δ further stabilizes the approximations.
We conclude that the inclusion of the quartic terms in m H and cubic terms in δ provides an approximation to the (unknown) exact result below the percent level (see also Fig. 2 or Ref. [72] which shows a comparison for gg → ZH).
Next we discuss the results for F box1 and F box2 for a range of values for the transverse momentum p T . In Fig. 9 we show the real and imaginary parts of F box1 and F box2 for p T between 120 GeV and 800 GeV. The colours correspond to the results from approach (B); here we also show the uncertainty band from the Padé method. The results from approach (A) are shown as faint uncertainty bands. They are only visible for small values of p T , where one observes deviations between the two approaches.
Above p T = 200 GeV the uncertainty from the Padé method is negligible. For p T = 150 GeV differences between the approaches are only visible for the real part of F box2 . The situation is similar for p T = 120 GeV for √ s ∼ > 400 GeV where the uncertainty bands are still small. Up to this value the results for F box1 and the imaginary parts of F box2 agree quite well. The real part of F box2 shows larger uncertainties for large values of √ s in approach (B); for approach (A), however, the uncertainties remain small. Note, that F box2 is numerically less important than F box1 . Fig. 9 shows that both ways to treat the internal boson mass leads (within uncertainties) to equivalent physical results. In view of the discussion above we expect that approach (B) only approximates the unknown exact result down to √ s ≈ 520 GeV. However, approaches (A) and (B) agree for even smaller values of √ s. It seems that the master integrals of approach (B) with non-analytic behaviour at the three-particle threshold are numerically suppressed.
In Fig. 10 we show the real and imaginary parts of F box1 and F box2 for fixed scattering angle θ = π/2 for √ s between the top quark threshold and 1200 GeV. The solid curves represent Padé results and the dashed curves the expansions. We observe that the expansions start to diverge 10 for the real parts for box1 . From Ref. [9] (see also Section 3.3 of Ref. [73]) we find that F (1,0) box1 is about O(1) if the scattering angle is fixed to θ = π/2 and for √ s a few hundred GeV. This is also the case for F (0,yt) box1 as can be seen from Fig. 10. For the pre-factors in Eq. (8) we have α t /α s ≈ 0.6 and thus it might very well be that the electroweak corrections provide sizeable contributions to the Higgs pair cross section. Of course, we should emphasize that in this paper only a certain diagram class has been considered; in particular, no triangle diagrams are included. Furthermore for this estimate no interference contributions are taken into account.

Conclusions
In this paper we take the first step towards the electroweak corrections to Higgs boson pair production. We consider the subset of diagrams where a Higgs boson is exchanged between the top quarks. Effects from Higgs self couplings are neglected.
We are interested in analytic calculations of the form factors in the high-energy limit; we perform expansions in m 2 t /s, m 2 t /t and m 2 t /u taking into account up to about 60 expansion terms. We study two methods for the treatment of the internal massive Higgs boson, which is a new feature as compared to the QCD corrections. In our first approach we assume that it is small as compared to the top quark mass, whereas in the second approach it is assumed that the internal Higgs boson is of the same order of magnitude as the top quark mass. In both cases we perform expansions in the respective small parameters. For physical values of the mass parameters both expansion methods agree at the percent level for smaller values of p T and at the permille level for larger values.
The approach with a small internal Higgs boson leads to master integrals which have been computed in the context of QCD corrections. The other approach leads to 140 new master integrals. We describe in detail our approach to compute them analytically using differential equations and the Mellin Barnes method.
We supplement the expansion for small m t by combinations of Padé approximations and the associated uncertainty estimates, which significantly increases the region of phase space where the analytic expansions can be used. We show that Padé approximants based on up to about 60 m 2 t expansion terms provide excellent result down to p T = 150 GeV and even for p T = 120 GeV results with moderate uncertainties are obtained. On the basis of a scalar (master) integral we validate that the uncertainty estimate covers the exact result.
The analytic result for C T is obtained from a consistency condition obtained from the system of m t -expanded t-differential equations for the 140 master integrals. On the other hand, for C S we first perform various shifts of integration contours and analytic continuations to bring the three-dimensional MB integrals into a better form, which can be reduced to, at most, two-dimensional integrals in terms of only Gamma functions by the Barnes lemmas. The resulting MB integrals are the solved by the analytical summations and PSLQ algorithm. Note that it is straightforward to directly compute C S and C T numerically and obtain a precision of about ten digits, which is sufficient for practical applications.