Numerical Multi-Loop Calculations via Finite Integrals and One-Mass EW-QCD Drell-Yan Master Integrals

We study a recently-proposed approach to the numerical evaluation of multi-loop Feynman integrals using available sector decomposition programs. As our main example, we consider the two-loop integrals for the $\alpha \alpha_s$ corrections to Drell-Yan lepton production with up to one massive vector boson in physical kinematics. As a reference, we evaluate these planar and non-planar integrals by the method of differential equations through to weight five. Choosing a basis of finite integrals for the numerical evaluation with SecDec3 leads to tremendous performance improvements and renders the otherwise problematic seven-line topologies numerically accessible. As another example, basis integrals for massless QCD three loop form factors are evaluated with FIESTA4. Here, employing a basis of finite integrals results in an overall speedup of more than an order of magnitude.


Introduction
In the early days of hadron collider physics and Quantum Chromodynamics (QCD), the Drell-Yan process was recognized to be of fundamental importance [1] and, over the years, it has received much attention. The electroweak corrections of relative order α were calculated in [2][3][4]. Another milestone of particular interest was set long ago when the corrections of relative order α 2 s were calculated in references [5][6][7][8]. The fact that these corrections were found to be relatively large in size motivates one to consider the subleading terms in the two-loop perturbative expansion of the Standard Model Drell-Yan production cross section as well. The most important class of subleading two-loop contributions are those of relative order αα s , what we shall hereafter refer to as mixed Electroweak-Quantum Chromodynamic (EW-QCD) corrections. The mixed Electroweak-Quantum Chromodynamic corrections of relative order αα s have already been studied in certain approximations. First, in reference [9], the virtual corrections without propagating W or Z bosons were calculated. More recently, the full set of two-loop corrections were considered in the single massive gauge boson resonance region [10,11]. Progress on real emission contributions and the infrared structure has been presented in [12][13][14]. Recently, the planar two-loop master integrals with up to two same-mass vector bosons in the loops were calculated in the Euclidean region [15]. The non-planar master integrals relevant to the problem are vertex integrals, which have been available in the literature for some time [16,17].
In this paper, we consider the two-loop master integrals for mixed EW-QCD corrections to Drell-Yan production with up to a single massive vector boson exchanged and focus on physical kinematics. The integrals may be conveniently computed using the method of differential equations [18][19][20][21][22][23][24]. We employ a normal form basis [25][26][27], integrate the -expanded differential equations in terms of multiple polylgarithms [28], and impose regularity conditions to fix the boundary constants. In some cases this is supplemented by explicit solutions worked out directly from Feynman parameters to all orders in the parameter of dimensional regularization, . Analytical solutions for all of our -expanded integrals are included with our arXiv submission through to weight four. 1 Although the weight four solution is one of the primary goals of our study, we found it interesting to calculate to one order higher in than necessary to facilitate checks in the second part of this work.
To cross-check an analytical solution for a multi-loop integral or to even produce the primary result for a phenomenological application, numerical evaluations using sector decomposition programs [31][32][33][34] are very valuable. One might hope that one could just pick any integral basis, run one of the available programs for a physical point in phase space, and then obtain a reasonably accurate and precise result in, say, a few days' time. However, for our mixed EW-QCD integrals we found this not to be the case in practice. Instead, picking a basis which is free of infrared and ultraviolet divergences allows us to obtain reliable and relatively fast results with SecDec 3 [33]. This method of finite integrals introduced by Erik Panzer and the authors [35,36] is generic and straightforward to automate. The approach has successfully been employed in the calculation of twoloop integrals for double Higgs production with exact top quark mass dependence [37,38] already. Here, we present a detailed numerical study of performance gains stemming from the use of this technique. We quantify the effect also for a very different setup: the evaluation of the massless three-loop form factor integrals with FIESTA 4, which performs particularly well for such Euclidean integrals. Again, we observe drastic improvements with regard to the numerical convergence.
The outline of this article is as follows. In Section 2, we comment on the Feynman diagrams relevant to the mixed EW-QCD corrections and define the integral families which we use to describe our master integrals. In Section 3, we define our normal form basis, present the differential equations they satisfy, and solve the integrals in terms of multiple polylogarithms. In Section 4, we present the results of our numerical analysis of the two-loop mixed EW-QCD Drell-Yan integrals and discuss the crucial performance enhancements obtained by employing a basis of finite integrals. We present our conclusions in Section 5. Finally, in Appendix A, we present a self-contained numerical study of the master integrals for massless three-loop form factors up to contributions of weight eight.

