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

TL;DR

This paper develops a systematic method for high-energy expansion of two-loop massive four-point diagrams, addressing non-planar integral challenges and providing explicit template integrals for Higgs pair production calculations.

Contribution

It introduces a systematic approach using analytic continuation and template integrals to evaluate two-loop integrals in high-energy limits, applicable to Higgs production and similar processes.

Findings

Resolved indefinite sign issue of Symanzik polynomial via analytic continuation.

Provided explicit template integrals for Higgs pair production.

Developed techniques for Mellin-Barnes integral solutions.

Abstract

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.

Equations

$$\log (z+i0)\stackrel{\mathrm{AC}}{=}\log (-z-i0)+i\pi ,$$

$$\frac{1}{1-m/M}=\sum _{n=0}^{\infty }{\left(\frac{m}{M}\right)}^n$$

$$X+\mathcal{O}\left((m_H^2{)}^{n_H},(m_t^2{)}^{n_t},{ϵ}^n\right)\equiv X+\mathcal{O}\left((m_H^2{)}^{n_H}\right)+\mathcal{O}\left((m_t^2{)}^{n_t}\right)+\mathcal{O}\left({ϵ}^n\right).$$

$$\int {\mathfrak{D}}^n{\alpha }^{\delta }\equiv \prod _{i=1}^n\left({\int }_0^{\infty }\frac{\mathrm{d}{\alpha }_i{\alpha }_i^{{\delta }_j}}{\Gamma (1+{\delta }_j)}\right).$$

$$\underset{ϵ,{\delta }_n,...,{\delta }_2,{\delta }_1\to 0}{\lim }X\equiv \underset{ϵ\to 0}{\lim }\underset{{\delta }_n\to 0}{\lim }\cdots \underset{{\delta }_2\to 0}{\lim }\underset{{\delta }_1\to 0}{\lim }X.$$

$$\int {\mathfrak{D}}^2{\alpha }^{\delta }\left({\alpha }_1^2+{\alpha }_1{\alpha }_2\right)={\int {\mathfrak{D}}^2{\alpha }^{\delta }|}_{{\delta }_1\to {\delta }_1+2}+{\int {\mathfrak{D}}^2{\alpha }^{\delta }|}_{{\delta }_1\to {\delta }_1+1,{\delta }_2\to {\delta }_2+1}$$

$${\alpha }_{i_1...i_n}\equiv {\alpha }_{i_1}+\cdots +{\alpha }_{i_n},{\delta }_{i_1...i_n\overline{i_{n+1}}...i_{n^{\prime }}}\equiv {\delta }_{i_1}+\cdots +{\delta }_{i_n}-{\delta }_{i_{n+1}}+\cdots +{\delta }_{i_{n^{\prime }}}.$$

$$\Gamma \left[x_1,\dots ,x_n\right]\equiv \prod _{i=1}^n\Gamma (x_i).$$

$${\alpha }_1\sim m,{\alpha }_2\sim M,{\alpha }_3\sim m,{\alpha }_4\sim M,$$

$$m\to \chi m,M\to M,{\alpha }_1\to \chi {\alpha }_1,{\alpha }_2\to {\alpha }_2,{\alpha }_3\to \chi {\alpha }_3,{\alpha }_4\to {\alpha }_4,$$

