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][4][5][6], and elliptic polylogarithms [7][8][9][10][11][12][13][14][15], see e.g. ref. [16] 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 [17][18][19][20] and FIESTA [21][22][23][24] software packages, while the latter is the basis for the AMBRE/MBnumerics project [25][26][27][28][29]. Both approaches provide an algorithmic procedure for removing UV and IR singularities, but they require very JHEP01(2020)024 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 [30,31] and 3VIL [32], 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 [33,34].
• Planar integrals with without sub-loop self-energies are tackled with a variant of the method introduced in ref. [35].
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 known in the literature, see section 3 and appendix A. TVID provides the integrated subtraction terms in the framework of dimensional regularization and then numerically evaluates the finite remainder integrals. More information on the implementation of the three-loop self-energy integrals in TVID 2.0 can be found in section 4, together with a discussion of potential problems and comparison to results in the literature. A manual for the installation and usage of TVID 2.0 is provided in appendix B.

Planar three-loop self-energy topologies and master integrals
This article focuses on the evaluation of the "planar-type" three-loop self-energy integrals that descend from diagrams containing two closed fermion loops. With "planar-type" we refer to topologies that can be considered descendants of the "master topology" U 8a in figure 1, by removing and/or doubling some propagators. For the set of master integrals, we choose only integrals without numerator terms, which generically are of the form Figure 2. Master integral topologies with doubled propagators considered in this paper. The dot indicates a propagator that is raised to the power 2.

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 [33], 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 figure 3.
The s 1 and s 2 integrals in eq. (3.3) 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 JHEP01(2020)024 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. (3.4) 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 [30,32]. Similarly, the two-loop function T 3a can be easily determined using the technique of ref. [33]. The basic one-loop function B 0 is known analytically [39][40][41] (see ref. [30] 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. [35]. 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 [39][40][41]. By moving to the center-of-mass frame for p and integrating over the solid angle of q, this becomes This formula suggests that it is convenient to adopt x and y as integration variables, leading to the two-dimensional integral . (3.8) Here the Heaviside Θ function is inserted to ensure that the integral runs only over kinematically allowed values of x and y.

JHEP01(2020)024
The integrand in eq. (3.8) has singularities at x = m 2 4 and y = m 2 5 , which lead to difficulties for numerical integration routines. In ref. [35] 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 [30] ∞ 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 (3.10) 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. (3.9).

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.

JHEP01(2020)024
The numerical of TVID uses quadruple precision floating points numbers to reduce rounding-off 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). 1 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 .

JHEP01(2020)024
• TVID 2.0 does not numerically evaluate the O( ) parts of the two-loop functions in figure 3 (labeled T ...,delta in appendix A). In principle, these O( ) terms are needed if one wishes to evaluate all the master integrals in figures 1 and 2 to O( 0 ). However, in the calculation of any physical observable the T ...,delta functions should drop out when including the appropriate counterterm contributions, so that their explicit numerical value should not be needed.
• 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 [24]. 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 table 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: a) One choice where p 2 < m 2 i , so that p 2 is below any threshold, and all master integrals are real. For these cases FIESTA has been run with the settings CurrentIntegratorSettings = {{"epsrel","1.000000E-05"},{"maxeval","5000000"}}; The results are shown in table 2.
b) A second choice where p 2 m 2 i , all master integrals develop an imaginary part. For these cases FIESTA has been run with the settings CurrentIntegratorSettings = {{"epsrel","1.000000E-04"},{"maxeval","5000000"}}; The results are shown in table 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 (table 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 (table 3). 2

JHEP01(2020)024 5 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. [35] 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 [42].
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 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.

A Subtraction of divergent terms
Before the master integrals in figures 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) , (A. 6) and Θ(t) is the Heaviside step function. The UV subtractions for the two-loop integrals can be taken over from ref. [42]: (A.14)
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. [30] be included.

External code elements. TVID includes the Gauss-Kronrod routine QAG from the
Quadpack library [47], translated into C++, and the C++ package doubledouble for 30 digit floating point arithmetic [48]. It further requires the package LoopTools [49], 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.
infile 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: ucall infile outfile -e 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 [47] 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 [48].

B.2 Algebraic part
The algebraic part of TVID runs in Mathematica 10 [46] 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).

JHEP01(2020)024
PrepInt takes as input any master integral in figures 1, 2 and 3, as well as the threeloop 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 tables 4 and 5. For example  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.

JHEP01(2020)024
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  In this output, the numbers in front of pm[n] denote the errors of the five finite remainder functions that are visible in the output Out 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.

JHEP01(2020)024
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 UNIXtype 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.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.