High-Energy Expansion of Two-Loop Massive Four-Point Diagrams

We apply the method of regions to the massive two-loop integrals appearing in the Higgs pair production cross section at the next-to-leading order, in the high energy limit. For the non-planar integrals, a subtle problem arises because of the indefinite sign of the second Symanzik polynomial. We solve this problem by performing an analytic continuation of the Mandelstam variables such that the second Symanzik polynomial has a definite sign. Furthermore, we formulate the procedure of applying the method of regions systematically. As a by-product of the analytic continuation of the Mandelstam variables, we obtain crossing relations between integrals in a simple and systematic way. In our formulation, a concept of"template integral"is introduced, which represents and controls the contribution of each region. All of the template integrals needed in the computation of the Higgs pair production at the next-to-leading order are given explicitly. We also develop techniques to solve Mellin-Barnes integrals, and show them in detail. Although most of the calculation is shown for the concrete example of the Higgs pair production process, the application to other similar processes is straightforward, and we anticipate that our method can be useful also for other cases.


Introduction
Perturbative quantum field theory has been describing particle physics phenomena very well, yet improving the perturbative series, i.e., calculating Feynman integrals, has been always challenging. It becomes more and more important to include higher order corrections to theoretical predictions as particle physics experiments, especially at the Large Hadron Collider (LHC), become more and more precise over the years. This means that a deeper understanding of higher order corrections is required to obtain meaningful theoretical predictions.
One of the most interesting problems in this field is higher order corrections to multiscale processes. The prime examples are the Higgs+jet, Higgs pair, and Higgs+Z production cross sections at the LHC. These 2 → 2 processes involve two kinematical variables, such as the center of mass energy and the transverse momentum, and the masses of internal or external particles serve as additional scales. Here, the word "multi-scale" is used when there are more than three scales. The bottlenecks of the calculation of multi-scale processes are the integration-by-parts (IBP) reduction and the evaluation of the resulting master integrals. The subject of this paper is the second issue, namely, the evaluation of multi-scale Feynman integrals.
Efforts to solve multi-scale Feynman integrals have persisted over the years. One of the milestones is the analytic computation of all the planar master integrals contributing to the Higgs → 3 partons process at two-loops [1]. Another milestone is the numeric evaluation of the Higgs+jet [2] and Higgs pair production [3,4] cross sections at two-loop level using the program SecDec [5,6]. An independent numerical evaluation of the Higgs pair production cross section is given recently [7]. It is worth mentioning some recent analytic calculations of three-scale four-point two-loop diagrams; the non-planar master integrals for µe scattering [8,9], the planar double box integral relevant to top pair production [10] and the planar master integrals relevant to di-photon and di-jet production [11]. These works show that even three-scale problems are difficult to solve. Recently, some of the non-planar master integrals for these processes in the limit m H = 0 have been solved [12], but there still remain unsolved master integrals.
It is a promising idea to reduce the number of scales entering integrals by expanding them in some small parameters. For a summary of this topic, see Ref. [13]. In this direction, the large-m t expansion of the Higgs+jet [14][15][16][17][18][19][20] and Higgs pair production [21][22][23][24][25][26][27][28] cross sections is very well investigated. However, it is not until recently that expansions in other parameters have been investigated. Concerning the Higgs+jet production cross section, the expansion in the small bottom quark mass m b m H [29,30] and in the small top quark and Higgs masses m t > m H [31] are performed. For the Higgs pair production cross section, the expansion in small m t for the planar master integrals [32] and in small Higgs transverse momentum [33] are performed. The rest of the master integrals of Higgs pair production in the small-m t expansion are obtained in Ref. [34] together with the results of this paper. Many of the non-planar master integrals are the same as, or related to those of Ref. [31] but we provide some new information needed for Higgs pair production. Furthermore, the method used in this paper-the method of regions-is completely different from the one in Ref. [31] at all steps of the calculation, so it provides a complementary understanding of the massive non-planer integrals. We would like to emphasize that the method of regions is a generic and systematic procedure to expand integrals, and thus the calculations shown in this paper can be applied to other integrals in a straightforward way.
The concept of dividing the domain of integration variables into several regions and expanding the integrand according to hierarchies in each region was introduced by Beneke and Smirnov [35]. The method is now called "expansion by regions" or "strategy of regions," and in this paper we call it the method of regions. A mathematical proof of the method of regions for a general integral is not yet known although many successful applications have been reported. In fact, the author of the most up-to-date textbook on this topic states in his book [36] The strategy of expansion by regions still has the status of experimental mathematics.
In the cases of off-shell large-momentum expansion and large-mass expansion, a mathematical proof based on a graph-theoretical language is known and it is called as "expansion by subgraphs" [13,37]. The procedure of expansion by subgraphs is implemented and can be performed in an automatic way [38,39]. The large-m t expansion mentioned above belongs to this category. A new proof of the method of regions was proposed by Jantzen [40] but its application is limited. The purpose of this paper is to show non-trivial examples where the method of regions works well, and our calculation shows the first application of the method to the high energy expansion of non-planar four-point integrals.
The remainder of the paper is organized as follows: in Section 2, we briefly summarize the method of regions. In Section 3 we introduce conventions, ideas, and techniques, which will be used in the following sections. In Section 4, 5 and 6, we apply the method of regions to the one-loop box diagram, the two-loop planar massive diagrams, and the two-loop nonplanar massive diagrams, respectively.

