Towards a new approximation for pair-production and associated-production of the Higgs boson

We propose that loop integrals with internal heavy particles can be evaluated by expanding in the limit of small external masses. This provides a systematically improvable approximation to the integrals in the entire phase space, and works particularly well for the high energy tails of kinematic distributions (where the usual $1/M$ expansions cease to be valid). We demonstrate our method using Higgs boson pair production as an example. We find that at both one-loop and two-loop, our method provides good approximations to the integrals appearing in the scattering amplitudes. Comparing to existing expansion methods, our method are not restricted to a special phase space region. Combining our efficient method to compute the two-loop amplitude with an infrared subtraction method for the real emission corrections, we expect to have a fast and reliable tool to calculate the differential cross sections for Higgs boson pair production. This will be useful for phenomenological studies and for the extraction of the Higgs self-coupling from future experimental data. Our method can also be applied to other processes, such as the associated production of the Higgs boson with a jet or a $Z$ boson.


Introduction
Higgs boson pair production and Higgs boson production associated with a jet (H + j) are both important processes at the Large Hadron Collider (LHC). Higgs boson pair production can be used to measure the trilinear self-coupling of the Higgs boson, which is essential to understand the electroweak symmetry breaking in the Standard Model (SM). It also happens that in the SM, contributions from different diagrams to the Higgs pair production cancel each other delicately, leading to a rather small production cross section [1,2]. In a new physics model beyond the SM, such cancellations are not generically expected, and the new physics contribution can potentially be much larger than the SM one. This will lead to large deviations in the total rate as well as kinematic distributions, which have been studied in many new physics models (see, e.g., [3] and references therein) and in the effective field theory framework [4,5]. It is expected that Higgs boson pairs will be detected during the high luminosity phase of the LHC, or even earlier if new physics effects enhance the production rate. This will give us hints about the nature of electroweak phase transition in and beyond the SM.
The H + j production process is also useful both for testing the SM and for probing new physics. The recoil of the additional jet gives rise to a non-zero transverse momentum (p T ) of the Higgs boson, which is a highly important observable actively being measured by the ATLAS [6,7] and the CMS [8,9] experiments at the LHC. The high energy tail of the p T distribution will allow us to probe possible loop corrections from heavy particles which may not be visible in the single Higgs boson production rate.
Due to the importance of these processes, lots of theoretical efforts have been devoted to their precision predictions. At the LHC, the main production mechanism for Higgs boson pairs and for H + j is gluon fusion, in which the Higgs bosons couple to gluons via a top-quark loop. Typical Feynman diagrams at the leading order (LO) are depicted in Figure 1. The LO contributions for both processes have already been known for 30 years, [1,2] for Higgs boson pair and [10,11] for H + j. However, the calculations beyond the LO in quantum chromodynamics (QCD) turned out to be highly complicated, due to the multiple energy scales involved in these processes. This makes the higher order amplitudes rather lengthy, and also forbids an analytic representation of the multi-loop integrals in terms of commonly known mathematical functions. In [12,13], the exact NLO QCD corrections to Higgs boson pair production were finally calculated, using a fully numerical method to evaluate the difficult two-loop integrals. The exact NLO QCD corrections to H + j production were also calculated with a similar method in [14]. The numerical nature of these calculations make them rather time-consuming and not very flexible to perform comprehensive phenomenological studies.
Because of the difficulties in the evaluation of multi-loop integrals with many scales, various approximation methods in certain kinematic limits have been proposed to study production processes of the Higgs boson. A common goal of these methods is to reduce the number of scales involved in the loop integrands, such that the resulting integrals are simpler to calculate. The earliest and easiest one is the Higgs effective field theory (HEFT) approach. In this approach, one takes the top quark m t to infinity, such that the top quark is integrated out leading to effective operators involving Higgs fields and gluon fields. In this way, the NLO QCD corrections only involve one-loop integrals (compared to two-loop ones in the exact calculation) and the scale m t does not appear in propagators anymore. These greatly simplify the higher order calculations. In fact, the single Higgs boson production cross section has been calculated in this approach to the next-to-next-to-next-to-leading order (NNNLO) in QCD [15]. In [16][17][18][19][20][21], the H + j process were calculated at the next-to-next-to-leading order (NNLO) in the HEFT way. As for the Higgs boson pair production, the NNLO corrections in the infinite m t limit were calculated in [22][23][24].
The HEFT approach provides reasonable approximations for the total cross section, as well as for differential cross sections in the kinematic regions where the Higgs bosons are not highly boosted. However, this approach has some drawbacks. First, with the increasing of the perturbative order, the 1/m t power corrections to the effective theory becomes non-negligible. This is the case for the inclusive single Higgs boson production, where the finite m t effect at the NLO [25] is similar in size to the NNNLO contribution [15], and has to be taken into account for phenomenology. More severely, if one considers the kinematic distributions, the HEFT is simply not valid when the energy of the Higgs boson is comparable to or larger than the top quark mass. These high energy tails of the differential distributions, on the other hand, are sensitive to new physics effects and are phenomenologically much more interesting. Therefore, the HEFT approach requires some refinements in these two aspects.
One improvement to the HEFT approximation is performing the 1/m t expansion of the loop integrals to higher powers. This corresponds to including higher dimensional operators in the effective field theory. In this approach, one takes the limit where m t is much larger than the energies of external particles, and performs a power expansion in terms of p µ /m t for the loop integrand, where p µ is some external momentum. The remaining integrals only involves one mass scale, and are very easy to evaluate. This has been done at the NLO for H + j production [26] and for Higgs boson pair production [27][28][29]. It was found that the expansion converges rather well for the total cross sections as well as in the low energy regions of differential distributions. For these observables, the 1/m t expansion therefore provides a fast method to obtain predictions with enough precision.
When the energies of external particles increase above the top quark mass, however, the 1/m t expansion quickly fails to converge. And one does not expect that this will work for the high energy tails of distributions. If the energy is high enough, on the other hand, one could exploit the opposite limit |s|, |t|, |u| ≫ m 2 t to simplify the calculation. Here s, t, u are the usual Mandelstam variables. In this limit, one can perform a double expansion in terms of m 2 t /s and m 2 h /m 2 t , where m h is the mass of the Higgs boson. This has been done for the two-loop amplitudes in H + j production in [30,31], and for the two-loop planar integrals in Higgs boson pair production in [32]. It should be noted that the m t → 0 limit is not regular, and the expansion in m 2 t /s is not a Taylor series. Instead, one generally encounters powers of ln(m 2 t /s) in the expansion. This method is efficient for the high energy tails of differential distributions, e.g., in the region of very large Higgs transverse momentum or the region of very large Higgs pair invariant mass.
Recently, a new method was proposed in [33] to calculate the NLO QCD corrections to Higgs boson pair production, in which the integrals are expanded in the limit of small transverse momentum (compared to √ s and m t ). The resulting integrals are functions of s and m 2 t , and can be calculated analytically. This method provides a rather good approximation for the bulk of the cross section √ s 750 GeV. However, it cannot be applied to the high energy tails of differential distributions, as the expansion becomes divergent in the region √ s 900 GeV or in the region of large transverse momentum.
Given the above expansion methods valid in several distinct kinematic regions, it is interesting to ask whether it is possible to construct an expansion which is valid in the whole phase space. The answer is yes, as will be outlined in this article, with the price that the resulting integrals are more complicated (but are still simpler than the original integrals, and are possible to be calculated with the help of differential equations). Our starting point is an expansion in the limit of small m h , without assumptions on the other scales s, t, u and m t . The resulting integrals involve only massless external legs and resemble those appearing in, e.g., the di-jet production process [34]. For these integrals, we first perform the standard integration by parts (IBP) reduction into a set of master integrals [35], and then construct differential equations for the master integrals [36,37]. Wherever possible, we convert the differential equations into a canonical form [38] via an appropriate basis choice, and obtain the solutions in terms of the Chen iterated integrals [39]. This is the case for 3 out of 4 integral families in Higgs boson pair production. As for the remaining integral family, we reduce the system of differential equations such that at most two master integrals are coupled at the leading order in the dimensional regulator. This allows a solution in terms of elliptic integrals.
The method developed along this line can be viewed as a unification of the existing expansion-based approaches, and is valid in a broader region of phase space. As such, it represents an improvement over the other methods, especially in the phenomenologically important intermediate regions. Our method is not restricted to Higgs boson pair production and H + j production. 1 It can be applied to any process involving a heavy quark loop, such as the top quark loop contribution to gg → HZ [41] and gg → ZZ [42], as well as the mixed QCD-electroweak corrections to e + e − → HZ [43,44].
The article is organized as follows. In Section 2, we introduce our method to derive the small Higgs mass expansion, and demonstrate its validity using the known one-loop results. In Section 3, we apply our method to the two-loop non-planar integrals in Higgs boson pair production, and show the comparison of our results in one of the topologies against the numerical results from sector decomposition. We conclude and discuss future developments of our method in Section 4. And finally in the Appendix, we list the basis of master integrals we use in our calculation.

