Evaluating multi-loop Feynman integrals numerically through differential equations

TL;DR

This paper introduces a numerical method for evaluating multi-loop Feynman integrals using differential equations, improving efficiency in complex quantum field theory calculations.

Contribution

It presents a novel approach combining differential equations and sector decomposition for efficient multi-loop integral evaluation in the physical region.

Findings

Successfully computed two-loop integrals for key scattering processes

Completed master integrals for non-planar diagrams in relevant particle interactions

Demonstrated the method's effectiveness in complex multi-loop calculations

Abstract

The computation of Feynman integrals is often the bottleneck of multi-loop calculations. We propose and implement a new method to efficiently evaluate such integrals in the physical region through the numerical integration of a suitable set of differential equations, where the initial conditions are provided in the unphysical region via the sector decomposition method. We present numerical results for a set of two-loop integrals, where the non-planar ones complete the master integrals for $gg\to\gamma\gamma$ and $q\bar{q}\to\gamma\gamma$ scattering mediated by the top quark.

Figures

Tables

Table 1. Table 1 : The full list of singularities other than infinity is shown, as well as whether it is a branching point(marked as ”Y”) or not(marked as ”N”). If such point is not a singular point of corresponding family, ”-” is shown. Note that we adopt u = − s − t 𝑢 𝑠 𝑡 u=-s-t to show the crossing symmetry. For s ​ t + 4 ​ u = 0 𝑠 𝑡 4 𝑢 0 st+4u=0 it becomes a branching point only when s > 0 , t > 0 formulae-sequence 𝑠 0 𝑡 0 s>0,t>0 , thus we mark it as ”Y/N”, similarly for the other two t ​ u + 4 ​ s = 0 𝑡 𝑢 4 𝑠 0 tu+4s=0 and s ​ u + 4 ​ t = 0 𝑠 𝑢 4 𝑡 0 su+4t=0 .

	F1	F2	F3
$s=0$	N	Y	N
$s=4$	Y	Y	Y
$s=-16$	-	N	-
$t=0$	N	N	N
$t=4$	Y	Y	Y
$u=0$	N	N	N
$u=4$	-	Y	Y
$t=u$	-	N	N
$st+4u=0$	Y/N	Y/N	Y/N
$tu+4s=0$	-	Y/N	Y/N
$su+4t=0$	-	Y/N	Y/N
$4t^2-s{\left(t-1\right)}^2=0$	N	N	N
$4u^2-s{\left(u-1\right)}^2=0$	-	-	N

Table 2. Table 2 : The comparison of our numerical initial conditions obtained from Nift nift with the analytical ones for the Feynman integral I 1 subscript 𝐼 1 I_{1} and I 2 sub superscript subscript 𝐼 2 sub I_{2}^{\text{sub}} . The two initial points are: IC1( s = − 1.33 , t = − 0.891 formulae-sequence 𝑠 1.33 𝑡 0.891 s=-1.33,t=-0.891 ) and IC2( s = − 1.63 , t = − 0.632 formulae-sequence 𝑠 1.63 𝑡 0.632 s=-1.63,t=-0.632 ). c 0 subscript 𝑐 0 c_{0} is the leading term of the ϵ italic-ϵ \epsilon expansion of these finite integrals.

			$c_0$	time(s)
$I_1$	IC1	Nift	$-0.059087788\left(6\right)$	1.93
		Ref. Caron-Huot:2014lda	$-0.059087788$	–
	IC2	Nift	$-0.056016652\left(5\right)$	1.74
		Ref. Caron-Huot:2014lda	$-0.056016650$	–
$I_2^{\text{sub}}$	IC1	Nift	$0.28729542\left(1\right)$	3.55
		Ref. vonManteuffel:2017hms	$0.28729543$	–
	IC2	Nift	$0.26181028\left(1\right)$	3.57
		Ref. vonManteuffel:2017hms	$0.26181029$	–

Table 3. Table 3 : The comparison of our numerical results with the analytical ones for the Feynman integral I 1 subscript 𝐼 1 I_{1} and I 2 sub superscript subscript 𝐼 2 sub I_{2}^{\text{sub}} at the point ( s = 5 , t = − 2 ) formulae-sequence 𝑠 5 𝑡 2 (s=5,t=-2) . The IC1 and IC2 denotes the two different choices of the initial conditions. c 0 subscript 𝑐 0 c_{0} is the leading term of the ϵ italic-ϵ \epsilon expansion of these finite integrals.

