Evaluation of multi-loop multi-scale Feynman integrals for precision physics

Modern particle physics is increasingly becoming a precision science that relies on advanced theoretical predictions for the analysis and interpretation of experimental results. The planned physics program at the LHC and future colliders will require three-loop electroweak and mixed electroweak-QCD corrections to single-particle production and decay processes and two-loop electroweak corrections to pair production processes, all of which are beyond the reach of existing analytical and numerical techniques in their current form. This article presents a new semi-numerical approach based on differential equations with boundary terms specified at Euclidean kinematic points. These Euclidean boundary terms can be computed numerically with high accuracy using sector decomposition or other numerical methods. They are then mapped to the physical kinematic configuration with a series solution of the differential equation system. The method is able to deliver 8 or more digits precision, and it has a built-in mechanism for checking the accuracy of the obtained results. Its efficacy is illustrated with examples for three-loop self-energy and vertex integrals and two-loop box integrals.

I. Introduction. With the discovery of the Higgs boson at the Large Hadron Collider (LHC), all building blocks of the Standard Model (SM) have been experimentally confirmed, with the only exception of the Higgs selfcoupling, which still awaits direct measurement. However, the SM does not account for important phenomena such as dark matter and the matter-antimatter asymmetry, so that physics beyond the SM is needed. It is reasonable to expect that this new physics couples to the electroweak and/or Higgs sector of the SM, since there are important model-building constraints for couplings to the strong force [1]. Therefore, possible evidence for such new physics can be explored in precision studies of electroweak and Higgs physics at the high-luminosity run of the LHC (HL-LHC) or one of several proposed future e + e − colliders: FCC-ee [2], CEPC [3], ILC [4,5], CLIC [6,7]. Through their high integrated luminosities of several ab −1 , these machines will be sensitive to very small deviations between the measured value and the SM expectation for a given observable. Thus they can probe extremely feebly coupled new particles or very large new physics scales of tens of TeV.
The SM predictions for these precision analyses are obtained by computing higher-order quantum corrections. At the HL-LHC, some of the most interesting precision studies are Higgs boson production and lepton pair (Drell-Yan) production. For the former, one of the largest sources of theoretical uncertainty stems from mixed QCD-electroweak corrections [8,9]. While some partial results at this order have been computed [10][11][12][13][14], contributions from electroweak diagrams with internal top quarks, both for 3-loop Higgs production and 2-loop Higgs+jet production, are still needed to complete this missing piece. For Drell-Yan production, 2-loop electroweak corrections for the full process pp → + − , not just on the Z-boson resonance, are important since LHC measurements cover a broad range of invariant mass [15,16].
Similarly, electroweak 2-loop corrections for several different pair production processes will be essential for the physics goals of future e + e − colliders [17]: e + e − → W + W − , e + e − → ZH, and e + e − → ff . Measurements of these cross-sections will allow us to determine the W-boson mass with high precision, constrain anomalous couplings between gauge bosons and/or the Higgs boson, and probe heavy neutral vector bosons (Z bosons). Currently, some results for mixed QCD-electroweak 2-loop corrections are available [18][19][20][21], but so far no complete electroweak 2-loop calculation for any pair production process has been carried out. Even higher-order corrections will be needed for studies of Z-boson production and decay at these future e + e − colliders, as well as the indirect prediction of the W boson mass from the Fermi constant. To match the expected experimental precision, 3-loop and partial 4-loop self-energy and vertex corrections will be required [17,22], which is one order of perturbation theory beyond the current state of the art [23].
It should be emphasized that these are loop corrections in the full SM, involving many massive particles inside the loops. The currently most advanced techniques for analytically computing such multi-loop Feynman in-tegrals first reduce them to a small set of master integrals, which then are solved by constructing suitable differential equations (DEs), see Ref. [24] for a recent review. Both of these steps require integration-by-parts (IBP) equation systems [25,26] that become intractable for multi-loop integrals with many masses. Instead, one must resort to numerical integration techniques.
The recent calculation of full 2-loop corrections to Zboson production and decay [23,27,28] made use of numerical evaluations based on sector decomposition (SD) [29][30][31][32][33] and Mellin-Barnes (MB) representations [27,[34][35][36][37][38]. However, these methods require large amounts of computing resources and do not always converge to the required level of accuracy, so that a straightforward extension to more loops and/or legs is not possible. Due to numerical cancellations between individual loop integrals, at least 8 digits precision are required in many cases for practical applications [28]. Another interesting approach, which allows to tackle a wider class of problems, evaluates a set of master integrals by numerically solving a DE system, either in terms of kinematic parameters [39,40] or in terms of an auxiliary flow variable [41][42][43].
This article introduces an efficient but still very general approach that can be applied to many challenging 2-and 3-loop problems with multiple mass and momentum scales [44]. The key elements are a system of DEs, with boundary terms evaluated at one or more Euclidean (space-like) kinematic points (which can be reliably determined to high precision with numerical methods). The DEs are then solved, using series expansions, to obtain the final result at the physical Minkowski (time-like) kinematic point. This approach,which is already fully automated in its main parts, will be described in more detail in the next section. In section III we will apply this technique to examples of SM self-energy and vertex Feynman integrals that occur in three-loop Z-decay corrections. The chosen examples are very difficult to evaluate with other analytical or numerical methods. A summary and outlook are given in the final section. The supplemental material included with this submission contains input parameters for the example integrals, additional examples for the application of our calculation technique, e.g. a 2-loop four-scale box diagram, and miscellaneous implementational details and remarks for a more detailed discussion of the method.
II. Description of the method. Solving Feynman integrals from DEs is an approach initiated in the last decade of the last century [45][46][47][48]. Many families of Feynman integrals admit a choice of master integrals for which the system of DEs has a particularly simple "canonical" form [49], which in many cases can be straightforwardly solved in terms of multiple polylogarithms.
More generally, not all Feynman integrals are of polylogarithmic type, and it can become increasingly difficult to find a closed set of analytic functions in terms of which the DEs can be solved. In such cases, it is often still possible to find precise numerical solutions. For example, the DE system can be integrated numerically [39][40][41][42][43]. In this work, we use the approach of iterated series expansions, which is both fully automatable and numerically efficient.
Let us give a brief overview of the method. Consider a basis of master integrals (MIs), F (x, ), depending on a single scale x. We work in dimensional regularization, with D = 4 − 2 space-time dimensions. We may then derive DEs of the form: whereM (x, ) is a block-triangular matrix. Each block is associated with a sector of integrals. If we denote such a sector by f i (x, ), we can decompose the DEs in the form: where M i (x, ) denotes the diagonal block ofM (x, ) corresponding to the sector i, and B i (x, ) g i (x, ) captures the off-diagonal terms. One can then expand the integrals and matrices in : and solve the system order by order in . For a given basis, the condition that M i (x, ) is finite in is not always manifest. We obtain such a form by rescaling individual master integrals with powers of , until the matrix is finite. For further ideas and software to help facilitate the choice of MIs, see Refs. [50][51][52][53][54][55][56][57].
The DE system in eq. (2) fixes the master integrals up to some boundary conditions. It turns out that, in the case of our automated DE approach, a convenient choice for the boundary terms are MIs which are finite in the dimensional regulator . We use the package Reduze [58][59][60][61] to identify these MIs. They can be evaluated efficiently for Euclidean kinematics using the method of SD, since only a small number of sectors is needed for finite integrals and no contour deformation is required to avoid Minkowskian thresholds. We employ the package pySecDec [30,31] for this purpose. The derivation of a DE system is done with the help of the IBP reduction program Kira [62][63][64][65]. With the boundary terms fixed numerically and the DE system derived analytically, we transport the Euclidean point to the Minkowski point with the aid of the method of series expansions of the DE system [66][67][68] [69] as implemented in DiffExp [70].
As demonstrated in Fig. 1 we may choose different Euclidean points to fix the boundary terms numerically. This allows us to obtain a numerical error estimate of our automated method by taking the difference of two generated results for the same final Minkowski point. A more detailed discussion of the error estimate is provided in the supplemental material.
Typically, the transport from the Euclidean boundary point to the physical Minkowski kinematics requires several steps since the convergence radius of the series expansion at the boundary point is not large enough to reach the target point. The program DiffExp automatically determines the convergence radius and the number of required transport steps.
In general, the complexity of the multi-loop computation increases with the number of loops and independent scales and the number of MIs involved. In our automated approach, the largest investment of computing resources is required for the IBP reduction with Kira and the numerical evaluation of the boundary terms with pySecDec. However, the former needs to be done only once for a given Feynman integral family, and the latter only once for a given choice of mass parameter values. The transport to the Minkowski region with DiffExp is very fast, so that one can easily evaluate results for multiple different kinematic points, as needed e.g. for phase-space integrations. Quantitative information on the run time for our approach is given in the supplemental material. There the reader can also find a description for how our method can be extended to problems with multiple timelike momentum scales. III. Results and discussion. To demonstrate the power and broad applicability of our method, in the following and in the supplemental material, we present examples for 3-loop self-energy and vertex integrals and 2-loop box integrals. As discussed in the introduction, these are all examples of key theory ingredients for the physics program of future e + e − colliders and/or the HL-LHC. The 3-loop integrals are needed for currently unknown third-order corrections to electroweak precision observables connected with Z-boson producton and decay, whereas two-loop box integrals are important to improve the precision of several 2 → 2 processes, such as W + W − , ZH or ff production [71].
The technique described in this article allows one to compute the desired integrals to, in principle, arbitrary order in the dimension regularization parameter = (4 − D)/2 with multi-digit precision. To achieve a certain order k , some boundary terms need to be evaluated to higher orders k > k in . The required order k is determined automatically from the IBP relations.
For the examples shown below, some simple boundaryterm integrals have to be computed to O( 7 ), whereas no more than O( 3 ) is needed for more complicated boundary terms. When evaluating the boundary terms with SD as implemented in pySecDec, the computing time grows approximately linear with the order in .
All of the following numerical examples are based on the input parameters given in the supplemental material.  (4) and (8), respectively. W, Z and t stand for the W-boson, Z-boson and top quark, respectively. Example 1. As part of the 3-loop O(α 2 α s ) corrections to electroweak precision observables, one encounters the following scalar non-planar self-energy integral with eight propagators and only one massive W-or Z-boson internal line [72] (see Fig. 2 left): where Dq n ≡ d D qn iπ D/2 and a = W, Z. This example, for the parameter point p 2 = M 2 Z and M a = M Z belongs to a group of integrals which are difficult to evaluate with SD due to threshold effects. Using pySecDec with 10 7 integration points we obtain a result with less than two digits precision: 1, 1, 1, 1, 1, 1, 1, 1 Increasing the number of integration points does not improve the accuracy substantially. On the other hand, pySecDec can deliver accurate results for Euclidean parameter points, p 2 < 0, which are used as boundary terms for our automated DE transport. We thus obtain stable and precise results at the physical point: Here and in all the following results, we show all significant digits, i.e. the numerical error only affects digits beyond the ones shown in the equations. The error estimation will be described in more detail in the supplemental material. The integral family I lhNp1 (4) involves 30 master integrals and is considered simple in the context of our method.
Example 2. The next example is a family of 3-loop vertex integrals with one massive top quark and two massive W-boson propagators, see Fig. 2 (right), defined as where p = p 1 + p 2 and p 2 1 = p 2 2 = 0. These integrals also appear in so far unknown O(α 2 α s ) corrections to Zpole electroweak precision observables, constituting their most difficult parts.
With pySecDec we are unable to obtain a numerical result for the Minkowski point p 2 = M 2 Z . The problem already starts with the contour deformation which is necessary for SD with Minkowski kinematics and which fails to complete in a reasonable time. Similar to the SD method, the MB technique fails to deliver high-accuracy results for the considered integrals for p 2 = M 2 Z .
Using our automated DE transport method, the calculation requires the numerical evaluation of 77 master integrals with Euclidean kinematics, p 2 < 0, for the boundary terms. For the purpose of the present example, they have been evaluated with pySecDec to 10-digit accuracy. After the transport to the physical point p 2 = M 2 Z , we get at least eight significant digits for integrals of the family (8) up to tensor rank-3 (i.e. −3 ≤ a 10 + a 11 + a 12 ≤ 0). We here give numerical result for one rank-3 case: 1, 1, 1, 1, 1, 1, 1, 1 Additional examples, a 3-loop self-energy diagram with many massive propagators and a two-loop box diagram with four scales, are discussed in the supplemental material. IV. Summary and Outlook. In this work, we have proposed an efficient and versatile approach for the evaluation of a wide class of massive multi-loop, multi-scale Feynman integrals numerically, with typically 8 or more digits precision. It is based on the method of DEs with boundary terms specified for Euclidean kinematics, which are transported to the physical Minkowski kinematics using series solutions of the DE. The Euclidean boundary term integrals avoid all threshold singularities and thus can be straightforwardly evaluated numerically. Our implementation combines the public programs Kira, Reduze, pySecDec and DiffExp in a way that allows us to automatically construct the required integral families and the transport from the Euclidean boundary point to the physical kinematic point.
In principle, the technique can be extended to higher numerical accuracy and to wider classes of integrals with more loops and more external legs. A major bottleneck are the integration-by-parts (IBP) reductions that are needed to construct the DE system. A significant speedup of this step is achieved when using numerical values for the relevant mass and kinematic parameters. In addition, the evaluation of the boundary terms for Euclidean kinematics can be time-consuming if a high level of precision is required. Fortunately, there are ongoing improvements to the sector decomposition (SD) and Mellin-Barnes (MB) methods, see e.g. Refs. [31,73] and [74].
It is worth mentioning that the 3-loop examples presented in this article are very difficult to solve with existing analytical techniques (e.g. using IBP and DE) and general numerical methods (such as SD or MB methods). The proposed new technique is sufficiently general to provide the foundation for the computation of the required 3-loop corrections needed for electroweak and Higgs precision studies at the HL-LHC and future e + e − colliders, which are key elements of the physics program of these machines [9,22] Note the tenth propagator in (10) is linearly dependent and can be written in terms of the first nine propgators.
We included it as an auxiliary propagator to the definition to improve the readability of the final results in (11)- (15). We consider a family of box integrals that are part of the O(α 2 ) corrections to massive eµ scattering [75], see Fig. 4: where The I box2l integral family involves 55 master integrals. The input parameters for the box2l topology are s = 2, t = 5, m 2 1 = 4 and m 2 2 = 16.

Tracing the run time of the calculations
For the calculations demonstrated in this paper it makes sense to categorize the run time in three parts. The first part includes the preparation of the DE system suitable for later use with DiffExp. For this step, the major bottleneck is the run time of Kira. For the most complicated example vtwPl we need 4 hours to construct the DE system on a 12 core, 2.7 GHz Intel Xenon processor cluster with 128 GB of RAM. The second part is the computation of Euclidean boundary terms, which we evaluate with the sector decomposition program pySecDec. Since the boundary terms are finite integrals by construction and are free of thresholds in the Euclidean region, the Monte-Carlo (MC) or quasi-Monte-Carlo (QMC) scaling rules apply straightforwardly. For QMC the theoretical bound on the error is proportional to O(1/N ), where N is the number of points. For traditional MC, the bound is O(1/ √ N ). For the most complicated example vtwPl, we ran the QMC integrator configured with a maximum of 10 7 points. We ran the computation on a machine equipped with a 16 core Threadripper Pro 3955WX. Two different points in the Euclidean region were chosen, and both took approximately 3 days to complete. The maximum relative error reported by pySecDec was approximately 2 · 10 −9 . We expect that the precision of our result can be further improved by running on a more powerful machine, or by allocating a longer running time. The final part of the computation is the transport of the boundary terms to the Minkowski regions. The run time of this part of the computation depends on the size of the differential equations, the order of the homogeneous parts of the differential equations, and the number of line segments. In our case we are interested in the evaluation of the Feynman integral in one particular point. With DiffExp this is accomplished in about 3.5 hours on a single CPU core for the most complicated example vtwPl. Note that once the computation with DiffExp is finished the integral is available for a fast evaluation (in terms of milliseconds) at an arbitrary point along the line of the transport from Euclidean to Minkowski point.

Tracing the error of the calculations
The numerical error estimation has two main ingredients, the numerical series expansion of the differential equation (DE) system with DiffExp and the numerical evaluation of the boundary terms with pySecDec. The numerical error from the series expansion can always be rendered negligible compared to the error from the boundary terms by evaluating the expansion to sufficiently high order. On the other hand, the numerical errors of the initial boundary terms can usually be directly mapped to the final result. However, instead of relying on the error estimate from pySecDec, we verify the accuracy by carrying out separate transports from two different initial boundary points to the same final Minkowski point and taking the difference as an error estimate. This cross-check was sufficient for the examples presented here. Also, if necessary, it is relatively straightforward to increase the accuracy of the boundary terms by using more Monte-Carlo integration points, since the SD integrand is well-behaved in the Euclidean region.

Iterative approach for multi-scale kinematics
In the case of problems with multiple time-like momentum scales, one may use our method in an iterative way. For instance, we may derive the DE system just in one scale m 1 , while setting all other scales {m i } \ m 1 to numerical values, to simplify the IBP reduction. Then one can carry out the transport with respect to m 1 from the numerical point a 1 , where the boundary conditions are supplied with all scales taken to be Euclidean (spacelike), to a new point of interest denoted as a 2 , where m 1 is time-like. Next we derive a new set of DEs with respect to the next scale m 2 , and again we set all other scales {m i } \ m 2 to numerical values, and transport from the point a 2 to a new point of interest a 3 . One can continue in this way until all (n) relevant momentum scales have been transported from the Euclidean to the time-like domain. The downside of the iterative approach is that one has to repeat the IBP reduction whenever we want to reach a new kinematic point that differs in any of the n−1 first transported variables. Thus new improvements in the IBP reduction programs are highly anticipated.