General Idea of the Method of Regions
The procedure of the method of regions is the following [13,35,36,40,41]: • Step 1: Assign a hierarchy to the dimensionful parameters.
• Step 2: Reveal the relevant scaling of the integration variable.
• Step 3: For each region, expand the integrand according to its scaling.
• Step 4: Integrate. Scaleless integrals such as ∞ 0 dx x a are set to zero.
• Step 5: Sum over the contributions from all the relevant regions.
The Step 2 is the crucial part of the method of regions, and an algorithm to reveal such scalings for a general integral is established based on the analysis of the convex hull [42,43]. One can use the algorithm, implemented in the Mathematica package asy2.1.m [43]. Although it is not proved that the algorithm works correctly for all the cases, no counterexample is known so far. Recently, a new idea to reveal relevant scalings is proposed based on the technique of power geometry, which is implemented in the Mathematica package ASPIRE [44]. In this paper we use asy2.1.m.
The practical bottleneck is Step 4, since the integration tends to be complicated even after the expansion if the original integral is very complicated. This is one of the reasons why testing the method of regions is difficult.
The method of regions was first applied to the momentum representation of the Feynman integrals, so the "regions" mean some domains of the loop momenta. Later, it was found that parametric representations such as the Feynman representation and the alpha representation are more convenient to apply the method [45]. Recently, it was proposed in Ref. [41] to use yet another parametric representation, Lee-Pomeransky representation [46], to apply the method of regions. For all the representations mentioned above, one has to follow Step 1 to 5 for practical calculation.

Conventions
We distinguish the exact equal sign " = " and the equal sign under a certain analytic continuation. For this purpose, we introduce a sign " AC = " and use it as, e.g., where i0 represents an infinitesimal positive imaginary number. We interpret log(z) as the principal value of the complex logarithm whose imaginary part lies in the interval (−π, π]. Both the left-hand side and the right-hand side of Eq. (3.1) are well-defined in the entire domain of z, but the equality is valid only in the upper half plane of z. This is how analytic continuation is performed, and that is why we add "AC" to the normal equal sign in Eq. (3.1). The equality of a series expansion like is in principle also regarded as an analytic continuation. However, when a hierarchy like m M is explicitly stated in the text, we use normal equal sign. We use a simplified expression of the Landau O notation for more than one variable as The Euler-Mascheroni constant is denoted as γ E . We use the alpha representation to calculate Feynman integrals. The integration measure and the analytic regularization parameters are defined as .
The analytic regularization parameters δ j play one essential role and three secondary roles: (i). The essential role is to regularize the contribution of individual regions which are divergent if we naively expand in α i . This means that individual contributions are regulator dependent, and the dependence on δ j cancel after we sum all the contributions and take the limits δ j → 0. In taking the limit, it is necessary to specify the order because some of them do not commute. We express the sequence of limits as lim ,δn,...,δ 2 ,δ 1 →0 (iii). By shifting δ j → δ j + 1, we can express polynomials of α i in the integrand. For example, when n = 2, This property is usually used to express the integrals with higher powers of propagators.
(iv). We use the property of Eq. (3.6) to express the higher order terms. [See the text below Eq. (4.20).] The sum of the variables will be expressed as The bar on an index indicates that the variable corresponding to the index is subtracted instead of added. Sometimes and δ j are treated in the same way, and in those cases we express as δ 0 . For example, δ 0012 = 2 + δ 1 − δ 2 . Also, we introduce the following compact notation for the product of Γ-functions In Step 3 of Section 2, we expand the integrand of Feynman integrals in terms of soft parameters. In order to control the expansion in a systematic way, we introduce an auxiliary soft-scaling parameter χ. For example, assume that we have four variables α 1 , ..., α 4 whose scalings are where m ∼ χ is a soft parameter and M ∼ 1 is a hard parameter. In this case, we apply a substitution to the integrand and expand in χ. After that, we can set χ = 1. In this paper we denote the scalings (3.9) as or simply (1, 0, 1, 0).

