TVID 2: Evaluation of planar-type three-loop self-energy integrals with arbitrary masses

We present TVID 2, a program to numerically evaluate an important class of planar three-loop self-energy master integrals with arbitrary masses. As with the predecessor version (TVID 1) the integrals are separated into a known piece, containing the UV divergencies, and a finite piece that is integrated numerically, implemented in C. The set of master integrals under consideration was found with self-energy diagrams containing two closed fermion loops in mind. Two techniques are employed in deriving the expressions for the finite pieces that are then numerically integrated: (a) Sub-loop dispersion relations in the case of topologies containing sub-bubbles, and (b) a modification of the procedure suggested by Ghinculov for integrals with only sub-loop triangles.


Introduction
The calculation of higher-order radiative corrections is important for the interpretation of precision measurements at the LHC and various e + e − machines, such as SuperKEKB and planned future Higgs and Z factories. Multi-loop contributions in the full Standard Model or models beyond the Standard Model (BSM) are particularly challenging due to the presence of many independent mass and momentum scales [1]. General loop integrals beyond the one-loop level cannot be solved analytically in terms of elementary functions. This observation prompted the investigation of new classes of special functions, such as harmonic polylogarithms [2], generalized harmonic polylogarithms [3], and elliptic polylogarithms [4], see e. g. Ref. [5] for a recent review. However, it is not clear if any multi-loop integral can be represented by these classes of functions, in particular beyond the two-loop level.
This motivates the development of numerical methods for multi-loop integrations. Two general approaches, which in principle can be applied to any number of loops and external legs, are known: sector decomposition and Mellin-Barnes representations. The former has been realized in the SecDec [6] and FIESTA [7,8] software packages, while the latter is the basis for the AMBRE/MBnumerics project [9]. Both approaches provide an algorithmic procedure for removing UV and IR singularities, but they require very significant computing resources, especially for integrals with physical internal thresholds that develop imaginary parts. Alternatively, more efficient numerical integration methods can be developed for limited classes of multi-loop integrals (see Ref. [1] for a review of some of these methods). At the three-loop level, an important step in this direction was achieved with the programs TVID [10,11] and 3VIL [12], which can evaluate the master integrals for arbitrary three-loop vacuum integrals.
This article reports on the new version 2.0 of TVID, which includes a large class of three-loop self-energy master integrals. This class consists of the master integrals necessary to evaluate three-loop self-energy diagrams containing two closed fermion loops. They are descendants of the planar three-loop self-energy topology (i. e. the ladder topology). It is shown that all these master integrals can be evaluated in terms of at most two-dimensional numerical integrals, by making use of the following two ideas: • Integrals with sub-loop self-energies are evaluated by using dispersion relations for the sub-loop. This technique was previously developed for the evaluation of two-loop self-energy master integrals [13,14].
• Planar integrals with without sub-loop self-energies are tackled with a variant of the method introduced in Ref. [15].
See section 2 for more details on the master integrals that fall into either of these two categories. The construction of the numerical integral representations for these master integrals is illustrated for a few characteristic examples in section 3. It should be noted that many of the master integrals are UV divergent, and these singularities must be removed before the numerical integration can be carried out. In TVID, this is achieved by subtracting terms that have the same singularity structure, but that lead to simpler integrals which are already Here ǫ = (4 − D)/2 and D is the number of space-time dimensions in dimensional regularization. Furthermore, the ν k are integer numbers which can be 0, 1 or 2 in our case.
To define our set of master integrals, we generated the diagrams with the topology of U 8a that occur in the three-loop self-energy diagrams with two closed fermion loops using FeynArts 3 [16]. We then performed an integral reduction based on integration-by-parts identities with the help of FIRE 5 [17]. The resulting set of irreducible three-loop master integrals is shown in Figs. 1 and 2. They are all of the form given in eq. (1) with ν k ∈ {0, 1, 2}. There are additional master integrals that factorize into products of one-loop and two-loop integrals. The complete set of the latter is shown in Fig. 3 (see also Ref. [18]).
We do not claim that this set of master integrals is minimal or optimal, but it is suitable for numerical evaluation in terms of two-dimensional numerical integrals, as will be demonstrated below. It should be emphasized that additional master integrals are needed for diagrams that do not conform to the planar master topology U 8a , such as non-planar three-loop self-energy diagrams. The integrals in Figs. 1 and 2 can be divided into two groups: • Integrals with one-or two-loop sub-loop self-energy. These can be evaluated efficiently using a dispersion relation for the sub-loop bubbles [10,13,14].
• Integrals without sub-loop self-energies. For these we employ a variant of the method proposed in Ref. [15]. This category comprises the master integrals U 7a , U 8a , U 7a1 and U 7a2 . All the remaining master integrals belong to the former category.