$$({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)\stackrel{\chi }{\sim }(1,0,1,0)$$

$$S=-s=-(q_1+q_2{)}^2,T=-t=-(q_1+q_3{)}^2,U=-u=-(q_2+q_3{)}^2.$$

$$m_H^2<m_t^2\ll S,T,U.$$

$$q_1^2=0,q_2^2=0,q_3^2=m_H^2,q_4^2=m_H^2.$$

$$m_H^2\sim \chi ,m_t^2\sim 1,S\sim 1,T\sim 1,U\sim 1,$$

$$I(S,T,U,m_t^2,m_H^2)=\sum _{n_H}(m_H^2{)}^{n_H}{c}_{n_H}(S,T,U,m_t^2).$$

$$q_1^2=0,q_2^2=0,q_3^2=0,q_4^2=0.$$

$$m_t^2\sim \chi ,S\sim 1,T\sim 1,U\sim 1,$$

$$I(S,T,U,m_t^2,m_H^2)=\sum _{n_H}(m_H^2{)}^{n_H}\sum _{n_t}(m_t^2{)}^{n_t}{c}_{n_H,n_t}(S,T,U,\log (m_t^2)).$$

$$S\stackrel{\mathrm{AC}}{=}{e}^{-i\pi +i0}(s+i0),T\stackrel{\mathrm{AC}}{=}sv,U\stackrel{\mathrm{AC}}{=}s(1-v),$$

$$h_0=H(0;v),h_1=H(1;v),h_2=H(2;v),h_{2,1}=H(2,1;v),$$

$$H\left(3;-\frac{u}{s}\right)=H\left(3;1-v\right)\stackrel{\mathrm{AC}}{=}-h_{2,1}+h_2{h}_1-\frac{1}{2}{h}_1^2{h}_0-\frac{{\pi }^2}{6}{h}_1+{\zeta }_3.$$

$$J_{a_1,a_2,a_3,a_4}=\int \frac{d^d\ell }{i{\pi }^{d/2}}\frac{1}{{\left(m_t^2-{\ell }^2\right)}^{a_1}{\left(m_t^2-(\ell +q_1{)}^2\right)}^{a_2}{\left(m_t^2-(\ell +q_1+q_2{)}^2\right)}^{a_3}{\left(m_t^2-(\ell -q_3{)}^2\right)}^{a_4}},$$

$$J_{1,1,1,1}\stackrel{\mathrm{AC}}{=}I$$

$$\mathcal{U}={\alpha }_{1234},\mathcal{F}=m_t^2{\alpha }_{1234}\mathcal{U}+S{\alpha }_1{\alpha }_3+T{\alpha }_2{\alpha }_4-m_H^2{\alpha }_{13}{\alpha }_4.$$

$$({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)\stackrel{\chi }{\sim }(0,0,0,0).$$

$$I$$

$$p_T^2=\frac{ut-m_H^4}{s},$$

$$t=-\frac{1}{2}\left[s-2m_H^2-\sqrt{s(s-4m_H^2-4p_T^2)}\right].$$

$$f_1(t,m_H^2)=f_2(p_T^2,m_H^2)$$

Full text

††preprint: TTP18-041

1 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 multi-scale processes. The prime examples are the Higgs+jet, Higgs pair, and Higgs+$Z$ production cross sections at the LHC. These $2\to 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 $\to$ 3 partons process at two-loops Bonciani:2016qxi . Another milestone is the numeric evaluation of the Higgs+jet Jones:2018hbb and Higgs pair production Borowka:2016ehy ; Borowka:2016ypz cross sections at two-loop level using the program SecDec Borowka:2015mxa ; Borowka:2017idc . An independent numerical evaluation of the Higgs pair production cross section is given recently Baglio:2018lrj . It is worth mentioning some recent analytic calculations of three-scale four-point two-loop diagrams; the non-planar master integrals for $\mu e$ scattering  Mastrolia:2017pfy ; DiVita:2018nnh , the planar double box integral relevant to top pair production Adams:2018kez and the planar master integrals relevant to di-photon and di-jet production Becchetti:2017abb . 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 Xu:2018eos , 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. Smirnov:2002pj . In this direction, the large-$m_{t}$ expansion of the Higgs+jet Boughezal:2013uia ; Chen:2014gva ; Boughezal:2015dra ; Boughezal:2015aha ; Caola:2015wna ; Chen:2016zka ; Neumann:2016dny and Higgs pair production deFlorian:2016uhr ; deFlorian:2013uza ; Grigo:2013rya ; deFlorian:2013jea ; Maltoni:2014eza ; Grigo:2014jma ; Grigo:2015dia ; Degrassi:2016vss 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}\ll m_{H}$ Melnikov:2016qoc ; Melnikov:2017pgf and in the small top quark and Higgs masses $m_{t}>m_{H}$ Kudashkin:2017skd are performed. For the Higgs pair production cross section, the expansion in small $m_{t}$ for the planar master integrals Davies:2018ood and in small Higgs transverse momentum Bonciani:2018omm are performed. The rest of the master integrals of Higgs pair production in the small-$m_{t}$ expansion are obtained in Ref. Davies:2018qvx together with the results of this paper. Many of the non-planar master integrals are the same as, or related to those of Ref. Kudashkin:2017skd 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. Kudashkin:2017skd 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 Beneke:1997zp . 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 Smirnov:2012gma

*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” Smirnov:1990rz ; Smirnov:2002pj . The procedure of expansion by subgraphs is implemented and can be performed in an automatic way Harlander:1997zb ; Seidensticker:1999bb . The large-$m_{t}$ expansion mentioned above belongs to this category. A new proof of the method of regions was proposed by Jantzen Jantzen:2011nz 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 non-planar massive diagrams, respectively.

2 General Idea of the Method of Regions

The procedure of the method of regions is the following  Beneke:1997zp ; Smirnov:2002pj ; Jantzen:2011nz ; Smirnov:2012gma ; Semenova:2018cwy :

• •

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 $\int_{0}^{\infty}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 Pak:2010pt ; Jantzen:2012mw . One can use the algorithm, implemented in the Mathematica package asy2.1.m Jantzen:2012mw . 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 Ananthanarayan:2018tog . 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 Smirnov:1999bza . Recently, it was proposed in Ref. Semenova:2018cwy to use yet another parametric representation, Lee-Pomeransky representation Lee:2013hzt , to apply the method of regions. For all the representations mentioned above, one has to follow Step 1 to 5 for practical calculation.

3 Notation and Technical Tools

3.1 Conventions

We distinguish the exact equal sign $``="$ and the equal sign under a certain analytic continuation. For this purpose, we introduce a sign $``\overset{\mathrm{AC}}{=}"$ and use it as, e.g.,

[TABLE]

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 $(-\pi,\pi]$. Both the left-hand side and the right-hand side of Eq. (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. (1). The equality of a series expansion like

[TABLE]

is in principle also regarded as an analytic continuation. However, when a hierarchy like $m\ll M$ is explicitly stated in the text, we use normal equal sign.

We use a simplified expression of the Landau $\mathcal{O}$ notation for more than one variable as

[TABLE]

The Euler–Mascheroni constant is denoted as $\gamma_{E}$.

We use the alpha representation to calculate Feynman integrals. The integration measure and the analytic regularization parameters are defined as

[TABLE]

The analytic regularization parameters $\delta_{j}$ play one essential role and three secondary roles:

1. (i).

The essential role is to regularize the contribution of individual regions which are divergent if we naively expand in $\alpha_{i}$. This means that individual contributions are regulator dependent, and the dependence on $\delta_{j}$ cancel after we sum all the contributions and take the limits $\delta_{j}\to 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

[TABLE] 2. (ii).

We use $\delta_{j}$ to regularize the Mellin-Barnes integral. [See the text below Eq. (142).] 3. (iii).

By shifting $\delta_{j}\to\delta_{j}+1$, we can express polynomials of $\alpha_{i}$ in the integrand. For example, when $n=2$,

[TABLE]

This property is usually used to express the integrals with higher powers of propagators. 4. (iv).

We use the property of Eq. (6) to express the higher order terms. [See the text below Eq. (42).]

The sum of the variables will be expressed as

[TABLE]

The bar on an index indicates that the variable corresponding to the index is subtracted instead of added. Sometimes $\epsilon$ and $\delta_{j}$ are treated in the same way, and in those cases we express $\epsilon$ as $\delta_{0}$. For example, $\delta_{001\bar{2}}=2\epsilon+\delta_{1}-\delta_{2}$.

Also, we introduce the following compact notation for the product of $\Gamma$-functions

[TABLE]

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 $\chi$. For example, assume that we have four variables $\alpha_{1},...,\alpha_{4}$ whose scalings are

[TABLE]

where $m\sim\chi$ is a soft parameter and $M\sim 1$ is a hard parameter. In this case, we apply a substitution

[TABLE]

to the integrand and expand in $\chi$. After that, we can set $\chi=1$. In this paper we denote the scalings (9) as

[TABLE]

or simply $(1,0,1,0)$.

3.2 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

[TABLE]

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-to-leading 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

[TABLE]

The high energy expansion in this case is two-fold. First, we treat $m_{H}^{2}$ as the soft parameter and $m_{t}^{2},S,T,U$ as the hard parameters. Afterwards, we treat $m_{t}^{2}$ 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

[TABLE]

In terms of $\chi$, the scaling we impose here is

[TABLE]

and the resulting series expression of an integral $I$ is expressed symbolically as111 In general, the coefficients $c_{n_{H}}$ depend on $\log(m_{H}^{2})$, but in our case not.

[TABLE]

A note of caution should be made regarding the dependence on $S,T,U$. Since the physical Mandelstam variables satisfy the relation $s+t+u=2m_{H}^{2}$, a similar relation should also hold for the positive Mandelstam variables. As a result, functions expressed in terms of $S,T,U$ are not unique. However this is not a problem, because at the end of the calculation, we express the result in a unique way. [See the text below Eq. (20).] We use the linear dependence of $S,T,U$ to make the second Symanzik polynomial positive definite. [See Subsection 6.1.] After the $m_{H}$-expansion, the external legs becomes massless legs.

In the second expansion, i.e. the $m_{t}$-expansion, the on-shell condition of the external momenta becomes

[TABLE]

In terms of $\chi$, the scaling we impose here is

[TABLE]

and the resulting series for an integral $I$ is expressed symbolically as

[TABLE]

The result should be expressed in a way suitable for the evaluation with the physical kinematics. In order to achieve that, the analytic continuation

[TABLE]

is applied, where $i0$ represents an infinitesimal positive imaginary number, (note that the massless on-shell condition (17) is adopted), and $0\leq v\leq 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) Remiddi:1999ew and we introduce an abbreviation for HPL as

[TABLE]

and so on. The argument of the HPL is always chosen to be $v$. For example,

[TABLE]

We use the Mathematica package HPL.m Maitre:2005uu ; Maitre:2007kp 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 replacements “crossing relations”. This subject is already well-established in the case of HPL Anastasiou:2000mf , 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.]

