Multi-loop calculations: numerical methods and applications

TL;DR

This paper reviews numerical techniques for multi-loop calculations in quantum field theory, focusing on sector decomposition methods and their application to complex two-loop integrals with multiple mass scales.

Contribution

It introduces recent advancements in sector decomposition methods and demonstrates their application to challenging two-loop integrals with multiple mass scales.

Findings

Enhanced numerical methods for multi-loop integrals.

Successful application to complex two-loop integrals.

Improved computational efficiency in multi-scale problems.

Abstract

We briefly review numerical methods for calculations beyond one loop and then describe new developments within the method of sector decomposition in more detail. We also discuss applications to two-loop integrals involving several mass scales.

Figures

Tables

Table 1. Table 1: Public tools for various steps of loop amplitude calculations beyond one loop.

\brStep to be performed	available public tools
\mrDiagram generation	Qgraf [22], FeynArts/FormCalc [29, 30]
Amplitude manipulations	Diana [31], FeynCalc [32, 33]
Reduction	Reduze [25, 26], Fire [34, 35], LiteRed [36, 37], Air [38]
Numerical evaluation	sector_decomposition [39], SecDec [40, 27], Fiesta4 [41],
of the loop integrals	Nicodemos [42], Ambre/MBnumerics [43, 44]
\br

Table 2. Table 2: Methods for the isolation of IR divergent real radiation at NNLO.

\brmethod	analytic integr. of	type/restrictions
	subtraction terms
\mrantenna subtraction [78]	yes	subtraction
$q_T$ subtraction [79]	yes	slicing; colourless final states
N-jettiness [80, 81]	yes	slicing
sector-improved residue subtraction[82, 83, 84, 85, 86, 87, 88]	no	subtraction
colourful subtraction [89, 90]	partly	subtraction; colourless initial states
\br

Equations

$$G$$

$${\mathrm{d}}^D{\kappa }_l$$

$$G$$

$$\mathcal{F}(\vec{x})$$

$$\mathcal{U}(\vec{x})$$

$$I$$

Full text

00footnotetext: ∗Speaker. To appear in the proceedings of the 4th Computational Particle Physics Workshop, October 2016, Hayama, Japan.

1 Introduction

Precision calculations are of primary importance to scrutinise the Standard Model (SM) of particle physics and in particular the Higgs sector, where experiments moved from the discovery phase to the phase of precision measurements of the Higgs properties. Small deviations from the expected values may be our only hints to physics beyond the SM for some time, and therefore precise theoretical predictions are mandatory.

In the last decade, predictions at next-to-leading order (NLO) in perturbation theory in the strong coupling constant $\alpha_{s}$ got a large boost due to advances in calculational methods, and, together with NLO matching to parton shower Monte Carlo programs, became the state of the art to describe the data. However, for the phases II and III of the LHC, and even more so at future colliders, the situation is drastically different: the experimental precision for many important SM processes already has reached a level where NLO QCD predictions fall short. Therefore, a lot of effort has been spent in the past years to come up with corrections going beyond NLO QCD, ideally not only for total cross sections, but also for differential distributions.

A measure of complexity for the calculation of higher order corrections in perturbation theory involves the number of loops in the virtual amplitude, the number of scales (Mandelstam invariants, masses) and the number of external legs. While the problem of infrared subtractions is more severe the more massless particles are involved, the difficulty to obtain analytic expressions for master integrals at two loops and beyond increases rapidly as the number of mass scales grows. Therefore numerical methods to calculate loop integrals seem particularly well suited for integrals with several mass scales.

For processes involving only massless particles, virtual two-loop 4-point amplitudes have been calculated about 15 years ago [1, 2, 3, 4, 5]. After that, the main bottleneck to be overcome to achieve NNLO predictions for processes involving massless two-loop 4-point amplitudes was the lack of an efficient subtraction scheme for the infrared singularities occurring in the real radiation part, where up to two particles can become unresolved. Filling this gap has been a very active field of research in the past years. Various methods have been devised and are still under active development; they will be listed briefly below.

The next big step in the field of NNLO QCD corrections for $2\to 2$ scattering processes was the availability of results for processes involving massive particles. Here the main problem currently resides in the virtual two-loop part of the calculation. Two problems are hampering progress here: (a) the reduction of the two-loop amplitudes to a minimal set of “master integrals” times coefficients gets increasingly complicated as the number of mass scales grows, and (b) the analytic calculation of the master integrals is extremely difficult, entering unexplored territory in terms of mathematical functions to express the occurring parameter integrals. Only very recently, analytic representations of two-loop integrals and amplitudes for $2\to 2$ scattering processes involving massive particles became available, see e.g. [6, 7, 8, 9, 10, 11, 12, 13, 14].

2 Methods and tools for two-loop calculations and beyond

The steps to perform for the calculation of a (multi-)loop amplitude can be roughly divided into four stages: (1) generation of algebraic expressions for the amplitude, (2) reduction of the amplitude to a set of “master integrals” times coefficients, (3) isolation of the ultraviolet and infrared poles and (4) evaluation of the master integrals and combination with the coefficients to obtain the amplitude. To calculate a full cross section, loop amplitudes and real radiation contributions need to be combined, which requires a suitable scheme for the isolation of infrared-divergent real radiation, which is highly non-trivial beyond one loop. It also requires the construction of a stable and fast Monte Carlo program to perform the phase space integration.