Examples
In the following subsections, our approaches for the numerical evaluation of the master integrals are described in more detail for a few characteristic examples from both categories.

Double-bubble integrals: U 5b
A basic one-loop self-energy sub-loop can be expressed in terms of a dispersion relation [13], e. g.
where ∆B 0 is the discontinuity of the one-loop function B 0 . In D = 4 dimensions, it is given by where λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + xz) and Θ(t) is the Heaviside step function.
U 5b contains two such sub-loop bubbles. Inserting the dispersion for each of them, one obtains Here T 3a is a two-loop self-energy function, see Fig. 3. The s 1 and s 2 integrals in eq. (4) diverge at the upper integral limit ∞, which can be attributed to the fact that U 5b is UV divergent. The expression can be rendered UV finite by subtracting suitable terms in the integrand that have the same UV singularity structure, but that are otherwise simpler than the full U 5b function. One way to achieve this purpose is by subtracting the first two terms in a Taylor expansion in p 2 : Here the prime in B ′ 0 etc. denotes a derivative with respect to p 2 . The integral in the last two lines of eq. 5 is now finite and can be evaluated numerically. U 5b (0, ...) and U ′ 5b (0, ...) are three-loop vacuum integrals, for which general methods for numerical evaluation are known [10,12]. Similarly, the two-loop function T 3a can be easily determined using the technique of Ref. [13]. The basic one-loop function B 0 is known analytically [19] (see Ref. [10] for expressions that use the same conventions as in this paper).

Planar master topology: U 8a
An simple method for numerically evaluating the master topology U 8a was presented in Ref. [15]. This integral is UV finite and thus can be computed in four dimensions. It can be written as where C 0 is the basic one-loop vertex function, which is known analytically in terms of logarithms and dilogarithms [19]. By moving to the center-of-mass frame for p and integrating over the solid angle of q, this becomes where This formula suggests that it is convenient to adopt x and y as integration variables, leading to the two-dimensional integral Here the Heaviside Θ function is inserted to ensure that the integral runs only over kinematically allowed values of x and y.
The integrand in eq. (9) has singularities at x = m 2 4 and y = m 2 5 , which lead to difficulties for numerical integration routines. In Ref. [15] this was addressed by using a deformation of the integration contours into the complex plane. Here we instead split the integrals into a residuum contribution and principal value integral, according to the prescription [10] ∞ This has the advantage that one does not need to worry about the complex contour crossing any other singular points.

Planar 7-propagator topology: U 7a
In princple, the 7-propagator integral can be treated with the same approach as described for U 8a in the previous sub-section. However, it turns out that the y integration is badly converging in this case. A better convergence behavior is achieved by using x and q 0 as integration variables, For the x integration again one can use the split into a residuum contribution and principal value integral according to eq. (10).