4 A First Example: One-Loop Box Diagram

We consider the massive one-loop Feynman integral family

[TABLE]

where the infinitesimal negative imaginary part of each denominator is implicit. The external momenta $q_{i}$ satisfy the on-shell conditions (14). We consider the box diagram shown in Fig. 3, $J_{1,1,1,1}$ and its alpha representation is

[TABLE]

where the first Symanzik polynomial $\mathcal{U}$ and the second Symanzik polynomial $\mathcal{F}$ are given by

[TABLE]

We make clear the analytic continuation in Eq. (24) because the right hand side is regularized by $\delta_{j}$ whereas Eq. (23) is explicitly $\delta_{j}$-independent. Also, we assume $m_{H}^{2}<0$ in order to ensure the convergence of the integral and perform the analytic continuation of the result to $m_{H}^{2}>0$ at the end, which turns out to be trivial.

4.1 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

[TABLE]

The expansion corresponding to this scaling is actually just the Taylor expansion of the integrand in $m_{H}^{2}$, and the original integral can be written as

[TABLE]

In particular, the leading term is identical to the box diagram with completely massless external legs.

The fact that the expansion in $m_{H}^{2}$ and the integration commute is reasonable because $m_{H}$ in the denominator of the integrand (appearing through $q_{3}^{2}=m_{H}^{2}$) is always accompanied by $m_{t}$ which regulates the integral and the limit $m_{H}\to 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

[TABLE]

or equivalently222 When we solve Eq. (28) in $t$ using the relation $s+t+u=2m_{H}^{2}$, there are two solutions. The other solution has “+” in front of the square root in Eq. (29), and it corresponds to $u$ in this case. Note that the amplitude is symmetric under $t\leftrightarrow u$, so we could choose the sign the other way around. We choose the sign such that $t=-p_{T}^{2}+\mathcal{O}(p_{T}^{2}/s,m_{H}^{2}/s)$.

[TABLE]

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

[TABLE]

which are related by Eq. (29). 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_{T}^{2},m_{H}^{2})$:

[TABLE]

On the other hand, the kinematic variable appearing in the Feynman integral is $t$ and the natural representation is $f_{1}(t,m_{H}^{2})$. Thus, we express the ingredients of Eq. (31) in terms of $f_{1}(t,m_{H}^{2})$ as

[TABLE]

where $t_{0}=t|_{m_{H}=0}$. Apparently, Eq. (31) becomes complicated when Eqs. (32), (33), (34) are substituted. However, taking into account the $m_{H}$-expansion of $f_{1}(t,0)$

[TABLE]

the expression becomes simpler

[TABLE]

where the $m_{H}$-dependent $t$ (29) is used. In general, for a given order of $m_{H}^{2}$, the difference between the strict $m_{H}$-expansion of $f_{2}(p_{T}^{2},m_{H}^{2})$ for fixed $p_{T}$ and the $m_{H}$-expansion of $f_{1}(t,m_{H}^{2})$ for the $m_{H}$-dependent $t$ is higher order in $m_{H}^{2}$.