Integral Families
We consider neutral-current and charged-current lepton pair production in quark-antiquark annihilation, 3) Figure 1. Selected planar and non-planar two-loop Feynman diagrams for Drell-Yan lepton production. The heavy wavy lines denote a Z boson, the light wavy lines denote a photon, and the curly lines denote a gluon. where p 2 1 = p 2 2 = p 2 3 = p 2 4 = 0. As usual, the Mandelstam invariants are For the physical part of the phase space we have and selected Feynman diagrams are depicted in Figure 1 above. We are interested in the two-loop master integrals required for the calculation of the O(α s α) corrections to these processes. In this paper, we restrict ourselves to the subset of these integrals, which involve at most one massive vector boson, W or Z, propagating in the loops. We will use the generic symbol m for the vector boson mass. The relevant loop integrals can be indexed using the four integral families shown in Table 1. Integration by parts reduction [39][40][41] allows us to rewrite any Feynman integral in these families as a linear combination of master integrals. The reductions for this paper were performed with the implementation Reduze 2 [29,[42][43][44] and we adapt its conventions for the labeling of sectors, etc. An ideal choice for the master integrals will depend on the application and we will discuss suitable options for different applications in the following sections.

Analytical Solution from Differential Equations
We employ the method of differential equations to derive analytical solutions for the master integrals. Altogether, we find that we need to consider forty-nine master integrals for the system of differential equations. Some of these integrals are related by a permutation of external legs. They are covered by crossed integral families, which we will denote by an overline (e.g. we write F:x for the crossed version of sector x from family F). Note that, because our four-particle integral topologies are naturally a function of the invariants s and t (and never u = −s − t), we may without loss of intelligibility suppress all explicit momentum labels; instead, we use directed external lines and explicit function arguments to help clarify our basis integral definitions when necessary. As usual, we use dots to denote doubled propagators, heavy lines to denote massive propagators, and explicitly write all numerator insertions in square brackets. The starting point for the construction of our basis are the following integrals: (s) f C:214 Building upon these definitions, we construct a convenient basis of integrals for the method of differential equations along the lines discussed in reference [45] 19 2 We would like to point out further recent work on the analysis of master integrals and the construction of a normal form basis [46][47][48][49][50][51][52][53][54][55][56][57]. Note that these definitions repeat nine massless integrals to allow for an independent treatment of the differential equations for the massless integrals, m 1 , . . . , m 12 , and the one-mass integrals, m 13 , . . . , m 58 . We employ the integration measure with = (4 − d)/2, which renders our integrals m i functions of the dimensionless variables only.
In this basis, we obtain a system of differential equations in normal form [25][26][27] for the vector m = (m i ), with the letters l k (x, y) of the symbol alphabet and matrices A (k) of rational numbers. Expanding the masters integrals m i about = 0, the differential equations fully decouple and can be integrated order-by-order in in terms of multiple polylogarithms. The alphabet is linear in both variables, which implies that the integration path can be chosen along the x and y axis, resulting in an expression written in terms of Goncharov G functions with either x or y in the final argument. This representation can introduce spurious singularities and therefore may not be optimal for numerical evaluations [58].
Instead, one can employ Li functions of fewer but possibly more involved arguments. Recently, it has been demonstrated explicitly that ln, Li 2 , Li 3 , Li 4 , and Li 2,2 functions are sufficient to reduce a general G function of at most weight four [59]. In order to arrive at real-valued functions with a convergent power series representation in a specific region of phase space, it can be useful to allow for Li 2,1 functions as well. A representation written in terms of such functions can be obtained either by recasting a solution written in terms of G functions or by directly integrating the differential equations via an ansatz built out of the required functions. We perform these calculations using in-house routines based on the symbol and coproduct calculus [60][61][62][63].
Note that, for phenomenological applications, we are interested in the solutions of our integrals through to weight four. In the next section, however, we consider alternative bases of master integrals for the purpose of numerical evaluation which partially involve weight five functions at the order required for physics applications. In order to validate our numerical analysis of the solution, we therefore find it useful to obtain results through to weight five.
We construct an ansatz for our solutions in terms of the functions using the Duhr-Gangl-Rhodes algorithm [62]. We require the symbol of each function to not introduce letters beyond those already present in our differential equation (3.5) and we construct suitable power products of our letters (3.6) accordingly. For the arguments of the logarithms ln we choose As admissible arguments for the classical polylogarithms Li n , we obtain the following sixtysix power products of letters which have the property that not only the argument z but also 1 − z factorize over the original symbol alphabet. Extending this set by the argument 1, forming pairs and selecting those pairs (z 1 , z 2 ), for which 1 − z 1 z 2 factorizes over the alphabet, gives us 618 admissible arguments for Li 2,1 , Li 2,2 , and Li 3,2 functions which we do not list here for the sake of brevity. Similarly, we find 3342 admissible triples (z 1 , z 2 , z 3 ) as arguments for our Li 2,2,1 functions. While we have no proof that Li 2,2,1 functions are strictly required, we were successful in constructing the genuine weight five part of the solution only when including them in addition to Li 3,2 and Li 5 functions. This observation is independent of further constraints to be described below, which motivate the inclusion of Li 2,1 functions. Not all of these functions are needed and we formulate further objectives for our functional basis. We prefer functions which are real-valued for physical kinematics This requirement is the reason why we include Li 2,1 functions in our basis. In addition to requiring real-valuedness, we prefer functions which possess a convergent power series representation. Insisting upon such a representation automatically avoids additional manipulations which might otherwise be necessary to carry out numerical evaluations. Indeed, we find that it is possible to use functions which are real-valued over the entire physical region of the phase space with a single important exception. The letter l 4 = 1 + y changes sign at s = m 2 and the functional form of our results will change in this region if we insist upon separating real and imaginary parts explicitly. This, however, is to be expected because there is a physical singularity at the point s = m 2 , which would be regulated by the width of the massive vector boson inside the loop. Our solution, in fact, contains just a single function which is sensitive to the cut starting at y = −1: ln(−l 4 ) + iπ for y < −1. (3.11) In contrast, no logarithms of the letter l 6 = x − y appear in our results and all other elements of our functional basis which depend of l 6 are real-valued for both x > y and x < y. Employing our functional basis, we construct an ansatz for the solution and match it against the differential equations (3.5). We employ regularity conditions and require real-valuedness in the Euclidean domain to fix the integration constants. In addition, we turn out to be a bit more non-trivial.
As an additional sanity check, we found it convenient in some cases to look at asymptotic expansions of integrals in particular limits in order to explicitly confirm their scaling behavior. For this purpose, we found the Mathematica package asy.m [64] to be extremely useful. For the convenience of the reader, we provide our analytical solutions through to weight four in the ancillary files sol-phys-AB.m and sol-phys-CD.m on arXiv.org.

Numerical Analysis
In an effort to check an analytical calculation such as the one presented in the previous section, it is desirable to numerically evaluate the loop integrals with an independent method. In other cases, a numerical evaluation might even be the method of choice, e.g. if it is not clear how to evaluate the integrals analytically. A standard method for numerical evaluation is sector decomposition and one might think that it is straightforward to numerically check the integrals f i in (3.1) using this approach. Unfortunately, this is not the case. Specifically, for double box integrals with additional numerators or dotted propagators we often encountered cases where none of the available software packages is able to provide reliable results. In some cases, issues related to the presence of problematic singularity structures in the integrand can be tracked down and, with dedicated effort, cured [65]. For what concerns the two-loop mixed EW-QCD topologies considered in this paper, the most challenging cases for the code are seven-line integrals which have −4 poles and doubled propagator denominators.
As discussed in earlier work by the authors and Erik Panzer [35,66], it is possible to employ a basis of finite Feynman integrals defined by allowing for higher-dimensional integrals with potentially higher powers of the propagator denominators to extract all poles in analytically before performing any numerical integrations. 3 The change of basis is performed with integration by parts reductions, whose computational complexity is nonnegligible but still lower than that of the reductions for the amplitude. This is true both for the processes discussed in this work and for all other phenomenologically relevant problems which we have studied.
The finite integrals could, in principle, be evaluated by direct numerical integration. However, their integrands may still contain structures which spoil the convergence of the numerical integrations, such as integrable but large variations near the boundaries. In our experiments with a limited number of Euclidean integrals, we find that an additional sector decomposition step improves the performance of the evaluation. Furthermore, handling branch cuts requires a dedicated treatment for physical kinematics. We therefore find it convenient to employ existing sector decomposition programs for the numerical integration of our finite integrals.
We employ the sector decomposition program SecDec 3 to evaluate our two-loop topologies using either conventional or finite master integrals. As is apparent from Table 2, working with a basis of finite integrals improves the numerical performance tremendously. Let us now explain in more detail how we used the code to produce Table 2. Although we made a good effort to use the program in an appropriate way, it is certainly possible that other researchers might manage to produce better results in less time. However, we at least took care to use identical program settings and hardware while carrying out our numerical studies of the conventional and finite integral bases. 4 Altogether, we ran for roughly sixty days on the four cores of a i7-4940MX processor. For our finite integral evaluations, just over an hour was required; the excessive run time is due to the fact that, in spite of allowing the VEGAS routine [67] provided by CUBA [68] to sample up to 5 × 10 8 points, numerous conventional basis integrals failed to achieve the very modest goal of epsrel = 1 × 10 −2 and epsabs = 1 × 10 −4 at the physical phase space point p = {s → 17, t → −7, m 2 → 6241/1681}. Analytical expressions for the finite basis integrals were obtained by using the explicit connection between the finite integral basis and the normal form integral basis. We found good agreement in all cases. At worst, we recorded a fractional difference of 6 × 10 −3 at weight four between our exact results evaluated at p and the numerical results output by the program. The average relative accuracy at weight four was found to be 9 × 10 −4 , quite acceptable given our modest accuracy goal.
Due to the fact that there is a substantial loss of precision in rotating from the finite integral basis to the normal form integral basis, we found that our initial epsrel = 1×10 −2 and epsabs = 1 × 10 −4 run of SecDec 3 did not allow us to directly check the normal form basis to similar part per mille precision. To achieve this, we had to carry out a second, higher precision run of the program with epsrel = 1 × 10 −4 and epsabs = 1 × 10 −8 , which took roughly eleven days of run time. Let us stress again that it would be hopeless to attempt this exercise using the conventional integral basis.  Table 2: Numerical performance of finite and conventional integral bases for two-loop mixed EW-QCD corrections to Drell-Yan lepton production in the physical region using SecDec 3. For each integral in the above, the run time is given in seconds and the fractional difference from the analytical solution is given for the expansion coefficient which first gives rise to weight four multiple polylogarithms. In the final row of the table, total run times and worst-case relative accuracies are recorded for both integral bases.
We find that, in several other cases, significantly more mileage can be squeezed out of publicly available sector decomposition programs simply by working with a well-behaved, finite integral basis. As a second example, we present the evaluation of three-loop form factors in massless QCD in Appendix A. We employ the other publicly available program under active development, FIESTA 4, and observe dramatic gains in both speed and numerical convergence when using a basis of finite integrals.

Conclusions and Outlook
In this paper, we considered the two-loop integrals relevant to the αα s corrections to Drell-Yan production with up to a single massive vector boson exchanged. As a reference, we calculated their Laurent expansion through to weight five in terms of multiple polylogarithms using the method of differential equations. Our representation in terms of real-valued functions allows for fast and precise numerical evaluations over the entire physical region of the phase space. We found it challenging to even check our analytical solutions using available sector decomposition programs. Employing a basis of finite integrals systematically improved the situation and rendered all integrals numerically accessible with SecDec 3, both in Euclidean and physical kinematics. Order of magnitude improvements both in program run time and integration error were also found for massless three-loop form factor integrals using FIESTA 4 when using finite instead of conventional master integrals. For the finite integrals, we allowed for shifts to higher numbers of spacetime dimensions and additional powers of denominators; the actual change of basis is performed with integration by parts reductions in a highly-automated way. An implementation of the algorithm to construct a basis of finite integrals is publicly available in the package Reduze 2.1 on HepForge. Numerical evaluations along the lines discussed in this paper are not only interesting for cross-checks of analytical solutions, but may actually serve as the primary method for phenomenological applications in especially complicated cases [37,38].

A Numerical Analysis for Massless Three-Loop Form Factors
In this appendix, we present the results of our FIESTA 4-based numerical study of the finite and conventional integral bases for the massless three-loop vertex functions (i.e. the bases employed in references [66] and [71] respectively). They show the advantages of using a basis of finite integrals in a setup independent of what was employed in Section 4. Due to the fact that the three-loop form factor master integrals have been calculated to weight eight [66,[72][73][74][75], we can also study how well our finite basis performs for the evaluation of higher order coefficients in the expansion.
By experimenting with the code and discussing some aspects of it with the developer, we settled on the FIESTA 4 settings ComplexMode = False, STRATEGY = STRATEGY X, BucketSize = 28, NegativeTermsHandling = None, and CurrentIntegratorSettings = {{"mineval","1000"}, {"maxeval","500000"}, {"nstart","50000"}} for the trivial phase space point q 2 = −1. Especially crucial is the setting NegativeTermsHandling = None, which seems to substantially improve the performance for Euclidean Feynman integrals in general. We also configured the code to run on all four cores of a i7-4940MX processor. Once again, we used the VEGAS Monte Carlo integration routine provided by CUBA, but, this time, we found that Mathematica 8.0.4 is a better choice than Mathematica 9.0.1 for the front end.
Let us begin by discussing Table 3, where all twenty-two irreducible integral topologies in both the finite and conventional integral bases are evaluated through to the terms of weight six. Reasonable results are obtained in both integral bases using the program settings described above. However, in going from the conventional basis to our chosen finite basis, the program run time drops by more than a factor of fifty: from just shy of twenty-two hours to about half an hour. There is, however, another crucial difference between the two sets of integrals which plays an even more important role when one goes from three to four loops: for the fixed program settings given above, one finds that the numerical convergence of the VEGAS algorithm is far superior for the finite basis integrals. 5 In fact, the average relative error of the weight six contributions with respect to the exact results is 2×10 −5 for the finite integral basis compared to 2×10 −4 for the conventional integral basis, an improvement of an order of magnitude. Finally, we see from Table 4 that, in going from weight six to weight eight with finite integrals, the run time increases by only a factor of three (with a slight loss of precision for fixed program settings). 6 Given the long run times with conventional integrals already at weight six, we refrained from trying to perform a similar weight eight evaluation of them. 43 s 5.05 × 10 −5 5 Let us point out that, in recent work on massless form factors [66], we were able to numerically evaluate a finite non-planar twelve-line master integral through to weight eight to four decimal digit accuracy using FIESTA 4 in about one day on a desktop computer. 6 This shows that, to obtain acceptably fast and stable numerical evaluations of master integrals at higher orders in , our integral basis is an effective alternative to the basis suggested in [76], where master integrals were systematically chosen to have the property that their coefficients remain finite in the → 0 limit.  Table 3: Numerical performance of finite and conventional integral bases for massless three-loop form factors with FIESTA 4. In addition to the run times, the fractional difference from the exact value is shown. In the final row of the table, total run times and average relative accuracies are recorded.  Table 4: Numerical performance of finite master integrals with FIESTA 4 for higher orders in the expansion. On the left, the weight six results of Table 3 are reproduced. On the right, results for weight eight are shown.