Expansion in terms of external Higgs masses
In this and the following sections, we will use Higgs boson pair production as the concrete example to demonstrate our method, namely, we consider two-loop contributions to the process where the kinematic invariants are They satisfy the usual relation s + t + u = 2m 2 h . As a result, the scattering amplitude depends on 4 energy scales which can be chosen as s, t, m h and m t , where the top quark mass enters through propagators.
The presence of multiple scales in the two-loop amplitude makes it rather difficult to calculate. On one hand, it is highly non-trivial to reduce the amplitude into a set of master integrals via the usual IBP method. On the other hand, many of the master integrals are not expected to have a representation in terms of (generalized) polylogarithms or even elliptic integrals. Given this situation, the various approximation methods mentioned in the introduction exploit different kinematic limits to reduce the number of scales in the problem. This simplifies both the reduction of the amplitude and the evaluation of the master integrals. For example, the large m t expansion corresponds to the limit m 2 t ≫ |s|, |t|, m 2 h ; the large energy expansion corresponds to |s|, |t| ≫ m 2 t ≫ m 2 h ; and the small p T expansion corresponds to |s|, m 2 t ≫ |t|, m 2 h . Note that all the above expansions can be obtained by first expanding around the limit m 2 h → 0, and then further dealing with the 3 remaining scales s, t and m 2 t . This small Higgs mass limit is therefore more general and is valid in a broader region of phase space.
In general, loop integrals may develop new singularities in the limit where one of the internal or external masses is taken to zero. If that's the case, the expansion around that limit will not be a normal power series. An example is the limit of small top quark mass discussed in [32], where the expansion involves powers of log(m t ) in addition to powers of m t . However, external Higgs bosons are special, since they only couple to massive particles directly. As a result, no new singularities arise in the limit m h → 0, and we can Taylor-expand a generic integral as 2 where The above expansion coefficients can be obtained by calculating integrals with only massless external legs, which are simpler than the original ones.
To demonstrate our method, we take the one-loop integrals appearing in Figure 1(a) as an example. The propagators are given by where k is the loop momentum. We define the family of integrals where d = 4−2ǫ is the space-time dimension in dimensional regularization. In order to perform the expansion, we need to take the derivative of the above integrals with respect to m 2 h . This can be accomplished via the following operator As an example, we have where we have suppressed the arguments of the integrals. The expansion of I 1,1,1,1 can then be written as Note that no reduction has been performed at this stage, which is important since the IBP reduction for the original m h -dependent integrals can become very complicated at the two-loop order. The reduction can be carried out for theĨ integrals when necessary, which is much easier to do. For example, after reduction, Eq. (9) can be simplified to In the above formula, we only show the expansion up to order m 2 h . There is no difficulty in extending the expansion to higher powers of m h . In practice we find that keeping terms up  to m 4 h or m 6 h already provides rather good approximations to the exact results. This will be clear from the numerical results shown in the following.
We choose the masses to be m t = 173.3 GeV and m h = 125.1 GeV, and vary the kinematic variables s and t to see the goodness of the approximation in different regions of phase space. We first show the result for the integral I 1,1,1,1 . For convenience we rescale the integral by an appropriate power of m t such that the result is a dimensionless number. In the upper plots of Fig. 2, we fix t = −(200 GeV) 2 and show the real part and the imaginary part of the order ǫ 0 coefficient as a function of √ s. It can be seen that the expansion up to order m 4 h already gives sub-percent accuracies for both the real and the imaginary parts for almost all values of √ s, ranging from the threshold region √ s 2m h , to the tt threshold √ s ∼ 2m t , to the high energy regime √ s ≫ 2m t . The only exception is the real part at around √ s ∼ 560 GeV, where it happens by coincidence that the value of the integral is close to zero. In such cases one needs to add the order m 6 h term, which leads to per-mille accuracy in all regions of phase space. Similar behavior can be observed in the lower plots of Fig. 2, where we fix √ s = 1000 GeV and show the integral as a function of √ −t. Here the approximation at order m 4 h gives betterthan-per-mille accuracy in the whole range and it is not necessary to include the O(m 6 h ) corrections.
The behavior of a single integral is perhaps not convincing enough. We now turn to investigate the partonic (differential) cross sections which are physically more relevant. We start by writing the amplitude as where the two tensor structures are given by [1] A µν where p T denotes the transverse momentum of the top quark and can be written as At one-loop, the two form factors F 1 and F 2 can be evaluated either exactly or using the small m h expansion. They can then be used to calculate the partonic differential cross section and also the partonic total cross sectionσ by integrating over p T . At this point, it is interesting to compare our approximation to the other methods, e.g., the 1/m t expansion in [27] and the p 2 T /s expansion in [33]. On the left side of Fig. 3, we show the partonic total cross sectionσ as a function of √ s. The black solid line is the exact result, while the other curves represent three different approximations. It is clear that the large-m t expansion only works in the region s < 4m 2 t , as expected. The p 2 T /s expansion is valid in a broader range, and provides a reasonable approximation to the exact result up to √ s 900 GeV. However, going beyond that, the p 2 T /s expansion quickly becomes divergent. On the other hand, our small-m h expansion works perfectly across the whole range. To see more clearly the behaviors of the small-p T expansion and our small-m h expansion, in the lower panel of the plot we show the relative error with respect to the exact result. We find that the qualities of the two approximations are similar for √ s < 500 GeV. Beyond that, the small-p T expansion becomes worse and worse, while the small m h expansion becomes better and better, and provides a better-than-per-mille approximation to the exact result.
To see more clearly the difference between the p 2 T /s expansion and the small-m h expansion at high energy, we show on the right side of Fig. 3  is at the level of 10 −5 in the whole range of p T . We also observe that the distribution peaks towards the right end, which means that the dominant contribution to the partonic total cross section comes from the high p T region. It is clear that the small-p T expansion cannot be a good approximation in this region, which is due to the fact that the condition p T ≪ m t is no longer fulfilled. This also explains why the small-p T expansion fails for the partonic total cross section at large √ s, as observed from the left plot.
The above discussions demonstrate the validity of the small-m h expansion in the entire phase space at the one-loop level. This makes us confident that the same will be true at higher loop orders. In the following section, we apply our expansion to the two-loop amplitude, with the goal to provide a fast and reliable method to evaluate the NLO QCD corrections to Higgs boson pair production.
3 Expansion at the two-loop order 3