$\left(s=5,t=-2\right)$		$c_0$	time(s)
$I_1$	IC1	$0.573661717-i0.45602298$	0.11
	IC2	$0.573662051-i0.45602316$	0.10
	Ref. Caron-Huot:2014lda	$0.573661756-i0.45602309$	–
$I_2^{\text{sub}}$	IC1	$-0.077764616+i0.34306744$	0.26
	IC2	$-0.077764595+i0.34306737$	0.23
	Ref. vonManteuffel:2017hms	$-0.077764620+i0.34306741$	–

Table 4. Table 4 : Comparison between numerical results obtained with our algorithm from two differential choices of initial conditions for the Feynman integral I 2 subscript 𝐼 2 I_{2} and I 3 subscript 𝐼 3 I_{3} at the point ( s = 5 , t = − 2 ) formulae-sequence 𝑠 5 𝑡 2 (s=5,t=-2) . c 0 , c 1 subscript 𝑐 0 subscript 𝑐 1 c_{0},c_{1} and c 2 subscript 𝑐 2 c_{2} denotes the first three coefficients in the Laurent series of ϵ italic-ϵ \epsilon . The results obtained from pySecDec  Borowka:2017idc is also shown for consistency check and the corresponding setup is not optimal.

$\left(s=5,t=-2\right)$		$c_0$	$c_1$	$c_2$	time(s)
$I_2$	IC1	$0.02188084-i0.00000002$	$-0.0870259+i0.05170117$	$-0.246416-i0.17602070$	0.26
	IC2	$0.02188080+i0.00000001$	$-0.0870262+i0.05170118$	$-0.246417-i0.17602072$	0.23
	pySecDec	$0.02187\left(3\right)+i0.00003\left(3\right)$	$-0.0869\left(3\right)+i0.0518\left(4\right)$	$-0.248\left(2\right)-i0.175\left(2\right)$	$𝒪\left({10}^4\right)$
$I_3$	IC1	$-0.0599222+i0.4204527$	$-1.2093294+i1.1271787$	$-3.737851+i0.435880$	0.74
	IC2	$-0.0599219+i0.4204528$	$-1.2093298+i1.1271798$	$-3.737851+i0.435879$	0.78
	pySecDec	$-0.05998\left(7\right)+i0.42048\left(8\right)$	$-1.2100\left(7\right)+i1.1262\left(7\right)$	$-3.737\left(3\right)+i0.430\left(3\right)$	$𝒪\left({10}^4\right)$

Equations

$$I={\left(\frac{e^{ϵ{\gamma }_E}}{i{\pi }^{\frac{d}{2}}}\right)}^L\int \prod _{i=1}^L{\mathrm{d}}^d{k}_i\frac{1}{{\prod }_{j=1}^N{D}_j^{a_j}},$$

$$\frac{\partial I(x;ϵ)}{\partial x_i}=J_i(x;ϵ)I(x;ϵ),i=1,\cdots ,K,$$

$$I=\sum _{i=0}^{+\infty }{c}_i{ϵ}^{p+i},$$

$$\frac{\Delta I}{I}=$$

$${\epsilon }_{\text{rel}}\left[\Delta I,I\right]=$$

$$\begin{array}{rl}& (s_a,t_a)\\ \to & (i\sqrt{-4s_a},t_a)\\ \to & (2+i\sqrt{-s_a}+i\sqrt{s_b-4},t_a)\\ \to & (2+i\sqrt{-s_a}+i\sqrt{s_b-4},(t_a+t_b)/2+0.1i)\\ \to & (2+i\sqrt{-s_a}+i\sqrt{s_b-4},t_b)\\ \to & (4+i\sqrt{4(s_b-4)},t_b)\\ \to & (s_b,t_b).\end{array}$$

$$\begin{array}{r}{D}_1=k_1^2-m^2,D_2=(k_1-p_1{)}^2-m^2,D_3=(k_1-p_1-p_2{)}^2,\\ D_4=k_2^2-m^2,D_5=(k_2-p_3{)}^2-m^2,D_6=(k_2-p_1-p_2{)}^2-m^2,\\ D_7=(k_1-k_2{)}^2.\end{array}$$