Implementation in TVID 2
The TVID 2 package has two components: • One component runs in Mathematica and performs the separation of the master integrals into UV-divergent subtraction terms and finite remainder functions, as described in section 3.1 and appendix A. This separation can be performed algebraically (keeping the momenta and masses as non-numerical symbols) or with numbers for the momenta and mass inserted from the beginning (which will speed up the evaluation in Mathematica).
• The second component carries out the numerical integration of the finite remainder functions (i. e. the functions labeled U ...,sub in appendix A). It is written in C and uses an adaptive Gauss-Kronrod algorithm for the integrals, which yields a relative precision of 9-10 digits for most cases (see below for exceptions to this statement). The input and output are handled through simple text files that contain a list of numerical parameter values. The Mathematica component of TVID 2 can directly call the numerical C component through an external system call.
The numerical of TVID uses quadruple precision floating points numbers to reduce roundingoff errors in the tails of the integrals. However, for some master integrals this turns out not to be sufficient. For these cases, we use an asymptotic formula for the integrand in the limit of large values of the dispersion variable s 1 and/or s 2 . The asymptotic formula is used for values of s 1 + s 2 > s cut , with a suitably chosen value for the parameter s cut . Specifically, s cut = c×p 2 , where c is a constant that depends on which function U ...,sub is being considered. Since s cut is proportional to p 2 , there could be a loss of precision for cases when p 2 is either much larger or much smaller than the masses in the integral.
The reader should take note of the following limitations of version 2.0 of TVID: • The program cannot handle IR-divergent integrals. IR divergencies may occur from certain configurations with multiple massless propagators or threshold singularities. TVID 2.0 furthermore does not check whether a certain parameter choice leads to an IR divergency; the user has to ensure that this is the case.
• There are additional cases (i. e. particular combinations of input parameters) that are IR finite but may require a special treatment within TVID to avoid numerical instabilities. A few of these are implemented in version 2.0, but there are probably many more that are currently missing. The authors encourage users to submit any such special cases when they discover them, and they will be considered for implementation in future versions of TVID.
• The finite remainder functions of U 6m1 , U 6m3 , U 6n2 , U 7a1 and U 7a2 are related to those of U 6m , U 6n and U 7a through mass derivatives: In version 2.0 of TVID these have been implemented using numerical differentiation (based on a five-point stencil) * . As a result the delivered precision for these functions is reduced to 6-7 digits.
• Due to the presence of C 0 functions in the integrands of U 7a , U 8a , U 7a1 and U 7a2 , their evaluation is much more time-consuming than that other master integrals. To mitigate this issue, TVID 2.0 uses only double precision floating point numbers for the evaluation of the C 0 functions, and the target precision of these master integrals is reduced to 8 digits for U 7a and 6 digits for U 8a , U 7a1 and U 7a2 . • The algebraic part in TVID 2.0, which runs in Mathematica, performs the separation of UV divergent subtraction terms (as detailed in appendix A) for individual master integrals or for expressions that contain any linear combinations of these. However, the resulting expressions can grow rather large, in particular if the masses are treated symbolically at this stage. The Mathematica code in TVID 2.0 is not optimized to deal with very large expressions, and the user may have to modify the PrepInt method in TVID to avoid excessively long computing times and related problems.
As a check and to calibrate the performance of TVID 2, we have performed comparisons with FIESTA 4.1 [8]. If one takes each master integral in isolation, the subtraction terms defined in appendix A would require the evaluation of some two-loop functions up to O(ǫ). As already mentioned above, these T ...,delta functions are currently not implemented in TVID 2.0, and they would also not be needed in the calculation of physical quantities. To circumvent this issue, we have such defined modified version of some master integrals, where certain terms have been subtracted that reflect the physical counterterm structure. These modified master functions are listed in Tab. 1. Note that the achievable numerical precision in some cases is somewhat reduced due to cancellations between different terms in these expressions. The benchmark tests have been performed on a single core of an Intel Xeon CPU with 3.7 GHz. Two parameter choices have been considered: * The reason for this choice is that the reduction formula for the mass derivative T 5a has rational coefficients with high polynomial degrees in both the numerators and denominators, which leads to many instabilities of the 0/0 type. Furthermore, the direct integration of U 7a1 and U 7a2 involves very large numerical cancellations between regions with positive and negative integrand, which leads to numerical instabilities.   The results are shown in Tab. 3.
As can be seen from the tables, there is generally excellent agreement between TVID 2 and FIESTA within integration errors. Only for the cases with 7 or more propagators (U 7a and U 8a ) that require contour deformation in sector decomposition (Tab. 3), the discrepancy between the two programs is larger than the integration error reported by FIESTA. This may in part be due to the contour deformation being unable to make the integrands sufficiently smooth in those cases. TVID 2 generally achieves 6-10 digit precision within run times ranging from less than 1 second to about 20 minutes. The precision and run time are not crucially affected by the presence of physical cuts (i. e. whether p 2 is below or above any of the thresholds of the integral). Note that the run times shown in the tables only reflect the time for the numerical integrations. Additionally, the preparation time for the integrals in the Mathematica module of FIESTA can be significant, in particular for the cases that require contour deformation (Tab. 3) † .