.1 Setup
We now turn to the NLO (two-loop) QCD corrections to Higgs boson pair production. The expansion in terms of m 2 h takes the form as Eq. (3) and can be performed using the derivative operator Eq. (7). We stress that this can be done at the amplitude level, without the need of reduction beforehand. We have carried out the expansion up to order m 4 h . The extension to higher powers in m h is straightforward. The expansion coefficients can be obtained by calculating integrals with massless external legs. After applying crossing symmetries, all the integrals can be classified into 6 integral families. They corresponds to the 6 topologies depicted in Fig. 4.
The integrals in topology C and D can be related to those in topology A and B via IBP relations. We therefore only need to consider 4 integral families. We first define the m h - where k 1 and k 2 are loop momenta, and {a i } denotes the collection of powers a i on the propagators D i . We then definẽ which are the main objects to be calculated in this section. The 4 relevant integral families are defined by their corresponding propagators as the following: A : The kinematic invariants are defined as in Eq. (2), with the exception that we now have p 2 3 = p 2 4 = 0 and s + t + u = 0. We choose s, t and m 2 t as independent scales and introduce the following dimensionless quantities Physically, we have s > 2m 2 h and t, u < 0. When s is above the tt threshold, namely s > 4m 2 t , some of the integrals will develop imaginary parts. For convenience, we will first work in the unphysical region which corresponds to This guarantees that all the integrals are real. After obtaining the expressions of the master integrals, we can perform an analytic continuation to the physical region and then evaluate them numerically. Among the 4 integral families, the two planar topologies A and B have already been discussed in [34,49]. In the following, we discuss the calculation of the master integrals in the two non-planar topologies E and F.