$$\begin{array}{r}{D}_1=k_1^2,D_2=(k_1-p_1{)}^2,D_3=(k_1-p_1-p_2{)}^2,\\ D_4=k_2^2-m^2,D_5=(k_2-p_1-p_2+p_3{)}^2-m^2,\\ D_6=(k_1-k_2{)}^2-m^2,D_7=(k_1-k_2-p_3{)}^2-m^2.\end{array}$$

$$\begin{array}{r}{D}_1=k_1^2,D_2=(k_1-p_1{)}^2,\\ D_3=k_2^2-m^2,D_4=(k_2-p_4{)}^2-m^2,D_5=(k_2-p_3-p_4{)}^2-m^2,\\ D_6=(k_1-k_2{)}^2-m^2,D_7=(k_1-k_2+p_2{)}^2-m^2.\end{array}$$

Full text

††preprint: CP3-18-71, MCNET-18-32

1 Introduction

With the advent of Run II of the Large Hadron Collider (LHC), a wealth of experimental measurements is expected to be performed at very high luminosities, probing the high energy scales extensively for the first time. To exploit the full potential of these experimental measurements, theoretical predictions for the scattering processes are required with an unprecedented accuracy and precision. In several cases, the foreseen experimental precision will demand the inclusion of higher order terms in the perturbative expansions of the gauge coupling constants of the standard model, de facto requiring the evaluation of multi-loop amplitudes. Even though for processes that can be mediated by heavy particles specific results have been obtained over the years Graudenz:1992pv ; Spira:1995rr ; Baernreuther:2012ws ; Czakon:2013goa ; Borowka:2016ehy ; Borowka:2016ypz ; Jones:2018hbb ; Baglio:2018lrj , a general algorithm to efficiently, analytically and automatically compute the corresponding amplitudes is still lacking and poses an enormous challenge. As of now, practical methods often rely on approximations Bolzoni:2010xr ; Cacciari:2015jma ; Anastasiou:2015ema ; Anastasiou:2016cez ; Brucherseifer:2014ama ; Berger:2016oht and/or expansionsBonciani:2018omm ; Davies:2018qvx .

In general, a multi-loop amplitude can be expressed in terms of a finite set of integrals, usually known as master integrals. Although various methods for calculating the master integrals have been proposed (see Ref. Smirnov:2012gma for a review), a fully general/universal one is not yet available (see Refs. Lee:2017qql ; Liu:2017jxz ; Borowka:2018dsa for recent developments). However, the master integrals can be shown to satisfy differential equations Kotikov:1990kg ; Remiddi:1997ny ; Gehrmann:1999as , which after the reduction to a canonical form Henn:2013pwa ; Argeri:2014qva , can be in some cases solved, iteratively. Although various results have been obtained in presence of the massive particles Caron-Huot:2014lda ; Bonciani:2016qxi ; vonManteuffel:2017myy ; Becchetti:2017abb ; Mastrolia:2017pfy ; Lee:2018rgs ; Ablinger:2018yae ; DiVita:2018nnh ; Chen:2018fwb , the final results are often represented as (iterated) integrals whose integrands consist of polylogarithms and other irrational functions, which still require numerical integration.

On the other hand, although solving differential equations numerically is a well-studied topic in applied mathematics, only a few phenomenological applications Boughezal:2007ny ; Czakon:2007qi ; Czakon:2008zk have been reported, till now. In such cases, the initial condition was obtained via expansions around a singular point, and finding out such expansions for other processes is highly non-trivial.

In the present work, we explore the possibility of evaluating Feynman integrals numerically through differential equations, where the initial conditions are provided using the sector decomposition method Binoth:2000ps . The basic idea is simple: obtain the initial conditions in the unphysical region, which is a fast and accurate procedure, and then use the differential equations to analytically continue the results into the desired physical region.

This work is organised as follows. In section 2 we describe the method in detail. In section 3 we illustrate the reach of our method by computing several two loop examples relevant for $gg\to\gamma\gamma$ and $q\bar{q}\to\gamma\gamma$ mediated by the top quark. We draw our conclusions in section 4.

2 Method

We define a (scalar) Feynman integral In $d=4-2\epsilon$ dimensions by

[TABLE]

where $L$ is the number of loops, $k_{i}$ is the loop momentum, $N$ is the number of propagators and $D_{j}=q_{j}^{2}-m_{j}^{2}+i0^{+}$ is the denominator of the $j$-th propagator, where $q_{j}$ is the linear combination of the loop momenta and the external momenta, and $m_{j}$ is the corresponding mass. The $a_{j}$ denotes the respective power of the denominator.

