TSIL: a program for the calculation of two-loop self-energy integrals

Abstract

TSIL is a library of utilities for the numerical calculation of dimensionally regularized two-loop self-energy integrals. A convenient basis for these functions is given by the integrals obtained at the end of O.V. Tarasov's recurrence relation algorithm. The program computes the values of all of these basis functions, for arbitrary input masses and external momentum. When analytical expressions in terms of polylogarithms are available, they are used. Otherwise, the evaluation proceeds by a Runge–Kutta integration of the coupled first-order differential equations for the basis integrals, using the external momentum invariant as the independent variable. The starting point of the integration is provided by known analytic expressions at (or near) zero external momentum. The code is written in C, and may be linked from C/C++ or Fortran. A Fortran interface is provided. We describe the structure and usage of the program, and provide a simple example application. We also compute two new cases analytically, and compare all of our notations and conventions for the two-loop self-energy integrals to those used by several other groups.

Program summary

Title of program: TSIL

Version number: 1.0

Catalogue identifier: ADWS

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/ADWS

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland

Programming language: C

Platform: Any platform supporting the GNU Compiler Collection (gcc), the Intel C compiler (icc), or a similar C compiler with support for complex mathematics

No. of lines in distributed program, including test data, etc.: 42 730

No. of bytes in distributed program, including test data, etc.: 297 101

Distribution format: tar.gz

Nature of physical problem: Numerical evaluation of dimensionally regulated Feynman integrals needed in two-loop self-energy calculations in relativistic quantum field theory in four dimensions.

Method of solution: Analytical evaluation in terms of polylogarithms when possible, otherwise through Runge–Kutta solution of differential equations.

Limitations: Loss of accuracy in some unnatural threshold cases that do not have vanishing masses.

Typical running time: Less than a second.

Introduction

The precision of measurements within the Standard Model requires a level of accuracy obtained from two-loop, and even higher-order, calculations in quantum field theory. The future discoveries of the Large Hadron Collider and a future $e^{+}{e}^{-}$ linear collider may well extend these requirements even further. For example, supersymmetry predicts the existence of many new particle states with perturbative interactions. The most important observables in softly broken supersymmetry are just the new particle masses. Therefore, comparisons of particular models of supersymmetry breaking with experiment will require systematic methods for two-loop computations of physical pole masses in terms of the underlying Lagrangian parameters.

In this paper, we describe a software package called TSIL (Two-loop Self-energy Integral Library)¹ that can be used to compute two-loop self-energy functions. In a general quantum field theory, the necessary dimensionally regularized Feynman integrals can be expressed in the form:

$${\mu }^{8-2d}\int {\mathrm{d}}^{d}k\int {\mathrm{d}}^{d}q\frac{{(k^2)}^{n_1}{(q^2)}^{n_2}{(k\cdot p)}^{n_3}{(q\cdot p)}^{n_4}{(k\cdot q)}^{n_5}}{{[k^2+x]}^{r_1}{[q^2+y]}^{r_2}{[{(k-p)}^2+z]}^{r_3}{[{(q-p)}^2+u]}^{r_4}{[{(k-q)}^2+v]}^{r_5}},$$

for non-negative integers $n_i,r_i$. (See the master topology M in Fig. 1.) The integrations are over Euclidean momenta in

$$d=4-2\mathrm{\epsilon }$$

dimensions. An algorithm has been derived by Tarasov [1], and implemented in a Mathematica program TARCER by Mertig and Scharf [2], that allows all such integrals to be systematically reduced to linear combinations of a few basis integrals. The coefficients are ratios of polynomials in the external momentum invariant and the propagator squared masses $x,y,z,u,v$. In the remainder of this section, we will describe our notations and conventions² for the two-loop basis integrals and some related functions, and describe the strategy used by TSIL to compute them.

First, we define a loop factor

$$C=(16{\pi }^2)\frac{{\mu }^{2\mathrm{\epsilon }}}{{(2\pi )}^d}={(2\pi \mu )}^{2\mathrm{\epsilon }}/{\pi }^2.$$

The regularization scale μ is related to the renormalization scale Q (in the modified minimal subtraction renormalization scheme based on dimensional regularization [6], or dimensional reduction [7] for softly-broken supersymmetric theories) by

$$Q^2=4\pi e^{-\gamma }{\mu }^2.$$

Logarithms of dimensionful quantities are always given in terms of

$$\overline{\ln }X\equiv \ln (X/Q^2).$$

The loop integrals are functions of a common external momentum invariant

$$s=-p^2$$