The conclusions of this subsection are the following:

1. (i).

As the result of the $m_{H}$-expansion, the integrals reduce to integrals with massless external legs. 2. (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$.

4.2 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 (18). The integral of interest is

[TABLE]

with

[TABLE]

Here we use the positive Mandelstam variables, $S,T$, to make all the terms in $\mathcal{F}$ positive. Otherwise, hard terms could cancel and result in a soft term, which breaks the method of regions. The use of positive Mandelstam variables in Eq. (38) is conceptually not new, since it corresponds to the integral in the $u$-channel where $s<0,t<0$ and $u=-s-t>0$. The absence of a negative term in the $u$-channel is reasonable because there is no physical cut in those kinematics.

By using the package asy2.1.m Jantzen:2012mw , we reveal five relevant scalings:333 More precisely, asy2.1.m reveals the scalings which lead to homogeneous and non-scaleless integrals.

[TABLE]

The scalings of regions 2 to 5 reflect the symmetries of the integral, $\alpha_{1}\leftrightarrow\alpha_{3}$ and $\alpha_{2}\leftrightarrow\alpha_{4}$. Eq. (37) is thus expressed as the sum of the contributions from these five regions: $\sum_{i=1}^{5}I^{(i)}$.

Region 1 (all-hard region)

The region where all the alpha variables scale as $\chi^{0}$, i.e., $(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})\overset{\chi}{\sim}(0,0,0,0)$ is special, and we call this region the “all-hard region”. We can make several general statements about this region within our high-energy expansion:

1. (a).

Every integral has one all-hard region. 2. (b).