The modern approach of multi-loop integrals consists in dividing the integrals into different topologies depending on their propagators. For each topology, a set of integration-by-parts (IBP) identities Chetyrkin:1981qh , relating different integrals, is generated exploiting the Poincaré invariance of the integrals. With such system of linear identities at hand, any integral with the same topology can be written as a linear combination of a finite subset of integrals, called the master integrals. Using the fact that derivatives of the master integrals with respect to the external kinematic variables and internal masses yield a linear combination of Feynman integrals in the same topologies, IBP relations can be used to reduce them back to the linear combination of the master integrals, leading to a system of first order partial differential equations.

Let us consider a vector of $M$ master integrals $I=(I_{1},I_{2},\cdots,I_{M})^{\mathrm{T}}$, depending on $K$ independent kinematic variables $x=(x_{1},x_{2},\cdots,x_{K})$ and $\epsilon$, one can express the set of equations as

[TABLE]

where $J_{i}$ is an $M\times M$ matrix, whose elements are rational functions of the kinematics $x$ and the dimension $d$. Each element of $J_{i}$ contains singularities originating from both the kinematics and the dimension $d$. The singularities from the kinematics are governed by the Landau equations Landau:1959fi , while the poles on $d$ must be rational numbers.

Although formally Eq. (2) is a set of partial differential equations, only one initial condition is needed to fix the solution and as a result such system can be integrated iteratively with respect to the kinematics, thereby making them similar to ordinary differential equations. Therefore, the method for initial value problems Stoer2002 can be applied straightforwardly to obtain the solution of the differential equation of the integrals. The main challenge is to obtain the suitable initial conditions and design subsequent integration contours to fully fix the solution.

In the previous studies, an expansion around singular points Czakon:2008zk ; Lee:2017qql ; Liu:2017jxz was suggested.111In Ref. Caffo:2002ch , the results of the integral at singular points were adopted, which conflicts the Lipschitz condition and becomes ill-defined, thus requiring a modification, which is equivalent to an expansion around singular points. However, such expansion is highly non-trivial, and the short distance to singular points would lead to loss of accuracy and efficiency.222Here the loss of accuracy means the accuracy on the target points are much lower than the accuracy on the initial conditions.

As the numerical algorithms are based on or related to the Taylor series expansion, the ideal initial conditions should be at the regular points, far away from all the singularities. However, the computation of the integrals at those points by analytical or semi-analytical methods is as complicated as obtaining results at any regular points in the physical region.

In this work, we propose to obtain the initial conditions for the differential equations through the sector decomposition method Binoth:2000ps . All the ultraviolet and infrared divergences of the integrals are isolated in terms of a Laurent series in $\epsilon$, by dividing the integration domain and performing variable transformations according to well-designed strategies Bogner:2007cr ; Kaneko:2009qx . The series can be expressed in the following form

[TABLE]

where $c_{0}$ represents the leading term, and the integer $p\in\mathbb{Z}$ is determined by the strength of the divergence of the integral. The numerical values of the coefficients $c_{i}$ are obtained after performing a multi-dimensional integration. In the unphysical region, where the $i0^{+}$ prescription is no longer needed,333Here we refer it as the unphysical region, but it also includes the physical region below all the thresholds of the internal particles so that the $i0^{+}$ prescription is not needed. especially in the Euclidean region, the integrands are sufficiently flat to achieve high precision through suitable multi-dimensional integration algorithm such as quasi-Monte Carlo algorithm Li:2015foa . At this point, one can exploit the analytic properties of the Feynman integrals: considering the integral as a complex function, the differential equations themselves can provide the analytical continuation from the unphysical to the physical region. As a consequence, the results of the integral in the physical region can be obtained as a Laurent series in $\epsilon$ as expressed in Eq. (3). On the other hand, with suitable contour deformations Soper:1999xk ; Anastasiou:2007qb , the sector decomposition method can also provide the results for the physical kinematics. Such a deformation, however, requires a rather complicated variable transformation. In addition, the integrands still having large oscillations exhibit poor convergence in numerical integration. Therefore, the direct computation via sector decomposition in the physical region tends to be computationally quite heavy.