using a Euclidean or signature ($-+++$) metric. (Note that the sign convention is such that for a stable physical particle with mass m, there is a pole at $s=m^2$.) On the physical sheet, s has an infinitesimal positive imaginary part. Since all functions in any given equation typically have the same s and $Q^2$, they are not included explicitly in the list of arguments in written formulas. A prime on a squared-mass argument of a function indicates differentiation with respect to that argument. Thus, for example,

$$f(x^{\prime },x,z^{\prime })={[\frac{{\partial }^2}{\partial y\partial z}f(y,x,z)]}_{y=x}.$$

We now define one-loop vacuum and self-energy integrals as:

$$\mathbf{A}(x)=C\int {\mathrm{d}}^{d}k\frac{1}{[k^2+x]},$$

$$\mathbf{B}(x,y)=C\int {\mathrm{d}}^{d}k\frac{1}{[k^2+x][{(k-p)}^2+y]}.$$

We also define two-loop integrals according to:

$$\mathbf{S}(x,y,z)=C^2\int {\mathrm{d}}^{d}k\int {\mathrm{d}}^{d}q\frac{1}{[k^2+x][q^2+y][{(k+q-p)}^2+z]},$$

$$\mathbf{I}(x,y,z)=\mathbf{S}(x,y,z){|}_{s=0},$$

$$\mathbf{T}(x,y,z)=-\mathbf{S}(x^{\prime },y,z)$$

$$\mathbf{U}(x,y,z,u)=C^2\int {\mathrm{d}}^{d}k\int {\mathrm{d}}^{d}q\frac{1}{[k^2+x][{(k-p)}^2+y][q^2+z][{(q+k-p)}^2+u]},$$

$$\mathbf{V}(x,y,z,u)=-\mathbf{U}(x,y^{\prime },z,u),$$

$$\mathbf{M}(x,y,z,u,v)=C^2\int {\mathrm{d}}^{d}k\int {\mathrm{d}}^{d}q\frac{1}{[k^2+x][q^2+y][{(k-p)}^2+z][{(q-p)}^2+u][{(k-q)}^2+v]}.$$

The corresponding Feynman diagram topologies are shown in Fig. 1. These integrals have various symmetries that are clear from the diagrams: $\mathbf{B}(x,y)$ is invariant under interchange of $x,y$; the “sunrise” integral $\mathbf{S}(x,y,z)$ and $\mathbf{I}(x,y,z)$ are invariant under interchange of any two of $x,y,z$; $\mathbf{T}(x,y,z)$ is invariant under $y\leftrightarrow z$; $\mathbf{U}(x,y,z,u)$ and $\mathbf{V}(x,y,z,u)$ are invariant under $z\leftrightarrow u$; and the “master” integral $\mathbf{M}(x,y,z,u,v)$ is invariant under each of the interchanges $(x,z)\leftrightarrow (y,u)$, and $(x,y)\leftrightarrow (z,u)$, and $(x,y)\leftrightarrow (u,z)$.

It is often convenient to introduce modified integrals in which appropriate divergent parts have been subtracted and the regulator removed. At one-loop order, we define the finite and ε-independent integrals:

$$A(x)=\underset{\mathrm{\epsilon }\to 0}{\lim }[\mathbf{A}(x)+x/\mathrm{\epsilon }]=x(\overline{\ln }x-1),$$

$$B(x,y)=\underset{\mathrm{\epsilon }\to 0}{\lim }[\mathbf{B}(x,y)-1/\mathrm{\epsilon }]=-\underset{0}{\overset{1}{\int }}\mathrm{d}t\overline{\ln }[tx+(1-t)y-t(1-t)s].$$

At two loops, we let

$$S(x,y,z)=\underset{\mathrm{\epsilon }\to 0}{\lim }[\mathbf{S}(x,y,z)-S_{\mathrm{div}}^{(1)}(x,y,z)-S_{\mathrm{div}}^{(2)}(x,y,z)],$$

where

$$S_{\mathrm{div}}^{(1)}(x,y,z)=(\mathbf{A}(x)+\mathbf{A}(y)+\mathbf{A}(z))/\mathrm{\epsilon },$$

$$S_{\mathrm{div}}^{(2)}(x,y,z)=(x+y+z)/2{\mathrm{\epsilon }}^2+(s/2-x-y-z)/2\mathrm{\epsilon }$$

are the contributions from one-loop subdivergences and from the remaining two-loop divergences, respectively. In addition,

$$I(x,y,z)=S(x,y,z){|}_{s=0},$$

$$T(x,y,z)=-S(x^{\prime },y,z).$$