Conclusions
Numerical integration is currently the most efficient way to evaluate the finite pieces of multi-loop integrals with arbitrary masses, which set a multitude of different scales. The program TVID 2 aims to provide an efficient and automizable procedure for the numerical evaluation of a large class of three-loop self-energy integrals.
The master topologies fall into one of two categories. Topologies containing one or two sub-loop self-energies are generally UV-divergent and are treated by subtracting simpler, known integrals with the same divergence structure. The remaining finite pieces are written as one-or two-dimensional integrals over analytically known functions with the help of dispersion relations, which in turn can be evaluated numerically. The topologies without sub-loops self-energies are UV-finite and can be performed as two-dimensional integrals containing one-loop triangle functions in integrand. The general procedure for this class has been described in Ref. [15] and was adapted to avoid complex contour deformation. Several technical subtleties when following this approach are discussed in this paper. TVID 2 contains also contains all the basic elements of the S2LSE package for two-loop self-energy master integrals [20].
In order to confirm the correctness and accuracy of the implementation we carried out a multitude of independent comparisons. To that end we utilized the publicly available package FIESTA and found excellent agreement for almost all the master integrals. If the external momentum squared is sufficiently large so that the master integrals develop a non-zero imaginary part, the precision and accuracy of FIESTA is significantly diminished, leading to † In fact, the preparation of the contour deformation for U 8a in FIESTA takes several days.   less perfect agreement between the two programs in some cases. Generally, TVID 2 achieves 6-10 digit precision for all the master integrals in run times of seconds to minutes. At the present time a few issues remain to be addressed before one has complete computational control over the entirety of three-loop self-energy-type integrals. The current version of TVID is not equipped to treat integrals that exhibit soft/collinear divergencies, which can occur in some master integrals for certain input parameter combinations. Furthermore, even though we cover a large subset of master integrals, there is a number of topologies missing, specifically the descendants of the "Mercedes star" and the non-planar topologies. These issues are delegated to future work.
Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. A. F. is grateful to J. Gluza and other members of the Polish National Science Centre grant no. 2017/25/B/ST2/01987 for discussions and exchanges related to this topic.

A Subtraction of divergent terms
Before the master integrals in Figs. 1, 2 and 3 can be evaluated numerically, one needs to remove their UV divergencies. This can be achieved by subtracting simpler integrals that have the same UV singularity structure, but that are known in the literature. The finite remainder parts are denoted by T xxx,sub and U xxx,sub in the equations below. These can be evaluated with the numerical part of TVID 2.
For the sake of brevity, the following shorthand notations are used in this section: denotes the B 0 function with the order-n Taylor expansion subtracted, where λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + xz) , and Θ(t) is the Heaviside step function.
Copyright: The TVID source code may be freely used and incorporated into other projects, but the authors ask that always a reference to this document and to Ref. [10] be included.
External code elements: TVID includes the Gauss-Kronrod routine QAG from the Quadpack library [24], translated into C++, and the C++ package doubledouble for 30 digit floating point arithmetic [25]. It further requires the package LoopTools [26], version 2.10 or higher, which must be installed by the user separately.
Code availability: The TVID source code is available for download at http://www.pitt.edu/~afreitas/.