Kinematics and High Energy Expansion
The assignment of the external momenta q 1 , ..., q 4 is illustrated in Fig. 1. We consider the 2 to 2 process but define all the external momenta as incoming. In addition to the usual Mandelstam variables s, t, u, (which we call physical Mandelstam variables), we introduce S, T, U as In our calculation we assume that S, T, and U are positive and thus call them positive Mandelstam variables. Sometimes two of the three Mandelstam variables are sufficient to express four-point functions, and indeed for planar integrals we do not use U [See Section 5]. However for non-planar integrals, all of S, T, U are required to make the second Symanzik polynomial positive [See Subsection 6.1].
In this paper, we consider the master integrals of Higgs pair production at next-toleading order, where the loops are induced by the top quark. Therefore there are two mass scales, the Higgs mass m H and the top quark mass m t . We consider the high energy limit where the following hierarchy is satisfied The high energy expansion in this case is two-fold. First, we treat m 2 H as the soft parameter and m 2 t , S, T, U as the hard parameters. Afterwards, we treat m 2 t as the soft parameter and S, T, U as the hard parameters.
In the first expansion, i.e. the m H -expansion, m H enters the integrals through the on-shell condition of the external momenta (3.14) In terms of χ, the scaling we impose here is 15) and the resulting series expression of an integral I is expressed symbolically as 1  In the second expansion, i.e. the m t -expansion, the on-shell condition of the external momenta becomes In terms of χ, the scaling we impose here is 18) and the resulting series for an integral I is expressed symbolically as The result should be expressed in a way suitable for the evaluation with the physical kinematics. In order to achieve that, the analytic continuation is applied, where i0 represents an infinitesimal positive imaginary number, (note that the massless on-shell condition (3.17) is adopted), and 0 ≤ v ≤ 1 in the physical kinematics. After expressing the result in terms of s and v, the expression is unique. The analytic continuation of T, U is trivial since their signs are consistent with those of the physical kinematics. The results are expressed in terms of harmonic polylogarithms (HPL) [47] and we introduce an abbreviation for HPL as and so on. The argument of the HPL is always chosen to be v. For example, We use the Mathematica package HPL.m [48,49] in dealing with HPL. When diagrams of the same topology but with a different assignment of external momenta are considered, it is convenient to relate them by applying replacements of the external momenta, which means replacements of the Mandelstam variables. We call these -7 - In the second expansion, i.e. the m t -expansion, the on-shell condition of the external momenta becomes In terms of χ, the scaling we impose here is 18) and the resulting series for an integral I is expressed symbolically as The result should be expressed in a way suitable for the evaluation with the physical kinematics. In order to achieve that, the analytic continuation is applied, where i0 represents an infinitesimal positive imaginary number, (note that the massless on-shell condition (3.17) is adopted), and 0 ≤ v ≤ 1 in the physical kinematics. After expressing the result in terms of s and v, the expression is unique. The analytic continuation of T, U is trivial since their signs are consistent with those of the physical kinematics. The results are expressed in terms of harmonic polylogarithms (HPL) [47] and we introduce an abbreviation for HPL as and so on. The argument of the HPL is always chosen to be v. For example, We use the Mathematica package HPL.m [48,49] in dealing with HPL. When diagrams of the same topology but with a different assignment of external momenta are considered, it is convenient to relate them by applying replacements of the replacements "crossing relations". This subject is already well-established in the case of HPL [50], but we propose a simpler way to obtain the crossing relations. In Fig. 2, the commutative diagram of the crossing and analytic continuation is given. Usually an integral is given as the bottom-left expression, where the result is expressed in terms of the physical kinematic variables. On the other hand, we proceed the crossing in the upper expression, where the result is expressed in terms of the positive Mandelstam variables. The analytic continuation in the upper expression is the simple replacement of the positive Mandelstam variables, whereas the analytic continuation of the bottom expression requires the precise knowledge of the branch cuts. We take the approach (a) of Fig. 2 because it is easy to implement in programs, and crosscheck the result using approach (b). One can interpret the simplification by introducing the positive Mandelstam variables as the resolution of singularities by increasing the dimension, or in another words, we lose some information when we map the three-variable function F into the two-variable function f . We would like to emphasize that the method to obtain the crossing relations explained above is a by-product of introducing the positive Mandelstam variables. The most important point in introducing the positive Mandelstam variables is that it makes the Symanzik polynomials positive and thus allows us to apply the method of regions safely. [See Subsection 6.1.]

A First Example: One-Loop Box Diagram
We consider the massive one-loop Feynman integral family where the infinitesimal negative imaginary part of each denominator is implicit. The external momenta q i satisfy the on-shell conditions (3.14). We consider the box diagram shown -8 - external momenta, which means replacements of the Mandelstam variables. We call these replacements "crossing relations". This subject is already well-established in the case of HPL [50], but we propose a simpler way to obtain the crossing relations. In Fig. 2, the commutative diagram of the crossing and analytic continuation is given. Usually an integral is given as the bottom-left expression, where the result is expressed in terms of the physical kinematic variables. On the other hand, we proceed the crossing in the upper expression, where the result is expressed in terms of the positive Mandelstam variables. The analytic continuation in the upper expression is the simple replacement of the positive Mandelstam variables, whereas the analytic continuation of the bottom expression requires the precise knowledge of the branch cuts. We take the approach (a) of Fig. 2 because it is easy to implement in programs, and crosscheck the result using approach (b). One can interpret the simplification by introducing the positive Mandelstam variables as the resolution of singularities by increasing the dimension, or in another words, we lose some information when we map the three-variable function F into the two-variable function f . We would like to emphasize that the method to obtain the crossing relations explained above is a by-product of introducing the positive Mandelstam variables. The most important point in introducing the positive Mandelstam variables is that it makes the Symanzik polynomials positive and thus allows us to apply the method of regions safely. [See Subsection 6.1.]

A First Example: One-Loop Box Diagram
We consider the massive one-loop Feynman integral family where the infinitesimal negative imaginary part of each denominator is implicit. The external momenta q i satisfy the on-shell conditions (3.14). We consider the box diagram shown in Fig. 3, J 1,1,1,1 and its alpha representation is where the first Symanzik polynomial U and the second Symanzik polynomial F are given by We make clear the analytic continuation in Eq. (4.2) because the right hand side is regularized by δ j whereas Eq. (4.1) is explicitly δ j -independent. Also, we assume m 2 H < 0 in order to ensure the convergence of the integral and perform the analytic continuation of the result to m 2 H > 0 at the end, which turns out to be trivial.

Expansion in the Higgs Mass
We first expand in m H . By using asy2.1.m, we find that there is only one relevant scaling The expansion corresponding to this scaling is actually just the Taylor expansion of the integrand in m 2 H , and the original integral can be written as In particular, the leading term is identical to the box diagram with completely massless external legs. The fact that the expansion in m 2 H and the integration commute is reasonable because m H in the denominator of the integrand (appearing through q 2 3 = m 2 H ) is always accompanied by m t which regulates the integral and the limit m H → 0 always exists.
In collider physics, it is common to use s and the transverse momentum p T as the kinematic variables, where t and p T are related as in Fig. 3, J 1,1,1,1 and its alpha representation is where the first Symanzik polynomial U and the second Symanzik polynomial F are given by We make clear the analytic continuation in Eq. (4.2) because the right hand side is regularized by δ j whereas Eq. (4.1) is explicitly δ j -independent. Also, we assume m 2 H < 0 in order to ensure the convergence of the integral and perform the analytic continuation of the result to m 2 H > 0 at the end, which turns out to be trivial.