Similarly, we define

$$U(x,y,z,u)=\underset{\mathrm{\epsilon }\to 0}{\lim }[\mathbf{U}(x,y,z,u)-U_{\mathrm{div}}^{(1)}(x,y)-U_{\mathrm{div}}^{(2)}],$$

where

$$U_{\mathrm{div}}^{(1)}(x,y)=\mathbf{B}(x,y)/\mathrm{\epsilon },$$

$$U_{\mathrm{div}}^{(2)}=-1/2{\mathrm{\epsilon }}^2+1/2\mathrm{\epsilon }$$

are again the contributions from one-loop sub-divergences and the remaining two-loop divergences. Also, we define

$$V(x,y,z,u)=-U(x,y^{\prime },z,u).$$

The master integral is free of divergences, so we define

$$M(x,y,z,u,v)=\underset{\mathrm{\epsilon }\to 0}{\lim }\mathbf{M}(x,y,z,u,v).$$

Thus, bold-faced letters A, B, I, S, T, U, V represent the original regularized integrals that diverge as $\mathrm{\epsilon }\to 0$, while the ordinary letters A, B, I, S, T, U, V, M are used to represent functions that are finite and independent of ε by definition. Note, however, that as defined above I, S, T, U, V are not simply the ε-independent terms in expansions in small ε. The following expansions are useful for converting between I, S, T, U, V and I, S, T, U, V:

$$\mathbf{A}(x)=-x/\mathrm{\epsilon }+A(x)+\mathrm{\epsilon }{A}_{\mathrm{\epsilon }}(x)+\mathcal{O}({\mathrm{\epsilon }}^2),$$

$$\mathbf{B}(x,y)=1/\mathrm{\epsilon }+B(x,y)+\mathrm{\epsilon }{B}_{\mathrm{\epsilon }}(x,y)+\mathcal{O}({\mathrm{\epsilon }}^2),$$

where

$$A_{\mathrm{\epsilon }}(x)=x[-1-\zeta (2)/2+\overline{\ln }x-{(\overline{\ln }x)}^2/2],$$

$$B_{\mathrm{\epsilon }}(x,y)=\zeta (2)/2+\frac{1}{2}\underset{0}{\overset{1}{\int }}\mathrm{d}t{(\overline{\ln }[tx+(1-t)y-t(1-t)s])}^2.$$

Here ζ is the Riemann zeta function. The function $B_{\mathrm{\epsilon }}(x,y)$ can be expressed in terms of dilogarithms [8], and is given by the coefficient of δ in Eq. (83) of Ref. [9]. Also,

$$\mathbf{B}(x^{\prime },y)=[(3-d)(s-x+y)\mathbf{B}(x,y)+(2-d)\{\mathbf{A}(y)+(s-x-y)\mathbf{A}(x)/2x\}]/{\mathrm{\Delta }}_{sxy},$$

where

$${\mathrm{\Delta }}_{abc}\equiv a^2+b^2+c^2-2ab-2ac-2bc.$$

From the preceding equations it follows that

$$\mathbf{I}(x,y,z)=-(x+y+z)/2{\mathrm{\epsilon }}^2+[A(x)+A(y)+A(z)-(x+y+z)/2]/\mathrm{\epsilon }+I(x,y,z)+A_{\mathrm{\epsilon }}(x)+A_{\mathrm{\epsilon }}(y)+A_{\mathrm{\epsilon }}(z)+\mathcal{O}(\mathrm{\epsilon }),$$

$$\mathbf{S}(x,y,z)=-(x+y+z)/2{\mathrm{\epsilon }}^2+[A(x)+A(y)+A(z)-(x+y+z)/2+s/4]/\mathrm{\epsilon }+S(x,y,z)+A_{\mathrm{\epsilon }}(x)+A_{\mathrm{\epsilon }}(y)+A_{\mathrm{\epsilon }}(z)+\mathcal{O}(\mathrm{\epsilon }),$$

$$\mathbf{T}(x,y,z)=1/2{\mathrm{\epsilon }}^2-[A(x)/x+1/2]/\mathrm{\epsilon }+T(x,y,z)+[A(x)-A_{\mathrm{\epsilon }}(x)]/x+\mathcal{O}(\mathrm{\epsilon }),$$

$$\mathbf{U}(x,y,z,u)=1/2{\mathrm{\epsilon }}^2+[1/2+B(x,y)]/\mathrm{\epsilon }+U(x,y,z,u)+B_{\mathrm{\epsilon }}(x,y)+\mathcal{O}(\mathrm{\epsilon }),$$