It should be mentioned that step (2) above is not mandatory. Reducing the set of integrals to a minimal “basis set” is usually beneficial, to reduce the number of integrals to calculate and to avoid large cancellations between linearly dependent integrals. However, it is also possible to evaluate the occurring integrals without reduction in a numerical approach, see e.g. [15, 16, 17]. Further, there are methods which aim to avoid the problems with IR singularities related to the split into real and virtual contributions by not performing such a partition at all [18, 19, 20, 21].

In Table 1 we give a list of some publicly available multi-purpose tools which have been developed to perform the specific tasks described above, focusing on the numerical evaluation of the loop integrals. Certainly this list is incomplete and omits a multitude of codes which may be more efficient, but are tailored to more specific classes of integrals or amplitudes. Efforts towards the development of a package that can provide all the steps listed in Table 1 by combining Qgraf [22], FORM [23, 24], Reduze [25, 26] and pySecDec [27] are described in [28].

Table 2 shows some of the subtraction schemes for infrared divergent real radiation at NNLO.

Concerning the reduction, we only listed the publicly available tools which, based on the integration-by-parts (IBP) method [45], can be used within a completely automated setup. Ideas how to reduce the computational complexity of IBP algorithms can be found in [46]. A fully automated system for amplitude generation and evaluation is also given by the Grace system [47, 48, 49].

Novel reduction methods (see e.g. [50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61]), based on ideas such as integrand reduction and maximal cuts, are very promising, but have not reached the level of automation yet which is provided by the tools listed in Table 1.

Numerous methods for the numerical calculation of multi-loop integrals have been developed in addition to the ones mentioned above, we list only a few more recent ones here: direct numerical integration in momentum space [16], dispersion relations [62], use of the loop-tree duality [63], a toolbox of various dedicated numerical techniques [64, 65], numerical solution of differential equations [66, 67], numerical extrapolation method [68, 69], numerical evaluation of Mellin-Barnes integrals [70, 71, 72], private implementations of sector decomposition [73, 74, 64, 75, 76, 77].

3 Sector decomposition

Now we will describe the program SecDec [91, 92, 40, 27] in more detail. The sector decomposition algorithm is described in [93, 94], which was inspired by earlier ideas as contained in [95, 73].

Higher order calculations in perturbation theory have in common that they involve multi-dimensional integrations over some parameters: Feynman (or Schwinger) parameters in the case of (multi-)loop integrals, or parameters related to the integration of subtraction terms over a factorised phase space in the case of infrared-divergent real radiation. Usually, these calculations are performed within the framework of dimensional regularisation, and one of the challenges is to factorise the poles in the regulator $\epsilon$.

The program SecDec  [91, 92, 40, 27] is designed to perform this task in an automated way, and to integrate the coefficients of the resulting Laurent series in $\epsilon$ numerically.

The original sector decomposition algorithm described in Ref. [93] is based on an iterative procedure, which may run into an infinite recursion. It was pointed out however [39] that the structure of Feynman integrals is such that a decomposition algorithm must exist which is guaranteed to stop, as the procedure can be mapped to a known problem in convex geometry. In Ref. [96], an algorithm was presented which cannot lead to infinite recursion and is more efficient than previously employed algorithms with this property. SecDec-3 and pySecDec  contain the implementation of a decomposition strategy (called $G_{2}$ in SecDec-3 and geometric in pySecDec ), based on a modification of the method of Ref. [96], which usually outperforms the original iterative strategy (called $X$, or iterative).

3.1 Feynman parameter integrals

Multi-loop Feynman integrals can be written in a generic form. For ease of notation, we limit ourselves to scalar integrals here. Integrals with loop momenta in the numerator, or inverse propagators, only lead to an additional function in the numerator, and can be treated in the same way. We refer to [94, 40, 97] for further details.

A scalar Feynman integral $G$ in $D$ dimensions at $L$ loops with $N$ propagators, where the propagators can have arbitrary, not necessarily integer powers $\nu_{j}$, has the following representation in momentum space:

[TABLE]

where the $q_{j}$ are linear combinations of external momenta $p_{i}$ and loop momenta $k_{l}$.

Introducing Feynman parameters in Eq. (1) leads to

[TABLE]

In the expressions above, $M$ is an $L\times L$ matrix containing Feynman parameters, $Q$ is an $L$-dimensional vector, where each entry is a linear combination of external momenta and Feynman parameters, and $J$ is a scalar expression containing kinematic invariants and Feynman parameters.

${\cal U}$ is a positive semi-definite function, which vanishes at the UV subdivergences of the graph. In the region where all invariants formed from external momenta are negative (“Euclidean region”), ${\cal F}$ is also a positive semi-definite function of the Feynman parameters $x_{j}$. If some of the invariants are zero, for example if some of the external momenta are light-like, an IR divergence may appear and ${\cal F}$ vanishes for certain points in parameter space. In the Euclidean region, the necessary condition ${\cal F}=0$ for an IR divergence can only be fulfilled if some of the parameters $x_{i}$ are zero. The endpoint singularities of both UV and IR nature can be regulated by dimensional regularisation and factored out of the functions ${\cal U}$ and ${\cal F}$ using sector decomposition.