Expansion in the Higgs Mass
We first expand in m H . By using asy2.1.m, we find that there is only one relevant scaling The expansion corresponding to this scaling is actually just the Taylor expansion of the integrand in m 2 H , and the original integral can be written as In particular, the leading term is identical to the box diagram with completely massless external legs. The fact that the expansion in m 2 H and the integration commute is reasonable because m H in the denominator of the integrand (appearing through q 2 3 = m 2 H ) is always accompanied by m t which regulates the integral and the limit m H → 0 always exists.
In collider physics, it is common to use s and the transverse momentum p T as the kinematic variables, where t and p T are related as Since these relations are m H -dependent, one should be careful when expanding in m H . In order to clarify the point, let us consider two expressions which are related by Eq. (4.7). We would like to analyze the cross section for fixed p T but not t, so let us consider the m H -expansion of f 2 (p 2 T , m 2 H ): On the other hand, the kinematic variable appearing in the Feynman integral is t and the natural representation is f 1 (t, m 2 H ). Thus, we express the ingredients of Eq. (4.9) in terms of f 1 (t, m 2 H ) as where t 0 = t| m H =0 . Apparently, Eq. (4.9) becomes complicated when Eqs. (4.10), (4.11), (4.12) are substituted. However, taking into account the m H -expansion of f 1 (t, 0) the expression becomes simpler where the m H -dependent t (4.7) is used. In general, for a given order of m 2 H , the difference between the strict m H -expansion of f 2 (p 2 T , m 2 H ) for fixed p T and the m H -expansion of f 1 (t, m 2 H ) for the m H -dependent t is higher order in m 2 H . The conclusions of this subsection are the following: (i). As the result of the m H -expansion, the integrals reduce to integrals with massless external legs.
(ii). The m H -expansion for fixed p T is obtained by the m H -expansion with fixed t, keeping the m H -dependence of t.

Expansion in the Top Quark Mass
After the expansion in m H , we have integrals which depend on m t , S, and T . We now consider the expansion in m t assuming the hierarchy (3.18). The integral of interest is The scalings of regions 2 to 5 reflect the symmetries of the integral, α 1 ↔ α 3 and α 2 ↔ α 4 . Eq. (4.15) is thus expressed as the sum of the contributions from these five regions:
(b). There is only one soft parameter in the all-hard region, which is m t .
(c). The contribution from the all-hard region can be expressed as the massless integral of the original topology. In particular, the leading order term is obtained by substituting m t = 0 into the original integral.
(d). The leading order contribution of the all-hard region is O(χ 0 ).
(e). The contribution from the all-hard region has no singularities in δ j .
Because of these properties, the contribution from the all-hard region can be calculated in two ways. The first is the procedure universal for any region, and the second is to use the momentum representation. We show them in order. First, we show the universal procedure. By expanding Eq. (4.15) in terms of χm 2 t , we obtain the contribution of this region as where U (1) = α 1234 and F (1) = Sα 1 α 3 + T α 2 α 4 . As stated above, the leading order term is the massless box integral and we name it T (1) δ 1 ,δ 2 ,δ 3 ,δ 4 , for later use: (4.20) The integrand of the higher order corrections has additional factors of α i which can be expressed by some shifts of δ j → δ j + 1. Indeed, Eq. (4.18) can be expressed as where P n x = Γ(x + n)/Γ(x) is the Pochhammer symbol. We call T (1) δ 1 ,δ 2 ,δ 3 ,δ 4 , a "template integral" since all the higher order terms can be expressed in terms of T (1) Since there is no singularity in δ j , we can safely set δ j = 0 in Eq. (4.21). Then, the leading order term is expressed as where the technique explained in Appendix A is used to set the integration contour to a straight line. The expression at → 0 is obtained by using the package MB.m [51], and the result is The higher order terms are given by T 1,0,0,0, , T 0,1,0,0, etc, and can be calculated in a similar way.
The other procedure to calculate the contribution of the all-hard region is the following. We return to the momentum representation (4.1) and expand each propagator in m t as Then, the contribution can be expressed in terms of the integral family as  The IBP-reduction of massless integrals is computationally easy even at the two-loop level, and thus very useful. In the calculation of two-loop integrals, we adopt this approach to calculate the contribution from the all-hard region.

Regions 2, 3, 4, 5
The contribution of Region 2 is obtained by applying the second scaling of Eq. (4.17) and expanding in χm 2 t , χα 3 , χα 4 , where U (2) = α 12 and F (2) = m 2 t α 12 U (2) + Sα 1 α 3 + T α 2 α 4 . The integration over α 1 , . . . , α 4 can be performed using the relation (A.3) and a variant of it, and the template integral of this region is The higher order terms in Eq. (4.28) which contain the inverse of U (2) can be expressed by the shift → − 1 in the template integral. Thus Eq. (4.28) is written as Recall that P n x = Γ(x + n)/Γ(x) is the Pochhammer symbol. As mentioned in Subsection 3.1, the result of the limits δ j → 0 depends on the order in which we take them. For example, when we take the sequence of limits with ascending values of j, we obtain whereas with descending values of j, we instead obtain (4.32) The order dependence is not problem, provided we use the same order throughout the calculation. The artifacts caused by the δ j will cancel after we sum the contributions from all the relevant regions.
Due to the symmetries of the diagrams, the template integrals of the other regions can be expressed in terms of T and when we take the ascending order of limits, we obtain . As mentioned, the result is δ j -independent. There are 24 possible ways to order δ 1 , δ 2 , δ 3 , δ 4 in taking the limit, and we have confirmed that the result is the same for all of the orderings.
Since the original integral is finite in the limit → 0, the poles of in the individual contributions from each region cancel.