$$\mathbf{V}(x,y,z,u)=\frac{1}{\mathrm{\epsilon }}[(s+x-y)(B(x,y)-1)+2A(x)+(s-x-y)A(y)/y]/{\mathrm{\Delta }}_{sxy}+V(x,y,z,u)+\{(s+x-y)[B_{\mathrm{\epsilon }}(x,y)-2B(x,y)]+2A_{\mathrm{\epsilon }}(x)-2A(x)+(s-x-y)[A_{\mathrm{\epsilon }}(y)-A(y)]/y\}/{\mathrm{\Delta }}_{sxy}+\mathcal{O}(\mathrm{\epsilon }),$$

$$\mathbf{M}(x,y,z,u,v)=M(x,y,z,u,v)+\mathcal{O}(\mathrm{\epsilon }).$$

The internal workings of the TSIL code use the functions $A,B,I,S,T,U,V,M$ rather than their bold-faced counterparts. This is because we find that renormalized expressions for physical quantities are more compactly written in terms of the non-bold-faced integrals. However, both types of functions are available as outputs, with the proviso that for I, S, T, U, V, M we keep only the pole and finite terms as $\mathrm{\epsilon }\to 0$, and for A, B only terms through order ε.

Tarasov's algorithm [1] allows any integral of the form of Eq. (1.1) to be expressed³ as a linear combination of the following two-loop basis integrals:

$$M(x,y,z,u,v),U(z,x,y,v),U(u,y,x,v),U(x,z,u,v),U(y,u,z,v),T(v,y,z),T(u,x,v),T(y,z,v),T(x,u,v),T(z,y,v),T(v,x,u),S(v,y,z),S(u,x,v),$$

plus terms involving the two-loop vacuum integrals $I(x,y,v)$ or $I(z,u,v)$, or quadratic in the one-loop integrals.

In particular, the V and V integrals can be expressed as linear combinations of the others (see Eqs. (3.22)–(3.28) and (6.21) of Ref. [4]), and so are not included in the basis. However, they are included as outputs, because some results are more compactly written in terms of them. Because $T(x,y,z)$ is divergent in the limit $x\to 0$, it is also sometimes useful to define the function:

$$\overline{T}(x,y,z)=T(x,y,z)+B(y,z)\overline{\ln }x.$$

For $x=0$, this function is well-defined and can be written in terms of dilogarithms (see Eqs. (6.18)–(6.19) of Ref. [4]). In that limit it can also be rewritten in terms of the other basis functions, see Eqs. (A.15)–(A.16) of Ref. [5], but is still available as an output of TSIL for convenience.

It remains to provide a means for the numerical computation of the basis integrals. For special values of masses and external momentum, it is possible to compute the two-loop integrals analytically in terms of polylogarithms [10]. This requires [11] that there is no three-particle cut of the diagram for which the cut propagator masses $m_1$, $m_2$, $m_3$ and the five quantities

$$s_{\mathrm{cut}},s_{\mathrm{cut}}-{(m_1\pm m_2\pm m_3)}^2,$$

(where $s_{\mathrm{cut}}=-p_{\mathrm{cut}}^2$ is the momentum invariant for the total momentum flowing across the cut) are all non-zero. Many analytical results for various such special cases have been worked out in Refs. [4], [9], [12], [13], [14], [15], [16], [17], [18], [19], [20], and Appendix B of the present paper. There are also expansions [21], [22], [23], [24], [25] in large and small values of the external momentum invariant, and near the thresholds and pseudo-thresholds [26], [27], [28], [29], [30], [31], [32], [33], [34], [35]. Analytical results in terms of polylogarithms for the A, B, I, S, T, U, M functions at generic values of s are reviewed in Section VI of Ref. [4]. These and some other analytical formulas for special cases have been implemented in the TSIL code. They consist of the master integral cases:

$$\begin{array}{cccc}M(x,x,y,y,0),\hfill & M(0,0,0,x,0),\hfill & M(0,x,0,x,x),\hfill & \text{Ref. [16]}\hfill \\ M(0,0,0,0,x),\hfill & M(0,x,0,y,0),\hfill & \hfill & \text{Ref. [9]}\hfill \\ M(0,0,0,x,x),\hfill & M(0,x,x,0,0),\hfill & M(0,x,x,x,0),\hfill & \text{Ref. [19]}\hfill \\ M(0,0,x,y,0),\hfill & \hfill & \hfill & \text{Appendix B}\hfill \\ M(0,x,x,0,x){|}_{s=x},\hfill & \hfill & \hfill & \text{Ref. [16]}\hfill \\ M(0,y,y,0,x){|}_{s=x},\hfill & \hfill & \hfill & \text{Appendix B}\hfill \\ M(0,x,y,0,y){|}_{s=x},\hfill & \hfill & \hfill & \text{Refs. [36], [37] (see also Appendix B)}\hfill \end{array}$$