B.1 Numerical part
The numerical part of TVID is programmed in C and evaluates the finite remainder functions defined in appendix A. It is called with the command ucall infile outfile where infile is the name of the input file, and outfile is the name of the file where the results shall be placed. By default, ucall is located in the subdirectory ccode. asks for the evaluation of U 4,sub (1, 2, 3, 4) and of U 5a,sub (20, 1, 1, 1.5, 2, 2). When ucall is completed, it fills outfile with a list of the numerical results, again separated by line breaks. For instance, the example above will return -5.555128856244808e1 0.0 0.306188821751692 -6.207131465925367 Here the first and second number in each row are the real and imaginary part of the result, respectively. The option -e allows the user to also receive information about the integration error: In this case, a third number is added to each row in outfile, which provides the integration error. For the example above, one obtains -5.555128856244808e1 0.0 0.771470006356965e-10 0.306188821751692 -6.207131465925367 0.564457354626832e-11 Internally, the numerical code uses the Gauss-Kronrod routine QAG from the Quadpack library [24] to evaluate the dispersion integrals. This routine has been translated into C++ from the original FORTRAN code, and amended to facilitate 30 digit floating point arithmetic from the package doubledouble [25].  Table 4: Symbols for basic finite remainder functions used in numerical and algebraic parts of TVID 1 [11].

B.2 Algebraic part
The algebraic part of TVID runs in Mathematica 10 [23] and performs the separation of divergent and finite pieces of the master integrals. The program is loaded in Mathematica with << mcode/i3.m The two main user functions are PrepInt and UCall (and the variant UCallE of the latter).
PrepInt takes as input any master integral in Figs. 1, 2 and 3, as well as the three-loop vacuum master integrals U 4 , U 5 or U 6 , or a linear combination thereof. It returns a series expansion in ǫ, whose coefficients contain the finite remainder functions listed in Tabs   Here the directive SetPrecision has been used to mitigate numerical rounding errors within Mathematica. PrepInt can also be called with symbols for the momentum and mass parameters, e. g. PrepInt[U4a[ps,m1s,m3s,m6s,m8s]], although this can lead to fairly large expressions.
UCall invokes the numerical code ucall (see previous subsection) to evaluate the finite remainder functions in the output of PrepInt. For the example above this leads to Technically, the executable ucall is called through an external operating system command, using the Mathematica function Run. The function UCall looks for the executable ucall in the subdirectory ccode of the TVID installation. If the user places ucall in a different directory, the variable $Directory in mcode/i3.m must be adjusted. For passing input and output to and from the executable, UCall uses the filenames specified in the variables $FileIn and $FileOut, respectively. In most cases, the user will not need to change any of these global variables.
The variant UCallE also returns integration errors for the various numerical master functions (by using the option -e when calling ucall). For the example above this yields

B.3 Installation
TVID 2 comes in a compressed tar archive. After saving it in the desired directory, it can be unpacked with the command tar xzf tvid.tgz The program contains the following subdirectory structure: ccode the C/C++ files for the numerical part of TVID 2; ccode/doubledouble the doubledouble for 30 digit floating point arithmetic [25]; mcode the Mathematica files for the algebraic part of TVID 2, as well as examples.
Before compiling, the user must ensure that LoopTools version 2.10 or higher is installed on the system. It may be necessary to adjust the LoopTools path in ccode/makefile.
To compile the numerical C part of TVID, execute the commands cd ccode make The make file provided has been tested on Scientific Linux 6. It makes use of the fcc script included in LoopTools, which should help to facilitate compilation on a range of UNIX-type operating systems. The authors cannot guarantee that the installation process is successful on any operating system, but they appreciate any helpful suggestions, comments and bug reports.