Higher Order Terms in m t
The method of regions can be used to obtain a series expansion up to arbitrary order in the soft parameters, m H or m t . For example, Eqs. (4.21), (4.30) includes the contributions up to O(χ). In general, the higher order terms can be expressed in terms of the template integrals. Therefore, in principle, it is straightforward to calculate higher order terms up to arbitrary order, once the template integrals have been obtained. However, the number of integrals to calculate increases rapidly with the order of χ, and this makes it hard to calculate higher order terms. Especially in the two-loop case, some of the template integrals contain multi-dimensional Mellin-Barnes integrals, which are not so easy to solve. Furthermore, the number of integrals increases more rapidly than in the case of the one-loop calculation. Therefore it is better to use another method to calculate the higher-order corrections.
The use of differential equations solves this issue [29,31,32,34]. We use the differential equation with respect to m 2 t to obtain higher order corrections in m 2 t . Since we know that the integral has the form Eq. (4.15) = the set of differential equations reduces to a set of linear relations of c n 1 ,n 2 (S, T ) which simplifies the problem a lot. In this sense, the leading order terms which we calculate in the previous subsections play the role of the boundary conditions of the differential equations.

Integrals with Fewer Lines
Once we have calculated the box integral, there are several shortcuts to calculate integrals with fewer lines such as the triangle integral and the self energy integral. Let us consider the s-channel triangle diagram, J 1,1,1,0 , as an example. The alpha representation of J 1,1,1,0 is obtained by setting α 4 → 0 in Eq. (4.3), since the forth propagator is absent in J 1,1,1,0 . If we use asy2.1.m to reveal the relevant scalings for J 1,1,1,0 , we obtain (0, 0, 0), (0, 0, 1), (1, 0, 0), (4.41) however, we do not have to do that. We do not have to derive the template integrals for J 1,1,1,0 because they are derived from the template integrals of J 1,1,1,1 .
Using the fact that the δ j -dependence of the alpha representation is expressed by the replacement of a i → 1 + δ j , the triangle integral J 1,1,1,0 can be expressed by the limit δ 4 → −1. Therefore by taking the limit δ 4 → −1 to the template integral of J 1,1,1,1 , one can obtain the template integral for J 1,1,1,0 . For example, 42) and this is the template integral for the region (0, 0, 1). Sometimes the limit vanishes due to a suppression factor 1/Γ(1 + δ 4 This fact is reasonable because the number of relevant regions for J 1,1,1,0 is smaller than for J 1,1,1,1 . Two of the four soft template integrals are non-vanishing after taking the limit, and they are the two soft template integrals of J 1,1,1,0 .

Two-Loop Planar Diagrams
We consider the following Feynman integral families J PL1 a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,a 7 = where the momenta of the lines are Recall that q 13 = q 1 + q 3 . We consider the integrals I PL1 = J PL1 1,1,1,1,1,1,1 and I PL2 = J PL2 1,1,1,1,1,1,1 whose diagrammatic representations are shown in Fig. 4. Conceptually there is no difference between the procedure of applying the method of regions to these integrals and the example discussed in Section 4. In particular, the property of the F-function that it is positive definite in the u-channel is the same [cf. the text below Eq. (4. 16)]. The only new ingredient is that now the template integrals are expressed by at most two-dimensional Mellin-Barnes integrals, which are not trivial to solve. However, their calculation is a subset of the calculation of the non-planar integrals, so we do not describe it here [cf. Subsection 6.4]. Therefore we briefly summarize the important ingredients of the two-loop planar integrals in this section.

Double Massive Box Diagram
By using the package asy2.1.m, we reveal thirteen relevant scalings: The template integrals of these regions are summarized in Appendix B, and can be found in the ancillary file to this paper [52]. The result of this integral at the leading order in m t , up to O(ϵ), is given by where the coefficients d n 1 ,n 2 are given by (5.9) -17 - The template integrals of these regions are summarized in Appendix B, and can be found in the ancillary file to this paper [52]. The result of this integral at the leading order in m t , up to O( ), is given by where the coefficients d n 1 ,n 2 are given by 90 + 24h 0 h 1 ζ 3 + 24iπh 1 ζ 3 − 56h 2 ζ 3 + 4(4h 23 + 7h 32 + 21h 5 − 13ζ 5 ) + 6π 2 ζ 3 . (5.9)
The template integrals of these regions are summarized in Appendix B, and can be found in the ancillary file to this paper [52]. The result of this integral at the leading order in m t , up to O( ), is where the coefficients d n 1 ,n 2 are now given by

Relevant Scaling
A new feature appears in Eq. (6.7): if we impose the relation S + T + U = 0, we obtain and the sign of F becomes indefinite. In such a case, it is not guaranteed that the method to reveal the relevant scalings works properly [43]. An idea to solve this problem is to perform a proper change of variables and decompose the integration domain such that F is positive-definite [43]. However in our case, this approach does not resolve the indefinite sign of F since there is no simple change of variables to make F positive-definite. The solution to this problem is to keep S, T, U as independent variables. It is obvious that Eq. (6.7) is positive definite in this case, and we can apply the method of regions, expand the integrand, and express the result in terms of Mellin-Barnes integrals in terms of the positive Mandelstam variables. The procedure to obtain expression (6.7) is the following: we first compute F respecting the original definition of the Mandelstam variables (3.12). At this point, there are some redundant terms in F such as (S +T +U )α 2 α 3 α 6 . We minimize the number of terms, under the condition that F remains positive definite. The resulting F is unique.
There are two commands in the package asy2.1.m to reveal the relevant regions. The first is AlphaRepExpand, which accepts a set of propagators and replacement rules as input. The other is WilsonExpand, which accepts the Symanzik polynomials as input. Here we must use WilsonExpand since the conventional routine used in AlphaRepExpand either eliminates U completely or keeps U completely, whereas we want to eliminate U partially, as explained above. There is an option Preresolve in AlphaRepExpand which makes it attempt some changes of variables to make F positive definite, but in our cases this option did not solve the problem.
With this setup and the hierarchy (3.18), we obtain the following fourteen relevant

Relevant Scaling
A new feature appears in Eq. (6.7): if we impose the relation S + T + U = 0, we obtain 8) and the sign of F becomes indefinite. In such a case, it is not guaranteed that the method to reveal the relevant scalings works properly [43]. An idea to solve this problem is to perform a proper change of variables and decompose the integration domain such that F is positive-definite [43]. However in our case, this approach does not resolve the indefinite sign of F since there is no simple change of variables to make F positive-definite. The solution to this problem is to keep S, T, U as independent variables. It is obvious that Eq. (6.7) is positive definite in this case, and we can apply the method of regions, expand the integrand, and express the result in terms of Mellin-Barnes integrals in terms of the positive Mandelstam variables. The procedure to obtain expression (6.7) is the following: we first compute F respecting the original definition of the Mandelstam variables (3.12). At this point, there are some redundant terms in F such as (S +T +U )α 2 α 3 α 6 . We minimize the number of terms, under the condition that F remains positive definite. The resulting F is unique.
There are two commands in the package asy2.1.m to reveal the relevant regions. The first is AlphaRepExpand, which accepts a set of propagators and replacement rules as input. The other is WilsonExpand, which accepts the Symanzik polynomials as input. Here we must use WilsonExpand since the conventional routine used in AlphaRepExpand either eliminates U completely or keeps U completely, whereas we want to eliminate U partially, as explained above. There is an option Preresolve in AlphaRepExpand which makes it attempt some changes of variables to make F positive definite, but in our cases this option did not solve the problem.
Region 1 (0, 0, 0, 0, 0, 0, 0) (all-hard region) As shown in Subsection 4.2, the contribution from the all-hard region can be expressed in terms of massless integrals of the same topology. The massless integral that is relevant here has been calculated in Ref. [53].