Analytic structures
Topology E is the simpler one in the two non-planar topologies, in that it is possible to cast the differential equations satisfied by the master integrals into a canonical form [38]. For that purpose we use a method similar to the one used in [50]. We start from the sub-topologies with the lowest number of propagators. We choose appropriate pre-canonical master integrals so that the differential equations are simple enough and are of the form where the vector f 0 ( x, ǫ) denotes the collection of the pre-canonical master integrals, and x is the collection of independent kinematic variables (in our case x = {µ, ν}.Ã 0i ( x) andB 0i ( x) are two square matrices which do not depend on the dimensional regulator ǫ. We then apply a linear transform T ( x, ǫ) on the vector f 0 (x, ǫ) such that the new vector f ( x, ǫ) = T ( x, ǫ) f 0 ( x, ǫ) satisfies a system of differential equations in the canonical form where the matricesÃ i ( x) are at most algebraic functions of the variables x i . We then proceed to topologies with the number of propagators higher by 1, and repeat the above process. Finally Figure 5: Pre-canonical master integrals in topology E. The thick lines represent massive propagators (from top quarks) and the thin lines denote massless propagators (from gluons). The labels s, t and u on the external lines represent the (squared) momenta flowing through those legs. The external lines without labels have light-like momenta.
we arrive at the top-level topology with 7 propagators. In Fig. 5 we list the diagrammatic representations of the pre-canonical integrals in topology E. The transformation of them into the canonical basis is discussed in the Appendix. Given the differential equation in the canonical form (23), the solution for the master integrals can be generically written as Chen iterated integrals [39]. For that purpose, it is convenient to rewrite the differential equations as where dA( x) = iÃ i ( x) dx i which can be expressed in the d-log form Here for a given k, A k is a constant matrix independent of the kinematic variables, while α k ( x) is an algebraic function of the kinematic variables and is called a "letter". The collection of all letters is called the "alphabet", which completely determines the class of functions appearing in the master integrals f ( x, ǫ). The general solution can be written as where P denotes path-ordering along the path connecting the boundary point x 0 and the destination x. The boundary conditions f ( x 0 , ǫ) can possibly be fixed by the analytic structure of the differential equations, or can be calculated directly. The formal solution (26) is exact in ǫ, while in practice one usually needs its expansion around the 4-dimensional limit ǫ = 0. In order to do that, it is convenient to normalize the master integrals to have the property of uniform transcendental weights. The concept of transcendental weight (we will simply call it "weight" in the following) is closely related to iterated integrals. The weight of an algebraic number is defined to be 0, the weight of π is defined to be 1, while the weight of the Riemann zeta value ζ n is n. Given a weight-n function g( x), the weight of the integral is defined to be n + 1, where α( x) is an algebraic function of the kinematic variables. With this definition, it is clear that the n-fold iterated integral of the form has transcendental weight n. Now considering the expansion of the master integrals around ǫ = 0, We will normalize the master integrals such that the components of the vector f (i) ( x) are all weight-i functions (or numbers). This is possible since they satisfy the canonical-form differential equation (23). For topology E, the prefactors for the normalization are collected in the Appendix. After normalization, the boundary conditions are simply given by lim µ,ν→∞ where the boundary µ, ν → ∞ corresponds to s, t → 0. The coefficient functions f (i) ( x) can then be written as iterated integrals order-by-order: Given the above formal solutions, it is still non-trivial to convert them to explicit functions such as logarithms, polylogarithms and multiple polylogarithms (MPLs) [56]. For this purpose, we will use the concept of "symbol" [54][55][56][57], which maps the iterated integrals to their integration kernels. Taking the function F ( x) in Eq. (28) as an example, it's mapped to the symbol The symbols of the iterated integrals in Eq. (31) can be written as The symbols satisfy a lot of algebraic relations which are of great help to simplify the complicated expressions. For example where c is a constant. After simplification, it is possible to find an explicit functional representation for each symbol. In particular, when all the letters α k ( x) appearing in a given symbol are rational functions, it is straightforward to represent the function as polylogarithms or MPLs, which are well-studied and allow fast numerical evaluations. For example However, letters in the alphabet for Higgs boson pair production contain square roots where β 1 = µ, β 2 = ν and β 3 = −µν/(µ + ν). They make it challenging to convert the formal solutions to explicit functional forms. Fortunately, up to weight 2, only the first two kinds of square roots appear. In particular, there are only 4 kinds of symbols appearing at weight 2: where h(β i , β j ) is some function of β i and β j . The functional representation for the last symbol is simple: while for the first 3, we can get rid of the square roots with appropriate changes of variables. For the first two symbols, we use We work in the region 0 < z i < 1 such that the resulting functional representation is singlevalued. This corresponds to β i > 0. 3 We can then express z i in terms of β i as The expressions for β i < 0 can be found by analytic continuation. For the third symbol, we parameterize where we take 0 < y ij < x ij < 1 which corresponds to β i > 0 and β j > 0. The inverse relation is given by Now we can employ the algebraic properties of the symbols to further simplify the expressions. For example and similarly for the remaining two symbols. These symbols are simple enough, such that their functional representations can be found via direct integration. The results are We now turn to the weight-3 and weight-4 parts of the solution. These will involve the third square root in Eq. (36). Although it is still possible to find explicit functional forms from the symbols, it is often rather difficult [40]. Therefore, we write them as one-fold integrals over the weight-2 functions So far, we have discussed the solutions valid in the unphysical region. In practice, we need to do an analytic continuation to the physical region s > 2m 2 h . Up to weight 2, this can be simply done using the analytic expressions in Eq. (38) and (44), with the branch choice according to s → s + iδ and m 2 t → m 2 t − iδ. The treatment of the weight-3 and weight-4 parts is more tricky, since they are represented as one-fold integrals. We need to carefully deform the integration contour to avoid possible singularities. After these, we can numerically evaluate all the master integrals for topology E in the physical region. The results are shown in the next subsection.

Numerical results for topology E
In this subsection, we perform a numerical study of the 7-propagator two-loop integral I 1,1,1,1,1,1,1,0,0 in topology E. The purpose is to check how well the small-m h expansion can approximate the exact result. We calculate the exact result using the method of sector decomposition implemented in pySecDec [61]. We perform the small-m h expansion up to order m 4 h , which can be extended to higher powers of m h straightforwardly.
As in the one-loop case, we first fix √ −t = 200 GeV and show the value of the integral as a function of √ s in Figure 6. The upper two plots show the coefficient of ǫ −1 , and the lower two plots show the coefficient at ǫ 0 . We observe similar behaviors as the one-loop case: the small-m h expansion provides a good overall approximation to the exact result in the whole range of √ s, from the threshold region √ s 2m h , to the tt threshold √ s ∼ 2m t , and to the high energy regime √ s ≫ 2m t . There are exceptional values of √ s where the relative errors grow, which is due to that the value of the integral is close to zero. One should however not be worried since in these phase space points, this integral is not expected to be the dominant contribution. Similar behaviors have been observed in the one-loop case, as was shown in Figure 2 and 3. Even if there is a concern, one could easily add the order m 6 h terms which will further improve the accuracy of the approximation.   We further investigate the behavior of our approximation as a function of the transverse momentum p T of the Higgs boson. The invariant t is related to p T by where the ± sign corresponds to the forward and backward scatterings, respectively. For convenience, we only show the results with the + sign in the following. We take two typical values of the partonic center-of-mass energy: √ s = 500 GeV which is in the bulk region of the partonic cross section, and √ s = 1000 GeV which is in the high energy region. The corresponding numerical results are shown in Figure 7 and 8, respectively. At √ s = 500 GeV, we find that the approximation at order m 4 h works rather well for the real part of the integral, with per-mille accuracy in the whole range of p T . For the imaginary part, the accuracy is about 1%, and if one needs to have a better approximation, the order m 6 h terms should be added. When the center-of-mass energy goes higher, at √ s = 1000 GeV, the quality of the approximation becomes better, with per-mille accuracy in all situations. This can be expected since in the high energy region all the scales are much larger than m h . Finally, we stress that although in this subsection we only studied the behavior of a single integral, similar behavior is expected for the full amplitude. This has been verified at the one-loop level. At the two-loop level, this can only be done with the results for topology F, which is the subject of the next subsection.

Towards a solution for topology F
Topology F is the most difficult one as the differential equations for the master integrals cannot be transformed into a canonical form. The first place where this shows up is the 6-propagator sub-topology depicted in Figure 9 (which has been discussed in [59]). We denote the 6propagator master integrals as f (µ, ǫ), and collect the master integrals with fewer propagators in g(µ, ǫ), where µ = −4m 2 t /s as before. Then the differential equation satisfied by f (µ, ǫ) can  be written as The transformation to a canonical form amounts to get rid of the order ǫ 0 coefficient matrix B(µ) in the above equation. However, for the topology in Figure 9, we find that differential equations for the two top-level master integrals f (µ, ǫ) = ǫ(1 + 4ǫ) µ 2Ĩ 1,0,1,1,1,1,1,0,0 , ǫ µ 2Ĩ 1,0,1,2,1,1,1,0,0 involves the non-diagonal coefficient matrix at order ǫ 0 : (49) s Figure 9: Elliptic sub-topology for topology F. The thick lines represent massive propagators (top quarks), while the thin lines represent massless propagators (gluons). The dashed external legs are light-like, while the solid external leg is massive.
which cannot be transformed away. We call this situation as a "two-coupled" system of differential equations. In this case the solution necessarily involves elliptic integrals. To see that, we turn the system of two first-order differential equations into a second-order differential equation for f where the rational functions a(µ) and b(µ) are related to the matrix B(µ), and the function c (n) 1 (µ) depends on A(µ), C(µ, ǫ) and g(µ, ǫ). The homogeneous part of the above equation (i.e., with c (n) 1 (µ) absent) can be solved in terms of elliptic integrals, upon which the inhomogeneous part of the solution can be added.
We now turn to the integrals with 7 propagators in topology F, which to our knowledge were not discussed in the literature. There are 4 top-level master integrals in this case, and we are facing a four-coupled system of differential equations to begin with. In order to reduce the system to smaller blocks, we employ the method of [60]. Briefly speaking, the method goes as follows. We consider an N-coupled system of differential equations with respect to a single kinematic variable x (the extension to multiple variables is straightforward). We start from one of the master integrals I, and study the derivatives d k I/dx k . In general, not all the derivatives are linearly independent, and we can find the maximal number of independent derivatives I, dI/dx, . . . , d (r−1) I/dx (r−1) . If r < N, we can choose these r integrals as new master integrals such that they are decoupled from the remaining N − r master integrals.
Applying the above method, we find the basis of the top-level integrals in topology F h( x, ǫ) = Ĩ 1,1,1,1,1,1,1,0,−1 The differential equation takes the form where the order ǫ 0 coefficient matrix dB( x) takes the illustrative form in which " * " denotes non-zero entries. It is then clear that the 7-propagator master integrals can also be solved in terms of elliptic integrals. That said, the solutions are still rather complicated, and we leave them for future works.

Conclusion and outlook
In this paper we propose a new method to evaluate loop integrals where the masses of the internal particles are larger than the external particles. This can be applied to the pair production process and the associated production processes of the Higgs boson, which are mainly mediated by top quark loops. Our method amounts to perform a Taylor expansion in terms of the small masses of external particles. The coefficients of the expansion are written in terms of loop integrals with fewer mass scales than the original integrals, and are therefore easier to evaluate. The main difference between our method and other expansion methods lies in the fact that the validity of our expansion is not restricted to a special phase space region. Instead, our method provides a systematically improvable approximation in the entire phase space. Our expansion works particularly well for the high energy tails of kinematic distributions where many other expansions cease to be valid. We demonstrate our method using Higgs boson pair production as an example. At the leading order (one loop), we compare the approximate and the exact results both at the level of a single master integral and at the level of differential cross sections. We find that our method leads to rather good approximations in both cases. At the next-to-leading order (two-loop), we expand the amplitude and classify the resulting loop integrals into 2 planar topologies and 2 non-planar topologies. We reduce these integrals to master integrals using IBP reduction, and derive differential equations satisfied by the master integrals. We find that the equations for the 2 planar topologies and the non-planar topology E can be casted into the canonical form, they can be solved in terms of Chen iterated integrals. The solutions up to weight 2 can be written in terms of logarithms and dilogarithms, while the weight 3 and weight 4 parts of the solutions are given as one-fold integrals. We present numeric results for an integral appearing in the original amplitude (with non-zero external masses), comparing the exact values from sector decomposition and the approximate values from our expansion. We observe similar behaviors as in the one-loop case, that our expansion up to order m 4 h leads to good approximations in the whole phase space, which can still be further improved by incorporating terms suppressed by more powers of m h .
To construct the approximation to the full two-loop amplitude, we still need to calculate the master integrals in the other non-planar topology (topology F). The differential equations for the topology cannot be transformed into a canonical form. We reduce the system of differential equations into smaller blocks, and find that they can be solved in terms of elliptic integrals. The full solution and the numeric study for this topology will be presented in another work. Combining our efficient method to compute the two-loop amplitude with an infrared subtraction method for the real emission corrections, we expect to have a fast and reliable tool to calculate the differential cross sections for Higgs boson pair production. This will be useful for phenomenological studies and for the extraction of the Higgs self-coupling from future experimental data.