and functions obtained by permutations of the arguments, and all subordinate integrals A, B, I, S, T, U, V, $\overline{T}$ that have topologies obtained by removing one or more propagator lines from the cases above. These include:

$$\begin{array}{cc}U(x,y,0,y),\hfill & \text{Ref. [15]}\hfill \\ U(x,y,0,0),U(0,0,0,x),\hfill & \text{Ref. [9]}\hfill \\ U(0,x,y,z),\hfill & \text{Ref. [4]}\hfill \\ U(x,0,0,y),\hfill & \text{Ref. [5]}\hfill \\ U(x,0,y,y){|}_{s=x},U(y,0,y,x){|}_{s=x},\hfill & \text{Refs. [38], [39] (see also Appendix B)}\hfill \end{array}$$

and

$$\begin{array}{cc}S(0,x,y),T(x,0,y),\overline{T}(0,x,y),\hfill & \text{Ref. [18]}\hfill \\ S(x,y,y){|}_{s=x},T(x,y,y){|}_{s=x},T(y,x,y){|}_{s=x}.\hfill & \text{Refs. [29], [38]}\hfill \end{array}$$

Also included are all of the functions at $s=0$, which can be easily expressed in terms of the A and I functions and derivatives of them, which are in turn known [14], [40] analytically in terms of logarithms and dilogarithms.

For the case of generic masses and s, another method is needed. Integral representations have been studied in Refs. [41], [42], [43], [44], [45], [46], [47], [48], [49], [50]. For TSIL, we instead use the differential equations method [51], [52] to evaluate the integrals numerically. For the S, T and U integrals, this strategy was proposed and implemented in [53], [54], [55], [56], [57]. The method was rewritten in terms of the S, T, U integrals and extended to M in Ref. [4]. To see how the method works, consider the functions listed in Eq. (1.40) and also $B(x,z)$ and $B(y,u)$ and the product $B(x,z)B(y,u)$. Let us denote these sixteen quantities by $F_i$ where $i=1,\dots ,16$. Considered as functions of s for fixed x, y, z, u, v, $Q^2$ they can be shown to satisfy a set of coupled first-order differential equations of the form

$$\frac{\mathrm{d}}{\mathrm{d}s}{F}_i=\sum _j{C}_{ij}{F}_j+C_i,$$

where the coefficients $C_{ij}$ and $C_i$ are ratios of polynomials in the squared masses and s (and of the analytically known A and I functions, for some of the coefficients not multiplying two-loop functions). The coefficients can be evaluated by using Tarasov's recurrence relations, and were listed in [4].

At $s=0$ the values of all of the functions and their derivatives with respect to s (and/or their expansions in small s) are known analytically in terms of dilogarithms. Therefore, one can integrate the functions⁴ by the Runge–Kutta method to any other value of s. In order to avoid numerical problems from integrating through thresholds and pseudo-thresholds, we use the suggestion of Ref. [56] by following a displaced contour in the upper-half complex s plane, as shown in Fig. 2, whenever s is greater than the smallest non-zero threshold or pseudo-threshold. This contour starts from $s=0$ (or, in some cases, a point very close by, as explained in the next section) and proceeds to a point $i{s}_{\mathrm{im}}$, from there to $s+i{s}_{\mathrm{im}}$, and from there to the desired value s. Here $s_{\mathrm{im}}$ is real and positive. Since s has an infinitesimal positive imaginary part on the physical sheet, this procedure also automatically produces the correct branch cut behavior for functions when s is above thresholds. When s is less than or equal to the smallest non-zero threshold or pseudo-threshold, the Runge–Kutta integration instead proceeds directly along the real s-axis. In typical physical problems, if the master integral is needed, then so will be all of its subordinate $B,S,T,U$ integrals. These are all obtained simultaneously as a result of the Runge–Kutta method. Furthermore, checks on the numerical accuracy can be made by varying the Runge–Kutta step size parameters and the choice of contour in the upper-half complex s plane.

Most of the practical difficulties in realizing this program have to do with numerical instabilities when the final value of s is at or near a threshold or pseudo-threshold, or when the starting point $s=0$ is itself a threshold. We describe these issues and the methods used by TSIL to successfully evade them in the next section.