Template Integrals for Regions 2 to 14
The template integral of each region is expressed as where 2 ≤ j ≤ 14. For simplicity, we omit the subscripts of T (j) δ 1 ,δ 2 ,δ 3 ,δ 4 ,δ 5 ,δ 6 ,δ 7 , and write T (j) when they are in the ordinary order. F (2) = m 2 t α 45 U (2) + S(α 2 α 3 α 5 + α 1 α 45 α 7 + α 4 α 5 α 7 ) + T α 1 α 3 α 6 + U α 2 α 4 α 6 . (6.12) The integration in α 7 can be performed using the relation (A.3). Then, we perform the following change of variables 13) and the template integral becomes This may be because there is no cancellation between two hard-scaling terms resulting in a soft-scaling term in Eq. (6.8). However, this observation is made in hindsight, since in principle it is not guaranteed that the regions are found correctly. One could also use ASPIRE [44] instead of asy2.1.m and obtain the correct scalings, if one similarly chooses suitable values for S and T .
There is another representation of T (13) which is more suitable if we require calculation beyond order 0 : . (6.43) This representation consists of two integrals with at most two kinematic parameters, and their calculation is simpler. It has some disadvantages, however; it holds only when δ 3 and δ 4 are non-negative integers, and Eq. (6.43) is shown in the special case that δ 3 = 0, δ 4 = 0, since that is the typical situation (here we may set these values because T (13) is regular in δ 3 and δ 4 ). Another disadvantage is that each integral produces singularities in which cancel in their sum, however as a by-product, higher-order derivatives of Γ-functions contribute to the 0 order.