There is only one soft parameter in the all-hard region, which is $m_{t}$. 3. (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. 4. (d).

The leading order contribution of the all-hard region is $\mathcal{O}(\chi^{0})$. 5. (e).

The contribution from the all-hard region has no singularities in $\delta_{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. (37) in terms of $\chi m_{t}^{2}$, we obtain the contribution of this region as

[TABLE]

where $\mathcal{U}^{(1)}=\alpha_{1234}$ and $\mathcal{F}^{(1)}=S\alpha_{1}\alpha_{3}+T\alpha_{2}\alpha_{4}$. As stated above, the leading order term is the massless box integral and we name it $\mathcal{T}^{(1)}_{\delta_{1},\delta_{2},\delta_{3},\delta_{4},\epsilon}$ for later use:

[TABLE]

The integrand of the higher order corrections has additional factors of $\alpha_{i}$ which can be expressed by some shifts of $\delta_{j}\to\delta_{j}+1$. Indeed, Eq. (40) can be expressed as

[TABLE]

where $\mathcal{P}_{x}^{n}=\Gamma(x+n)/\Gamma(x)$ is the Pochhammer symbol. We call $\mathcal{T}^{(1)}_{\delta_{1},\delta_{2},\delta_{3},\delta_{4},\epsilon}$ a “template integral” since all the higher order terms can be expressed in terms of $\mathcal{T}^{(1)}_{\delta_{1},\delta_{2},\delta_{3},\delta_{4},\epsilon}$ by shifting $\delta_{j}$.

Since there is no singularity in $\delta_{j}$, we can safely set $\delta_{j}=0$ in Eq. (43). Then, the leading order term is expressed as

[TABLE]

where the technique explained in Appendix A is used to set the integration contour to a straight line. The expression at $\epsilon\to 0$ is obtained by using the package MB.m Czakon:2005rk , and the result is

[TABLE]

The higher order terms are given by $\mathcal{T}^{(1)}_{1,0,0,0,\epsilon}$, $\mathcal{T}^{(1)}_{0,1,0,0,\epsilon}$ 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 (23) and expand each propagator in $m_{t}$ as

[TABLE]

Then, the contribution can be expressed in terms of the integral family

[TABLE]

as

[TABLE]

This corresponds to Eq. (43). Applying the IBP-reduction, all the integrals appearing Eq. (48) including higher order terms can be expressed by three master integrals

[TABLE]

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. (39) and expanding in $\chi m_{t}^{2},\chi\alpha_{3},\chi\alpha_{4},$

[TABLE]

where $\mathcal{U}^{(2)}=\alpha_{12}$ and $\mathcal{F}^{(2)}=m_{t}^{2}\alpha_{12}\mathcal{U}^{(2)}+S\alpha_{1}\alpha_{3}+T\alpha_{2}\alpha_{4}$. The integration over $\alpha_{1},\dots,\alpha_{4}$ can be performed using the relation (143) and a variant of it, and the template integral of this region is

[TABLE]

The higher order terms in Eq. (50) which contain the inverse of $\mathcal{U}^{(2)}$ can be expressed by the shift $\epsilon\to\epsilon-1$ in the template integral. Thus Eq. (50) is written as

[TABLE]

Recall that $\mathcal{P}_{x}^{n}=\Gamma(x+n)/\Gamma(x)$ is the Pochhammer symbol.

As mentioned in Subsection 3.1, the result of the limits $\delta_{j}\to 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

[TABLE]

whereas with descending values of $j$, we instead obtain

[TABLE]

The order dependence is not problem, provided we use the same order throughout the calculation. The artifacts caused by the $\delta_{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 $\mathcal{T}^{(2)}_{\delta_{1},\delta_{2},\delta_{3},\delta_{4},\epsilon}$ with exchanged $\delta_{j}$ as

[TABLE]

and when we take the ascending order of limits, we obtain

[TABLE]

Sum of all Regions

Summing Eqs. (45), (53), (58), (59), (60), we obtain the leading term of Eq. (23):

[TABLE]

As mentioned, the result is $\delta_{j}$-independent. There are 24 possible ways to order $\delta_{1},\delta_{2},\delta_{3},\delta_{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 $\epsilon\to 0$, the poles of $\epsilon$ in the individual contributions from each region cancel.

4.3 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. (43), (52) includes the contributions up to $\mathcal{O}(\chi)$. 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 $\chi$, 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  Melnikov:2016qoc ; Kudashkin:2017skd ; Davies:2018ood ; Davies:2018qvx . We use the differential equation with respect to $m_{t}^{2}$ to obtain higher order corrections in $m_{t}^{2}$. Since we know that the integral has the form

[TABLE]

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.

4.4 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 $\alpha_{4}\to 0$ in Eq. (25), 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

[TABLE]

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 $\delta_{j}$-dependence of the alpha representation is expressed by the replacement of $a_{i}\to 1+\delta_{j}$, the triangle integral $J_{1,1,1,0}$ can be expressed by the limit $\delta_{4}\to-1$. Therefore by taking the limit $\delta_{4}\to-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,

[TABLE]

and this is the template integral for the region $(0,0,1)$. Sometimes the limit vanishes due to a suppression factor $1/\Gamma(1+\delta_{4})$. For example,

[TABLE]

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}$.

5 Two-Loop Planar Diagrams

We consider the following Feynman integral families

[TABLE]

where the momenta of the lines are

[TABLE]

Recall that $q_{13}=q_{1}+q_{3}$. We consider the integrals $I^{\mathrm{PL1}}=J^{\mathrm{PL1}}_{1,1,1,1,1,1,1}$ and $I^{\mathrm{PL2}}=J^{\mathrm{PL2}}_{1,1,1,1,1,1,1}$ whose diagrammatic representations are shown in Fig. 4.

The $m_{H}$-expansion can be performed in the same way as in Subsection 4.1, and the alpha representations of the integrals after the $m_{H}$-expansion are given by

[TABLE]

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 $\mathcal{F}$-function that it is positive definite in the $u$-channel is the same [cf. the text below Eq. (38)]. 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:

[TABLE]

The template integrals of these regions are summarized in Appendix B, and can be found in the ancillary file to this paper anci . The result of this integral at the leading order in $m_{t}$, up to $\mathcal{O}(\epsilon)$, is given by

[TABLE]

where the coefficients $d_{n_{1},n_{2}}$ are given by

[TABLE]

Single Massive Double Box Diagram

By using the package asy2.1.m, we reveal ten relevant scalings:

[TABLE]

The template integrals of these regions are summarized in Appendix B, and can be found in the ancillary file to this paper anci . The result of this integral at the leading order in $m_{t}$, up to $\mathcal{O}(\epsilon)$, is

[TABLE]

where the coefficients $d_{n_{1},n_{2}}$ are now given by

[TABLE]

6 Two-Loop Non-Planar Diagram

For the two-loop massive non-planar diagrams, we consider the following two Feynman integral families

[TABLE]

where the momenta of the lines are given by

[TABLE]

and

[TABLE]

where the momenta of the lines are

[TABLE]

Recall that $q_{2\bar{3}}=q_{2}-q_{3}$ and $q_{23}=q_{2}+q_{3}$. We consider $I^{\mathrm{NPL1}}=J^{\mathrm{NPL1}}_{1,1,1,1,1,1,1}$ as an example in this section. The template integrals of $I^{\mathrm{NPL2}}=J^{\mathrm{NPL2}}_{1,1,1,1,1,1,1}$ is provided in Appendix B. The Feynman diagrams corresponding $I^{\mathrm{NPL1}}$ and $I^{\mathrm{NPL2}}$ are illustrated in Fig. 5.

In Subsection 4.1 we showed that the expansion in $m_{H}$ can be obtained by the naive Taylor expansion of the integrand. This holds also in this case, and the $m_{H}$-expansion is straightforward. Therefore we again consider only integrals with massless external legs. The alpha representation of our non-planar integral is

[TABLE]

where

[TABLE]

6.1 Relevant Scaling

A new feature appears in Eq. (84): if we impose the relation $S+T+U=0$, we obtain

[TABLE]

and the sign of $\mathcal{F}$ becomes indefinite. In such a case, it is not guaranteed that the method to reveal the relevant scalings works properly Jantzen:2012mw . An idea to solve this problem is to perform a proper change of variables and decompose the integration domain such that $\mathcal{F}$ is positive-definite Jantzen:2012mw . However in our case, this approach does not resolve the indefinite sign of $\mathcal{F}$ since there is no simple change of variables to make $\mathcal{F}$ positive-definite.

The solution to this problem is to keep $S,T,U$ as independent variables. It is obvious that Eq. (84) 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 (84) is the following: we first compute $\mathcal{F}$ respecting the original definition of the Mandelstam variables (12). At this point, there are some redundant terms in $\mathcal{F}$ such as $(S+T+U)\alpha_{2}\alpha_{3}\alpha_{6}$. We minimize the number of terms, under the condition that $\mathcal{F}$ remains positive definite. The resulting $\mathcal{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 $\mathcal{F}$ positive definite, but in our cases this option did not solve the problem.

With this setup and the hierarchy (18), we obtain the following fourteen relevant scalings444 In fact, it turns out that the correct scalings (86) can be obtained by using Eq. (85) or by using AlphaRepExpand, provided suitable values for $S$ and $T$ are chosen (e.g. $S=1,T=1$), such that $U\not=0$. This may be because there is no cancellation between two hard-scaling terms resulting in a soft-scaling term in Eq. (85). 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 Ananthanarayan:2018tog instead of asy2.1.m and obtain the correct scalings, if one similarly chooses suitable values for $S$ and $T$.

[TABLE]

In the case of the planar integrals Davies:2018ood , the scalings consist of only 0 and 1. Here, we additionally have a scaling $(0,\frac{1}{2},\frac{1}{2},0,0,\frac{1}{2})$ which is particular to these non-planar integrals.

The contribution of the all-hard regions can be calculated by using the massless integral family and IBP-reduction [cf. the text below Eq. (46)]. Therefore we do not need a template integral for Region 1, and we show the calculation of Region 1 separately. The template integrals for Regions 2 to 14 are calculated in the next subsection, and can be found in the ancillary file to this paper anci .

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. Tausk:1999vh .

6.2 Template Integrals for Regions 2 to 14

The template integral of each region is expressed as

[TABLE]

where $2\leq j\leq 14$. For simplicity, we omit the subscripts of $\mathcal{T}^{(j)}_{\delta_{1},\delta_{2},\delta_{3},\delta_{4},\delta_{5},\delta_{6},\delta_{7},\epsilon}$ and write $\mathcal{T}^{(j)}$ when they are in the ordinary order.

Region 2

$(0,\frac{1}{2},\frac{1}{2},0,0,\frac{1}{2},1)$, Region 3 $(\frac{1}{2},0,0,\frac{1}{2},1,\frac{1}{2},0)$

The Symanzik polynomials of Region 2 are given by

[TABLE]

The integration in $\alpha_{7}$ can be performed using the relation (143). Then, we perform the following change of variables

[TABLE]

and the template integral becomes

[TABLE]

The integration in $\beta_{1},\beta_{2},\beta_{3}$ is now straightforward. To integrate over the remaining variables, it is easiest to integrate in $\alpha_{1}$ first using the relation (144) and then integrate in $\alpha_{4},\alpha_{5}$. Finally we obtain the template integral as

[TABLE]

The template integral of the Region 3 is obtained by the replacement $\alpha_{1}\leftrightarrow\alpha_{2}$, $\alpha_{3}\leftrightarrow\alpha_{4}$, $\alpha_{5}\leftrightarrow\alpha_{7}$, $T\leftrightarrow U$ of $\mathcal{T}^{(2)}$

[TABLE]

Region 4 $(0,0,0,1,1,1,0)$,

Region 5 $(0,0,1,0,0,1,1)$

The Symanzik polynomials of Region 4 are given by

[TABLE]

The integrations in $\alpha_{4},\alpha_{5},\alpha_{6}$ can be done using the relation (143). Then the template integral is

[TABLE]

We introduce a Mellin-Barnes integral to separate $\alpha_{12}\alpha_{37}+\alpha_{3}\alpha_{7}$ into two factors, $\alpha_{13}\alpha_{7}+\alpha_{2}\alpha_{37}$ and $\alpha_{1}\alpha_{3}$, then the integration in $\alpha_{1}$ and $\alpha_{2}$ can be done using the relation (144). The remaining integration is also straightforward and we obtain

[TABLE]

The template integral of Region 5 can be obtained in a similar manner and the result is

[TABLE]

Region 6 $(0,1,0,0,0,1,1)$,

Region 8 $(1,0,0,0,1,1,0)$

The Symanzik polynomials of Region 6 are given by

[TABLE]

The integrations in $\alpha_{2},\alpha_{6},\alpha_{7}$ are straightforward using the relation (143) and the template integral is given by

[TABLE]

This looks similar to Eq. (96) but has one extra massive integral, and thus it is necessary to introduce two Mellin-Barnes integrals, giving

[TABLE]

The template integral of Region 8 can be obtained in a similar manner and the result is

[TABLE]

Region 7 $(0,1,1,0,0,0,1)$,

Region 9 $(1,0,0,1,1,0,0)$

The Symanzik polynomials of Region 7 are given by

[TABLE]

The template integral is obtained in a similar way as that of Region 6 and the resulting two-dimensional integral is given by

[TABLE]

The template integral of Region 9 can also be obtained in a similar manner and the result is given by

[TABLE]

Region 10 $(1,1,0,0,0,0,1)$,

Region 11 $(1,1,0,0,1,0,0)$

The Symanzik polynomials of Region 10 are given by

[TABLE]

The template integral is obtained in a similar way as that of Region 5 and the resulting one-dimensional integral is given by

[TABLE]

The template integral of Region 11 can be obtained in a similar manner and the result is given by

[TABLE]

Region 12 $(1,1,0,0,1,1,1)$

The Symanzik polynomials of Region 12 are given by

[TABLE]

We introduce Mellin-Barnes integrals four times, and the template integral is given by

[TABLE]

Region 13 $(1,1,1,1,0,0,1)$, Region 14 $(1,1,1,1,1,0,0)$

The Symanzik polynomials of Region 13 are given by

[TABLE]

It is necessary to introduce the Mellin-Barnes integral twice in order to separate the terms proportional to $S,T$, or $U$, respectively, and we obtain

[TABLE]

The template integral of Region 14 can be obtained in a similar manner and the result is

[TABLE]

In the limit where all the regularization parameters go to zero, Eq. (117) becomes

[TABLE]

which means that the two-dimensional integral does not contribute at $\epsilon^{0}$-order. Thus, the representation (117) is most useful when the required order is $\epsilon^{0}$.

There is another representation of $\mathcal{T}^{(13)}$ which is more suitable if we require calculation beyond order $\epsilon^{0}$:

[TABLE]

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 $\delta_{3}$ and $\delta_{4}$ are non-negative integers, and Eq. (120) is shown in the special case that $\delta_{3}=0,\delta_{4}=0$, since that is the typical situation (here we may set these values because $\mathcal{T}^{(13)}$ is regular in $\delta_{3}$ and $\delta_{4}$). Another disadvantage is that each integral produces singularities in $\epsilon$ which cancel in their sum, however as a by-product, higher-order derivatives of $\Gamma$-functions contribute to the $\epsilon^{0}$ order.

6.3 Analytic Continuation

As mentioned in the text below Eq. (141), the Mellin-Barnes integrals in the template integrals are assumed to be regularized by choosing suitable values of $\delta_{j}$. However, the quantity we need is the one where $\delta_{j}\to 0$ for all $j$. Therefore we need to analytically continue the Mellin-Barnes integrals in terms of $\delta_{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.555 In terms of the dimension of the Mellin-Barnes integrals, the most involved is Region 12. However, the analytic continuation turns out to be rather simple in this region.

We take the limit of the ascending order of $\delta_{j}$,

[TABLE]

As mentioned below Eq. (141), 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 $\mathrm{Re}(z_{1})=-1/5$, $\mathrm{Re}(z_{2})=-1/3$. Then, we may choose

[TABLE]

to regularize the template integral $\mathcal{T}^{(7)}$.666 These values can be changed, provided they do not cross the integration contours.

We try to set as many parameters as possible to zero in ascending $j$ order in Eq. (122). The limit of $\delta_{1,2}\to 0$ in Eq. (121) is now trivial.

The analytic continuation of $\delta_{3}$ from $-13/30$ to 0 makes the first rightmost left poles of $\Gamma(z_{1}+\delta_{\bar{3}6}-1)$, $\Gamma(z_{2}-\delta_{3})$, and $\Gamma(z_{12}-\delta_{03})$ into right poles, which must be compensated by adding their residues. This analytic continuation procedure is automatized in the Mathematica package MB.m Czakon:2005rk . After the analytic continuation in terms of $\delta_{3}$, the integral depends on $\delta_{4},...,\delta_{7}$ and $\epsilon$, and we repeat the same procedure for $\delta_{4}$, then, $\delta_{5}$ and so on. In this way, we obtain a combination of integrals for which the arguments of the $\Gamma$-functions in the integrand contain only $z_{1}$ and $z_{2}$ such as

[TABLE]

The methods to solve these integrals are explained in the next subsection.

6.4 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 hepforge 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. Davies:2018ood , and we fit or extend them to our integrals.

Three- and Four-Dimensional Mellin-Barnes Integrals

The template integral of the Region 12 (114) is expressed as a four-dimensional Mellin-Barnes integral. Thus, the contribution from Region 12 contain a four-dimensional integral of the form

[TABLE]

We use the relation

[TABLE]

to reduce the number of the $\Gamma$-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. (125) can be done only after the limits $\delta_{j}\to 0$ and $\epsilon\to 0$ since the $\Gamma$-functions of the denominator and the numerator have different dependence on $\epsilon$, 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

[TABLE]

or integrals with a single argument of the form

[TABLE]

we first reduce them to a one-dimensional integral using the generalized Barnes lemma Davies:2018ood

[TABLE]

where ${}_{3}F_{2}$ is the generalized hypergeometric function  hypgeo1 ; hypgeo2 . A useful corollary of Eq. (130) 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 $\epsilon$ so we do not need to consider them.

One-Dimensional Mellin-Barnes Integrals

For integrals with no argument such as 777 After the analytic continuation described in Subsection 6.3, we may have an expression where some of the poles merge. The following procedure can be applied also in these cases.

[TABLE]

we evaluate them numerically using Mathematica and apply the PSLQ algorithm PSLQ1 ; PSLQ2 to fit them to a basis of constants which consists of all possible products of

[TABLE]

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

[TABLE]

Various combinations of $X_{1},X_{2}$ appear since the template integrals contain them. We obtain the series expansion of Eq. (133) 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$: 888 Below Eq. (2), it was stated that we use the normal equal sign for series representations when the hierarchy is obvious. However, here we have to introduce additional assumptions for $X=T/S,U/S,T/U,$ since hierarchies between the positive Mandelstam variables have not been fixed up to this point. Thus we use the sign “$\overset{\mathrm{AC}}{=}$” in Eq. (134), indicating that a certain analytic continuation should be performed in order to ensure $X\ll 1$.

[TABLE]

Now we apply analytic continuation and obtain

[TABLE]

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 (135), the series in $v$ is directly fit to $h_{N}$. In the case of (136), 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 (137), 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\leq v\leq 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 (78), there are about 50 one-dimensional Mellin-Barnes integrals which are treated in this way.

6.5 Combining the Results

Summing the contributions from all the relevant regions, we obtain for our sample integral

[TABLE]

where

[TABLE]

This result is consistent with Ref. Kudashkin:2017skd 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:

[TABLE]

where $\Delta(1/\delta_{j})$ has poles in $\delta_{j}$, and these poles are cancelled by the contributions from the other regions.

6.6 Other Master Integrals

Seven-Line Integrals

After the IBP-reduction described in Ref. Davies:2018qvx , we find that there are four more diagrams which have seven internal lines. We consider $J^{\mathrm{NPL1}}_{2,1,1,1,1,1,1}$, $J^{\mathrm{NPL1}}_{1,1,2,1,1,1,1}$, $J^{\mathrm{NPL1}}_{1,1,1,2,1,1,1}$, $J^{\mathrm{NPL1}}_{1,1,1,1,1,2,1}$ which are used to calculate the Higgs pair production cross section. (For the detail, see Ref. Davies:2018qvx .) These integrals can be expressed by the proper shifts of $\delta_{j}$: $\delta_{1}\to\delta_{1}+1$ for $J^{\mathrm{NPL1}}_{2,1,1,1,1,1,1}$, for example.

Six-Line Integrals

We can use the method described in Subsection 4.4 to compute the integral with fewer lines. For example, $J^{\mathrm{NPL1}}_{1,1,1,1,1,1,0}$ is obtained by shifting $\delta_{7}\to\delta_{7}-1$ and repeat the same procedure described above. Note that the analytic continuation of $\delta_{j}$ may change due to the shift. For example, the template integral of Region 12 (114) can be regularized with $\delta_{i>0}=0$, whereas $J^{\mathrm{NPL1}}_{1,1,1,1,1,1,0}$ requires a non-zero $\delta_{7}$ to regularize that template integral.

7 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 non-planar 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_{i}^{2}\ll m^{2}\ll 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}$.

Acknowledgements

The author is grateful to Joshua Davies, Matthias Steinhauser, and David Wellmann for valuable collaborative work, and Hjalte Frellesvig, Christopher Wever for useful discussions. The author sincerely thanks Matthias Steinhauser, Joshua Davies, Vladimir Smirnov, and Yukinari Sumino for carefully reading the manuscript and a lot of useful comments.

Appendix A Mellin-Barnes Integrals

Here we explain elementary aspects of the Mellin-Barnes integrals. The basic tool is the identity

[TABLE]

where the integration contour $C$ satisfies the following three properties:

• $\circ$

$C$ runs from $-i\infty$ to $+i\infty$.

• $\circ$

The poles of $\Gamma(-z)$ lie on the right side of $C$. (We call them right poles.)

• $\circ$

The poles of $\Gamma(\lambda+z)$ lie on the left side of $C$. (We call them left poles.)

The separation of the left poles and right poles is possible if $\lambda\not\in\{0,-1,-2,-3,...\}$. For example, when $\lambda=-7/2$, some of the left poles lie on the right side of the left-most right pole, but there exists a contour which separates the left poles and the right poles in the proper manner. [See Fig. 6 (a).] In our calculation, $\lambda$ is a combination of $\epsilon$ and $\delta_{j}$ and integration variables of other Mellin-Barnes integrals, and $A,B$ are combinations of the positive Mandelstam variables and the alpha parameters. In particular, $A>0$ and $B>0$ are always satisfied. We choose $\delta_{j}$ to ensure the condition $\lambda\not\in\{0,-1,-2,-3,...\}$.999 If some of the left poles and right poles merge for any choice of $\delta_{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.

When $\mathrm{Re}(\lambda)>0$, it is possible to set $C$ as a straight line along the imaginary axis [See Fig. 6 (b).]

[TABLE]

The integrand converges rapidly to zero at $\mathrm{Im}(z)=\pm\infty$ 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 $\alpha_{i}$-integration should be taken into account. For example, the relations

[TABLE]

require $\mathrm{Re}(a)>-1$ and $\mathrm{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 $\Gamma$-functions in the left hand sides of Eqs. (143), (144) is possible.

After the integration with respect to all $\alpha_{i}$ is solved, we deform $C$ to a straight line by tuning $\delta_{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 $\delta_{j}$. In order to take the limit $\delta_{j}\to 0$, it is necessary to apply a final analytic continuation of $\delta_{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 Czakon:2005rk .

Appendix B Template Integrals for the Two-Loop Master Integrals of Higgs Pair Production

The results for the template integrals listed here can be found in the ancillary file to this paper anci , in a computer-readable format.

Template Integrals of NPL2

The template integrals of $I^{\mathrm{NPL2}}$ are

[TABLE]

Template Integrals of PL1

The template integrals of $I^{\mathrm{PL1}}=J^{\mathrm{PL1}}_{1,1,1,1,1,1,1}$ defined in Eq. (66) are

[TABLE]

Template Integrals of PL2

The template integrals of $I^{\mathrm{PL2}}=J^{\mathrm{PL2}}_{1,1,1,1,1,1,1}$ defined in Eq. (67) are

[TABLE]

Appendix C A Corollary of the Generalized Barnes Lemma

The generalized Barnes lemma (130)

[TABLE]

contains the generalized hypergeometric function ${}_{3}F_{2}$. We try to reduce ${}_{3}F_{2}$ to a product of $\Gamma$-functions on a case-by-case basis, using the relations given in Refs. hypgeo1 ; hypgeo2 . In this appendix we present a formula which is, to the knowledge of the author, not published. It yields a useful corollary of the generalized Barnres lemma (130).

Consider the case where

[TABLE]

where $\mathbb{N}$ is the set of positive integers. Note that $n_{1}\not=0,n_{2}\not=0$ because otherwise the left pole and the right pole merge.

Substituting Eq. (178) into Eq. (177), the arguments of ${}_{3}F_{2}$ become

[TABLE]

and here we express this type of ${}_{3}F_{2}$ in terms of the $\Gamma$-function only. To this end, consider

[TABLE]

where $x_{1},x_{2},x_{3}$ can contain $n_{1},n_{2}$. The products are resolved by the partial fraction decomposition

[TABLE]

and we obtain

[TABLE]

We can apply the partial fraction decomposition further and obtain

[TABLE]

Here the equation in the case of $x_{1}-x_{2}\not\in\mathbb{Z}$ is shown, but the case of $x_{1}-x_{2}\in\mathbb{Z}$ is similar. It is also possible to set $x_{2}\to x_{1}+n,n\in\mathbb{Z}$ at the end. The infinite sum of $m$ in Eq. (189) is now possible, and we obtain

[TABLE]

Substituting the replacement $x_{1}\to-b_{2}+b_{3}+n_{2},x_{2}\to-b_{1}+b_{3}+n_{1},x_{3}\to b_{3}-c+n_{1}+n_{2}$ into Eq. (190), we obtain

[TABLE]

In this expression, there are only finite sums, so they are evaluated in a straightforward way. Among the six variables in Eq. (191), 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 $\Gamma$-function, we can easily compute derivatives of, or analytically continue the result.