The basic concept of sector decomposition is the following: We consider a two-dimensional parameter integral which contains a singular region where both $x$ and $y$ vanish:

[TABLE]

Our aim is to factorise the singularities for $x\to 0$ and $y\to 0$. Therefore we divide the integration range into two sectors where $x$ and $y$ are ordered:

[TABLE]

Now we substitute $y=x\,t$ in sector (1) and $x=y\,t$ in sector (2) to remap the integration range to the unit square and obtain

[TABLE]

This way the singularities are factorised into monomials, while the remaining denominator goes to a constant if the integration variables approach zero. For more complicated integrands, this procedure can be iterated until a complete factorisation is achieved.

However, after the UV and IR singularities have been extracted as poles in 1/$\epsilon$, for non-Euclidean kinematics integrable singularities related to kinematic thresholds remain. These singularities imply that ${\cal F}$ is vanishing inside the integration region for some combinations of Feynman parameter values and values of the kinematic invariants. However, the integrals can be evaluated by deforming the integration contour into the complex plane [18], as explained in detail in Refs. [92, 98].

3.2 Program structure

The program consists of two main parts, an algebraic and a numerical part. The algebraic part constructs the integrand from the list of propagators or from the graph labels, performs the sector decomposition procedure to factorise the poles in the regulator $\epsilon$, the subtractions and the expansion in $\epsilon$, and prepares the contour deformation in the case of non-Euclidean kinematics. In SecDec-3, all the algebraic steps are performed in Mathematica. In the new version [27], the algebraic part has been completely restructured and implemented in python, therefore the new version is called pySecDec .

The numerical part consists of C++ functions which are integrated numerically with the Cuba library [99]. The new program pySecDec  produces C++ code using Form [23, 24], and in addition produces C++ libraries such that the finite parametric functions representing an integral after the algebraic procedure can be linked to other programs. The basic workflow is shown in Figs. 2 and 2.

3.3 Recent program developments

In addition to the new possibilities of usage, there are various new features in pySecDec  compared to SecDec-3.0:

• •

the functions can have any number of different regulators, not only the dimensional regulator $\epsilon$, needed for example in analytic regularisation within Soft-Collinear Effective Theory [100];

• •

numerators of loop integrals can be defined in terms of contracted Lorentz vectors or inverse propagators or a combination of both;

• •

the distinction between “general functions” and “loop integrands” is removed in the sense that all features which are not loop-integral-specific are also available for general polynomial functions;

• •

the inclusion of “user-defined” functions which do not enter the decomposition has been facilitated and extended;

• •

the treatment of poles which are higher than logarithmic has been improved;

• •

a procedure has been implemented to detect and remap spurious singularities which cannot be cured by contour deformation;

• •

a symmetry finder has been added which can detect possible isomorphisms between sectors.

Version 1 of pySecDec [27] is available at http://secdec.hepforge.org/.

3.4 Phenomenological application

The numerical approach based on SecDec  has been applied to calculate massive two-loop integrals entering $gg\to HH$ at NLO, retaining the full top quark mass dependence [17, 101, 102]. The calculation is based on the setup described in Refs. [28, 103, 101]. The amplitude generation leads to about 10000 integrals before any symmetries are taken into account, which have been reduced to ${\cal O}(300)$ integrals using Reduze[25, 26]. A complete reduction could not be obtained for the non-planar 4-point integrals. The inverse propagators appearing in unreduced integrals were rewritten in terms of scalar products and directly computed with SecDec.

For the total cross section at $\sqrt{s}=14$ TeV, we found a reduction of about 14% when including the full top quark mass dependence as compared to the Born-improved HEFT approximation, where in the latter the NLO corrections are calculated in the $m_{t}\to\infty$ limit and rescaled with the full Born level result. Fig. 3 shows results for the Higgs boson pair invariant mass distribution. For further details we refer to [17, 101, 102].

4 Conclusions

We have given a brief overview on numerical methods to calculate integrals (and cross sections) beyond one-loop order, before focusing on the program SecDec, in particular the new version pySecDec . We pointed to its application within a context that goes beyond the calculation of individual master integrals, for example the possibility to use it as a library to evaluate two-loop amplitudes where the analytic expressions for the master integrals are not known.

Acknowledgments

G.H. would like to thank the organizers of CPP2016 for the great workshop, and dedicate these proceedings to the memory of Shimizu-Sensei. We also would like to thank Nicolas Greiner and Tom Zirke for collaboration on parts of the projects presented here. This research was supported in part by the Research Executive Agency (REA) of the European Union under the Grant Agreement PITN-GA2012316704 (HiggsTools) and the ERC Advanced Grant MC@NNLO (340983). S. Borowka gratefully acknowledges financial support by the ERC Starting Grant “MathAm” (39568).