An alternative path can be followed, by choosing the initial conditions in the unphysical region first and then by carefully choosing the contour of the integration. To preserve the physical $i0^{+}$ prescription, the general idea is that along the contour except the target point, the integral do not require $i0^{+}$ prescription, and the target point is approached following the $i0^{+}$ prescription. In general, constructing such contour is highly non-trivial, and we give an example of a contour for those master integrals later in sec.3. The contour is constructed carefully after the study of the branch cuts of the integrals and we leave the automation of the choice of the contours for the future.

As argued in the Ref. Czakon:2008zk , explicit methods are sufficient to solve the system of differential equations. They can be broadly organised into three classes: one-step (Runge-Kutta methods), multi-step, and extrapolation methods. In practice, the final choice of a method in a specific problem depends on several criteria, including efficiency and availability. In this work, we focus on one-step methods, mainly due to the following reasons:

One-step methods only require one initial condition, in contrast with multi-step methods.444Some implementations of multi-step methods only apparently require one initial condition as one-step methods are used to provide other initial conditions. This offers a great advantage since providing multiple initial conditions is a problem and may enhance uncertainty. In addition, it grants more freedom on the choice of the integration contour, as a piecewise contour can be adopted naively, e.g., the contour in section 3. 2. 2.

The one-step methods are linear and with simple numerical coefficients. This yields negligible overhead time and very good numerical stability.

We find that in order to achieve optimal efficiency, it is desirable to also introduce an adaptive step-size control, as implemented in the Runge-Kutta-Fehlberg method. In this method, for each step, two estimates of the results are obtained, and the difference $\Delta I$ of them is calculated. Now, we have to define a relative error based on the $\Delta I$ and $I$, and then the adaptive step-size control is obtained through the comparison of this relative error to the desired local accuracy. Since $\Delta I$ and $I$ are Laurent series in $\epsilon$, the definition of the relative error is ambiguous. Now, for the purpose of defining a relative error, we observe the following facts. Firstly, as the integrals contribute in a non-trivial way to the evaluation of the amplitude, the uncertainty of each integral at the target point should be determined from the required precision on the value of the total amplitude itself. However, this determination can be very complicated in practical applications. Secondly, the uncertainties for the intermediate points should be based on the uncertainty of the final target point only, which is not known a priori, hence impossible to apply. Thirdly, we observe that while calculating the amplitude, the uncertainties from different orders of $\epsilon$ usually mix together. Keeping these points in mind, we introduce the notion of the relative error of the integral, a quantity which is independent of any prefactor and based on the whole master integral rather than its individual terms in the $\epsilon$ expansion. Considering a master integral $I$ with the difference $\Delta I$ described previously, we define the relative error $\varepsilon_{\text{rel}}[\Delta I,I]$ based on the ratio $\frac{\Delta I}{I}$ as following:

[TABLE]

And the maximum value of relative errors $\epsilon_{\textrm{rel}}[\Delta I,I]$ in the whole family is compared to the desired local accuracy, to control the step-size.

3 Results

In the following, we demonstrate our method with three different planar and non-planar two-loop integral families, which appear in di-photon, di-jet production mediated by the heavy quarks. The diagrams are given in Fig. 1, where $p_{1},p_{2}$ are incoming and $p_{3},p_{4}$ are outgoing. The thin lines represent the massless particles, while the thick lines represent massive particles. All external lines are on-shell $p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=p_{4}^{2}=0$, and the kinematic variables are defined as $s=(p_{1}+p_{2})^{2},t=(p_{1}-p_{3})^{2},u=(p_{1}-p_{4})^{2}$, which satisfy $s+t+u=0$. We normalise the invariants by the squared internal mass $m^{2}$, effectively setting $m^{2}=1$, and the $m^{2}$ dependence can be recovered later, by power counting.

As explained before, for each master integral, we adopt Niftnift to obtain the numerical results in the Euclidean region by the sector decomposition method, where the final numerical integration is performed with the quasi-Monte Carlo algorithm. We perform the IBP reduction with the C++ version of FIRE5 Smirnov:2014hma together with LiteRed Lee:2012cn ; Lee:2013mka , to obtain the corresponding differential equations in $s$ and $t$, treating them as independent variables. We perform the numerical integration of the differential equations with odeint odeint , and the Runge-Kutta-Fehlberg 7(8)-th order method Fehlberg:68 ; Stoer2002 is chosen, based on our experimentation on one-loop integrals.

We consider the target physical region defined by $s>4,t<0,u<0$, and choose the initial conditions lying in the region with $s<0,t<0,u<4$. We perform the evolution from the initial point $(s_{a},t_{a})$ to the target point $(s_{b},t_{b})$ along the following contour formed by six line segments:

[TABLE]

In particular, we consider the target point with $(s,t)=(5,-2)$, and we choose two different points in the Euclidean region as the initial points: one is marked as IC1, with $(s,t)=(-1.33,-0.891)$; another is marked as IC2, with $(s,t)=(-1.63,-0.632)$. The difference between the results obtained from those two different initial conditions provides an estimate of the uncertainties. We list all branch points on the physical Riemann sheet in table 1, and we verified that the above contour never crosses branch cut, as can be seen in fig. 2. Alternatively, instead of determining the branch points and the branch cuts, along the contour the sector decomposition method can be adopted to calculate the numerical values of the Feynman integrals directly, since we require that along the contour the $i0^{+}$ prescription is not needed. Such numerical values provide another cross check on the results obtained from the numerical integration of differential equations.

All the timings reported here are based on a laptop with Intel Core i5-6200U CPU and the time cost consists of the evaluation of all the master integrals in the whole family. We require the relative error on the initial conditions less than $10^{-7}$, and the relative error tolerance in each step of the differential equations is set to $10^{-10}$.

We begin with the family $F_{1}$, where the analytical results in $d=4$ dimension have been reported in Ref. Caron-Huot:2014lda . We choose the denominators as555Technically, to perform the IBP reduction, two extra denominators should be chosen. However, we choose all master integrals to be scalar master integrals without any numerator, hence the results are independent of the exact form of the auxiliary denominators in the IBP reduction. We neglect the two extra denominators here for simplicity. The above comment also applies to the other two families.:

[TABLE]

We denote the integrals in this family as $I(F_{1},a_{1}a_{2}a_{3}a_{4}a_{5}a_{6}a_{7})$, where $a_{i}$ is the corresponding propagator power, as described in Eq.  (1). Working in $d=4-2\epsilon$ dimension, after the IBP reduction, we obtain 29 master integrals. In Table 2, we show the initial conditions of one of the top level master integrals $I_{1}=I(F_{1},1111111)$. As mentioned before, the relative uncertainty on the initial conditions are required to be less than $10^{-7}$, and our results are consistent with analytical onesCaron-Huot:2014lda within such uncertainty. Using those two initial conditions, we evaluate these integrals for the benchmark value($s=5,t=-2$) in the physical region, and the results of $I_{1}$ are shown in Table 3. We also report the numerical value obtained from the analytical expression in Ref. Caron-Huot:2014lda . We find that the uncertainty of our numerical results compared to the analytical one is less than $10^{-6}$. Moreover, the difference between the results obtained using the initial conditions from IC1 and IC2 is also of the same order, providing a good estimate on the uncertainty. We note that to reach such high precision takes only 0.1s.

The next example is the family $F_{2}$, shown in Fig. 1(b), with the following denominators:

[TABLE]

There are 36 master integrals in this family, and some of them involve infrared divergences. The most complicated integrals in this family, i.e. the seven-propagator master integrals, are still unknown in literature 666Partial results has been reported in Ref. Xu:2018eos recently.. Instead, for comparison, we show numerical results for one non-planar integral in the lower sector, defined by $I_{2}^{\text{sub}}=I(F_{2},1011111)$777An alternative numerical evaluation for this topology has been reported in Ref. Bonciani:2018uvv (shown in fig. 1(d)), which has been studied in Ref. vonManteuffel:2017hms and in fact is independent of $t$. In Table 2, we show our numerical initial conditions obtained from Nift as well as the analytical one on $I_{2}^{\text{sub}}$. The uncertainties on the initial conditions are less than $10^{-7}$ and the computing time is well under control.

In Table 3, we show our numerical results as well as the analytical one on $I_{2}^{\text{sub}}$ in the physical region with $s=5$. Similarly to $I_{1}$, the uncertainty from our approach is less than $10^{-6}$. The time cost is several times larger than F1, but still less than 1 second. At the same time, we also obtain the results for the seven-propagator integral $I_{2}=I(F_{2},1111111)$. As no analytical results are known for this, we use pySecDec Borowka:2017idc to obtain the results for cross check. In table 4 we show the results for all the coefficients starting from $\epsilon^{-2}$ to $\epsilon^{0}$ in $\epsilon$ expansion. By estimating the uncertainties of our method through the difference between the two results, the relative error is at $\mathcal{O}(10^{-6})$. This is much more accurate than directly evaluating it via the sector decomposition method in the physical region.