Analytic Continuation
As mentioned in the text below Eq. (A.1), the Mellin-Barnes integrals in the template integrals are assumed to be regularized by choosing suitable values of δ j . However, the quantity we need is the one where δ j → 0 for all j. Therefore we need to analytically continue the Mellin-Barnes integrals in terms of δ j . We describe the procedure of the analytic continuation showing the case of Region 7, which is one of the most involved cases, as a concrete example. 5 We take the limit of the ascending order of δ j , lim ,δ 7 ,δ 6 ,δ 5 ,δ 4 ,δ 3 ,δ 2 ,δ 1 →0 T (7) . (6.44) As mentioned below Eq. (A.1), in general the integral contours of z 1 , ..., z 4 are assumed to be straight lines parallel to the imaginary axis [cf. Fig. 6 (b)]. For Region 7, we have only z 1 , z 2 as integration variables. We choose Re(z 1 ) = −1/5, Re(z 2 ) = −1/3. Then, we may choose , δ 5 = 1 500 , δ 6 = 17881 15540 , δ 7 = − 64003 62160 , = − 2507 10360 (6.45) to regularize the template integral T (7) . 6 We try to set as many parameters as possible to zero in ascending j order in Eq. (6.45). The limit of δ 1,2 → 0 in Eq. (6.44) is now trivial. The analytic continuation of δ 3 from −13/30 to 0 makes the first rightmost left poles of Γ(z 1 + δ3 6 − 1), Γ(z 2 − δ 3 ), and Γ(z 12 − δ 03 ) into right poles, which must be compensated by adding their residues. This analytic continuation procedure is automatized in the Mathematica package MB.m [51]. After the analytic continuation in terms of δ 3 , the integral depends on δ 4 , ..., δ 7 and , and we repeat the same procedure for δ 4 , then, δ 5 and so on. In this way, we obtain a combination of integrals for which the arguments of the Γ-functions in the integrand contain only z 1 and z 2 such as The methods to solve these integrals are explained in the next subsection.

Solving the Mellin-Barnes Integrals
The usual idea to solve the Mellin-Barnes integral is to apply the first and the second Barnes lemma and variants of them. The Mathematica package barnesroutines.m [54] performs this procedure in an automatic way, and solves some of the Mellin-Barnes integrals we encounter. Unfortunately, not all of them are solved by this package, and we describe here how to treat such cases. The essential points are mentioned in Ref. [32], and we fit or extend them to our integrals.

Three-and Four-Dimensional Mellin-Barnes Integrals
The template integral of the Region 12 (6.37) is expressed as a four-dimensional Mellin-Barnes integral. Thus, the contribution from Region 12 contain a four-dimensional integral of the form to reduce the number of the Γ-functions whose argument contains z 4 from 6 to 4. Now one can apply the first Barnes lemma to solve the z 4 -integral. The resulting three-dimensional integral can be easily reduced to a sum of two-dimensional integrals since the z 3 -integral can be solved by the variants of the first and second Barnes lemmas. Thus, we have two-dimensional integrals, and the way to solve them will be explained below.
Note that the reduction in Eq. (6.48) can be done only after the limits δ j → 0 and → 0 since the Γ-functions of the denominator and the numerator have different dependence on , thus the cancellation does not occur before the limits have been taken. The content of this subsection is not formulated in an algorithmic way and has been done manually.

Two-Dimensional Mellin-Barnes Integrals
In the cases of integrals with no argument such as 49) or integrals with a single argument of the form where 3 F 2 is the generalized hypergeometric function [55,56]. A useful corollary of Eq. (6.51) is presented in Appendix C. The resulting one-dimensional integrals can be solved by the method below. Integrals with two different arguments are difficult to solve. However in our case, such integrals only appear at higher orders in so we do not need to consider them.

One-Dimensional Mellin-Barnes Integrals
For integrals with no argument such as 7 we evaluate them numerically using Mathematica and apply the PSLQ algorithm [57,58] to fit them to a basis of constants which consists of all possible products of {1, γ E , π 2 , ζ 3 , ζ 5 } (6.53) up to a transcendental weight of five. The results turn out to only require constants to weight four. Typically, 50-70 digits of the numerical result are sufficient to obtain the correct answer, which we verify with 100 more digits.
The integrals with one argument typically have the form Various combinations of X 1 , X 2 appear since the template integrals contain them. We obtain the series expansion of Eq. (6.54) by taking the residue of the left poles or the right poles. By adjusting which poles we consider, we can choose the series in terms of either T /S, U/S, or T /U : 8 (log X) n 1 X n 2 , (6.55) Now we apply analytic continuation and obtain Recall that h 1 = − log(1 − v). We fit the series with HPL and express the result in terms of h N . In the case of (6.56), the series in v is directly fit to h N . In the case of (6.57), we first fit the series to HPL with the argument of (1 − v) and then express them in terms of h N . In the case of (6.58), we first fit the series to HPL with the argument of v/(1 − v) and then express them in terms of h N . Taking into account that 0 ≤ v ≤ 1 and the brach cut of h N >0 lies on the real axis of v > 1, we never cross the branch cut in the above procedure. The information of the branch cut is encoded in the analytic continuation of log X.
We already cover all of the combinations of X = X 1 /X 2 , so the calculation of the crossed diagrams can be done with the same procedure. For our sample integral (6.1), there are about 50 one-dimensional Mellin-Barnes integrals which are treated in this way.

Combining the Results
Summing the contributions from all the relevant regions, we obtain for our sample integral This result is consistent with Ref. [31] after a proper analytic continuation. One remarkable feature of the result is that it contains terms proportionals to 1/m t . Higher m t -order correction also contain odd-power terms. These odd-power terms come from Region 2 and 3: lim ,δ 7 ,δ 6 ,δ 5 ,δ 4 ,δ 3 ,δ 2 ,δ 1 →0 where ∆(1/δ j ) has poles in δ j , and these poles are cancelled by the contributions from the other regions.

Six-Line Integrals
We can use the method described in Subsection 4.4 to compute the integral with fewer lines. For example, J NPL1 1,1,1,1,1,1,0 is obtained by shifting δ 7 → δ 7 − 1 and repeat the same procedure described above. Note that the analytic continuation of δ j may change due to the shift. For example, the template integral of Region 12 (6.37) can be regularized with δ i>0 = 0, whereas J NPL1 1,1,1,1,1,1,0 requires a non-zero δ 7 to regularize that template integral.