Finally, we consider family $F_{3}$, shown in Fig. 1(c). This family 888Results in the Euclidean region has been reported recently in Ref. Xu:2018eos . contains 51 master integrals, and in particular five of them belong to the seven-propagator sector, indicating more complicated differential equations than the family $F_{1}$ and $F_{2}$. The propagators are given by:

[TABLE]

We use the same points IC1 and IC2 to obtain the initial conditions. We show the numerical results for $I_{3}=I(F_{3},1111111)$ in the Table 4 and further checked with pySecDec. The computing cost of our method is still less than one second, and the precision of our results is still at $\mathcal{O}(10^{-6})$.

The computing cost on multi-dimensional integration for obtaining the initial conditions varies from 17 seconds to 2 minutes depending on the complexity, which is much less than the time spent for IBP reduction, hence negligible in practical application. The number of steps for the numerical integration of the differential equations ranges from 61 to 133, thereby indicating that the discretisation error associated with the differential equations is at most around $10^{-8}$. As explained before, the dominant uncertainties come from the uncertainties on the initial conditions, and we verified it by adjusting the relative error tolerance on the initial conditions and/or the differential equations. Further information, including numerical results for all master integrals in the three integral families are available as ancillary files with the arXiv submission.

Remarkably, one does not need to start from the Euclidean region each time. Once the results at one physical point are obtained according to previous procedure, they can be adopted as the new initial condition, for other physical points. As the branch points and branch cuts in the physical region are well-understood, comparing to the general cases, much simpler contours can be adopted.

4 Conclusion and discussion

In this paper, we have presented a method to compute the Feynman integrals numerically. The main idea is to integrate the differential equations numerically, with the initial condition in the Euclidean region provided through the sector decomposition method. We have compared numerical results achieved by our method with the available analytical ones, for several two-loop examples, and shown $\mathcal{O}(10^{-6})$ accuracy can be reached within one second. Using the above method, we have provided numerical results of several two-loop integrals, whose analytical expressions are currently unknown. Those two-loop integrals complete the two-loop master integrals for $gg\to\gamma\gamma$ and $q\bar{q}\to\gamma\gamma$ scattering mediated by the top quark, and thus our results can be applied directly to investigate the role of top quark in di-photon production at the LHC.

The differential equations of the integral encode the full $\epsilon$ dependence, while the sector decomposition can provide any higher order terms in $\epsilon$. Therefore, the results of the integral at any order of $\epsilon$ expansion can be achieved within our method, which is usually desirable and required in practical applications.

Although in this paper we restricted it to the case with real masses only, our method can be applied to the case with complex masses. Clearly, the sector decomposition method works with complex masses. On the other hand, the integration contour still doesn’t cross any branch cut if the width is small. Such complex mass scheme, is crucial and essential to describe the threshold behaviour for processes involving unstable particles. In that case, we hope our method will provide an important role to obtain the precise prediction of relevant processes.

Our method builds up on the idea that the differential equations of the master integrals can be integrated numerically, providing an initial condition in the unphysical region and a suitable integration contour.

However, it is not always possible to obtain the differential equations of the master integrals as the IBP reduction usually fails in case of the integrals having a large number of scales. In this context, for example, the recent proposal Mastrolia:2018uzb ; Frellesvig:2019kgj of the use of intersection theory could overcome this problem. The initial conditions are obtained by using the sector decomposition method, which is quite efficient in the Euclidean region. While it can be a problem for massless cases, for processes with massive loop propagators, usually such region can be found.

Finally, we note that as both the IBP reduction and sector decomposition can be done systematically and automatically, our method could play a potential role towards an automated approach and framework to multi-loop computations. This, of course, only if an algorithm for the automatic determination of the integration contours could be identified. Work in this direction is in progress.

Acknowledgements.

We thank Fabio Maltoni and Pierpaolo Mastrolia for useful discussions and comments on manuscript. The authors would also like to thank Johannes Henn, Roman N. Lee and Andreas von Manteuffel for discussions and providing additional numerical results on cross checks. XZ has received funding from the European Union’s Horizon 2020 research and innovation programme as part of the Marie Skłodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). MKM has been supported by the Research Fellowship grant awarded by the Belgian Federal Science Policy Office (BELSPO) and by the UniPD STARS Grant 2017 “Diagrammalgebra”.