Summary
Asymptotic expansion is useful to extract information from multi-scale Feynman integrals, which are difficult to solve exactly. The method of regions plays an essential role in this extraction. The crucial part of the method of regions is to reveal the relevant regions correctly, and a naive application of the conventional method fails in the case of nonplanar integrals for which the second Symanzik polynomial does not have a definite sign. We solve this problem by performing an analytic continuation of the Mandelstam variables such that the second Symanzik polynomial is positive definite. We show the applicability of the method of regions by the explicit calculation of the master integrals of the Higgs pair production cross section at two-loop order, in the high energy limit. It is straightforward to extend our calculation to other four-point two-loop integral which satisfy q 2 i m 2 S, T, U where q i are the external momenta and m is the mass of the internal lines. We anticipate that our idea to make the second Symanzik polynomial positive definite works in more general cases.
In addition to solving the issue of the sign of the second Symanzik polynomial, we formulate the procedure of the calculation in a systematic way. The contribution from each region is expressed in terms of Mellin-Barnes integrals, and a way to solve them is presented. The procedure presented here to solve the Mellin-Barnes integrals beyond the Barnes lemmas is not applicable to the general case, although it is sufficient to solve our master integrals completely. The automatization of this part is a future project.
As a by-product of introducing the positive Mandelstam variables, it becomes easier to obtain the crossed integrals, since the crossing of the positive Mandelstam variables does not cross any branch cut.
We compute the first few terms of the series in m t , and the higher order terms can be obtained by the use of the m t -differential equations. In this sense, our results can be considered as the boundary conditions of the differential equations with respect to m t .

A Mellin-Barnes Integrals
Here we explain elementary aspects of the Mellin-Barnes integrals. The basic tool is the identity where the integration contour C satisfies the following three properties: • C runs from −i∞ to +i∞.  are always satisfied. We choose δ j to ensure the condition λ ̸ ∈ {0, −1, −2, −3, ...}. 9 When Re(λ) > 0, it is possible to set C as a straight line along the imaginary axis [See Fig. 6 (b).] The integrand converges rapidly to zero at Im(z) = ±∞ rapidly so that we can shift the end points of the integration. It is sometimes difficult to fix the contour C immediately when the Mellin-Barnes representation is introduced because additional convergence conditions concerning the α iintegration should be taken into account. For example, the relations After the integration with respect to all α i is solved, we deform C to a straight line by tuning δ j to separate the right poles and the left poles as in Fig. 6 (b). We obtain an expression which is valid for certain values of δ j . In order to take the limit δ j → 0, it is necessary to apply a final analytic continuation of δ j from their initial values into a neighborhood around 0. This time, we allow poles to cross C. For this procedure, we use the package MB.m [51]. 9 If some of the left poles and right poles merge for any choice of δj, it is necessary to compensate the contributions of the merged poles by adding or subtracting the residue of the poles. However, it turns out that the merger of poles does not happen in our cases.
-32 - The integrand converges rapidly to zero at Im(z) = ±∞ rapidly so that we can shift the end points of the integration. It is sometimes difficult to fix the contour C immediately when the Mellin-Barnes representation is introduced because additional convergence conditions concerning the α iintegration should be taken into account. For example, the relations require Re(a) > −1 and Re(a+b) < −1, and if a or b contains the Mellin-Barnes integration variable z, C should be chosen accordingly. This kind of convergence condition appears every time an integration is performed, and it can happen that an analytic continuation of the regularization parameters is necessary between two successive integrations. Therefore it is necessary to deform C smoothly in each step of the integration assuming a flexible contour like Fig. 6 (a) in order to ensure convergence. Note that after the integration, the analytic continuation of Γ-functions in the left hand sides of Eqs. (A.3), (A.4) is possible. After the integration with respect to all α i is solved, we deform C to a straight line by tuning δ j to separate the right poles and the left poles as in Fig. 6 (b). We obtain an expression which is valid for certain values of δ j . In order to take the limit δ j → 0, it is necessary to apply a final analytic continuation of δ j from their initial values into a neighborhood around 0. This time, we allow poles to cross C. For this procedure, we use the package MB.m [51].
Consider the case where a 1 = n 1 − b 1 , a 2 = n 2 − b 2 , n 1 , n 2 ∈ N, where N is the set of positive integers. Note that n 1 = 0, n 2 = 0 because otherwise the left pole and the right pole merge. Substituting Eq. (C.2) into Eq. (C.1), the arguments of 3 F 2 become We can apply the partial fraction decomposition further and obtain Here the equation in the case of x 1 − x 2 ∈ Z is shown, but the case of x 1 − x 2 ∈ Z is similar. It is also possible to set x 2 → x 1 + n, n ∈ Z at the end. The infinite sum of m in Eq. (C.8) is now possible, and we obtain (C.9) Substituting the replacement x 1 → −b 2 +b 3 +n 2 , x 2 → −b 1 +b 3 +n 1 , x 3 → b 3 −c+n 1 +n 2 into Eq. (C.9), we obtain (C.10) In this expression, there are only finite sums, so they are evaluated in a straightforward way. Among the six variables in Eq. (C.10), two of them (n 1 , n 2 ) should be positive integers, but the other four (b 1 , b 2 , b 3 , c) can take any value, provided the left poles and right poles do not merge. Once the integral is expressed in terms of the Γ-function, we can easily compute derivatives of, or analytically continue the result.