A Numerical Routine for the Crossed Vertex Diagram with a Massive-Particle Loop

TL;DR

This paper introduces a numerical routine for efficiently evaluating the two master integrals of the crossed vertex diagram with a massive top-quark loop, enabling precise calculations across all phase space regions.

Contribution

It provides a new numerical method using differential equations and power series expansions for the crossed vertex diagram with massive particles, facilitating fast and accurate integral evaluations.

Findings

High-precision numerical evaluations achieved

Method applicable across all phase space regions

Routine simplifies complex integral calculations

Abstract

We present an evaluation of the two master integrals for the crossed vertex diagram with a closed loop of top quarks that allows for an easy numerical implementation. The differential equations obeyed by the master integrals are used to generate power series expansions centered around all the singular points. The different series are then matched numerically with high accuracy in intermediate points. The expansions allow a fast and precise numerical calculation of the two master integrals in all the regions of the phase space. A numerical routine that implements these expansions is presented.

Equations

$$x=-\frac{S}{m^2}=-s,$$

$$\int {\mathcal{D}}^d{k}_1{\mathcal{D}}^d{k}_2\frac{1}{D_1^{a_1}{D}_2^{a_2}{D}_3^{a_3}{D}_4^{a_4}{D}_5^{a_5}{D}_6^{a_6}{D}_7^{a_7}}.$$

$$D_{1,..7}$$

$${\mathcal{D}}^d{k}_i=\frac{d^d{k}_i}{4{\pi }^{\frac{d}{2}}\Gamma (1+ϵ)}{\left(\frac{m^2}{{\mu }^2}\right)}^{ϵ},$$

$$\frac{d{\mathcal{T}}_9}{dx}$$

$$\frac{d{\mathcal{T}}_{10}}{dx}$$

$${\mathcal{T}}_1$$

$$\frac{d^2{\mathcal{T}}_9}{d{x}^2}+p(x)\frac{d{\mathcal{T}}_9}{dx}+q(x){\mathcal{T}}_9=\Omega (x),$$

$$p(x)$$

$$q(x)$$

$$\Omega (x)$$

$$H(-r,-r,x)$$

$$H(r,0,x)$$

$${\mathcal{T}}_{10}=\frac{x}{4m^2}\frac{d{\mathcal{T}}_9}{dx}+\frac{1}{2m^2}{\mathcal{T}}_9.$$

$${\mathcal{T}}_9\sim \log x\text{mbox}forx\to 0.$$

$$\frac{d^2{\mathcal{T}}_9^{(0)}}{d{x}^2}+p(x)\frac{d{\mathcal{T}}_9^{(0)}}{dx}+q(x){\mathcal{T}}_9^{(0)}=0.$$

$$p(x)$$

$$q(x)$$

$${\mathcal{T}}_9^{(0)}(x)=x^{\alpha }\sum _{n=0}^{\infty }{a}_n{x}^n,$$

$${\left(\alpha +\frac{3}{2}\right)}^2=0,$$

$${\mathcal{T}}_9^{(0)}(x)=\frac{1}{\sqrt{x}}\sum _{n=-1}^{\infty }{a}_n{x}^n+\frac{\log x}{\sqrt{x}}\sum _{n=-1}^{\infty }{b}_n{x}^n,$$

$$a_0$$

$$a_1$$

$$a_2$$

$$\Omega (x)=\sum _{n=-2}^{\infty }{k}_n{x}^n+\log x\sum _{n=-2}^{\infty }{r}_n{x}^n,$$

$$k_{-2}$$

$$k_{-1}$$

$$k_0$$

$$\widetilde{{\mathcal{T}}_9}(x)=\sum _{n=-1}^{\infty }{p}_n{x}^n+\log x\sum _{n=-1}^{\infty }{q}_n{x}^n.$$

$$p_{-1}=0,q_{-1}=0,$$

$$p_0$$

$$p_1$$

$$p_2$$

$${\mathcal{T}}_9(x)=\frac{1}{\sqrt{x}}\sum _{n=-1}^{\infty }{a}_n{x}^n+\frac{\log x}{\sqrt{x}}\sum _{n=-1}^{\infty }{b}_n{x}^n+\sum _{n=0}^{\infty }{p}_n{x}^n+\log x\sum _{n=0}^{\infty }{q}_n{x}^n.$$

$$a_{-1}=0,b_{-1}=0,$$

$${\mathcal{T}}_9(x)=\sum _{n=0}^{\infty }{p}_n{x}^n+\log x\sum _{n=0}^{\infty }{q}_n{x}^n,$$

Full text

IFT-UAM/CSIC-18-118 HU-EP-18/37

**A Numerical Routine for the

Crossed Vertex Diagram with a Massive-Particle Loop

Roberto Bonciani*Email: roberto.bonciani@roma1.infn.it, *Dipartimento di Fisica, Sapienza - Università di Roma, 00185, Rome, Italy

and INFN Sezione di Roma, 00185, Rome, Italy* Giuseppe Degrassi†††Email: degrassi@fis.uniroma3.it, *Dipartimento di Matematica e Fisica, Università di Roma Tre, 00146 Rome, Italy

and INFN, Sezione di Roma Tre, 00146 Rome, Italy* Pier Paolo Giardino‡‡‡Email: pier.giardino@uam.es Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, 28049, Madrid, Spain Ramona Gröber§§§Email: ramona.groeber@physik.hu-berlin.de *Institut für Physik, Humboldt-Universität zu Berlin, 12489 Berlin, Germany

and Institute for Particle Physics Phenomenology, Department of Physics, Durham University, Durham, DH1 3LE, UK*

We present an evaluation of the two master integrals for the crossed vertex diagram with a closed loop of top quarks that allows for an easy numerical implementation. The differential equations obeyed by the master integrals are used to generate power series expansions centered around all the singular points. The different series are then matched numerically with high accuracy in intermediate points. The expansions allow a fast and precise numerical calculation of the two master integrals in all the regions of the phase space. A numerical routine that implements these expansions is presented.

Key words: Feynman diagrams, Multi-loop calculations

1 Introduction

In the last years, we witnessed an impressive progress in the analytic calculation of multi-loop Feynman diagrams. This progress was mainly due to a procedure which is by now standard and consists in the reduction of the dimensionally regularized scalar integrals to the Master Integrals (MIs) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], and their calculation using the differential equations [14, 15, 4, 16, 17].

With this procedure, it was possible to calculate massless quantum corrections to important processes in collider physics, that are now known to three and four loops [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. These corrections can be usually expressed in terms of generalized polylogarithms (GPLs) [28, 29, 30, 31]. While sometimes higher-order massive corrections can also be expressed in terms of GPLs [32, 33, 34, 35, 36], they reveal in general a more complicated structure. This is for instance the case of the two MIs of the equal-mass two-loop sunrise. The related system of first-order linear differential equations cannot be decoupled and it admits solutions in terms of complete elliptic integrals of the first and second kind [37, 38, 39, 40, 41, 42, 43]. This is also the case of three- [44, 45] and four-point functions [46, 47, 48] that were considered recently and whose solutions are expressed as iterated integrals over elliptic kernels multiplied by polylogarithmic terms. The study of these new functions has just started [49, 50, 51, 52, 53].

In this article, we focus on the calculation of the two MIs of the vertex crossed topology with a closed heavy-quark loop. These two MIs were studied in detail in Ref. [45], where the authors worked out completely their solution in terms of repeated integrations over elliptic kernels. They enter the calculation of several processes at the two-loop level in perturbation theory, as the production of top-antitop pairs in hadronic collisions [54, 55, 56, 57, 58, 47, 48, 59, 60], di-photon or di-jet production [61] and they are part of the coefficients of the $p_{T}$ expansion of the double Higgs production cross section, as discussed in Ref. [62].

Our goal is to present a Fortran numerical routine that can be easily used to evaluate the MIs for every real value of the dimensionless parameter $x=-S/m^{2}$, which the MIs depend on, with double precision. The approach we use is a semi-analytical approach to the solution of the differential equations, namely the expansion of the differential equation near singular points. It was proposed in Ref. [63] for the sunrise with three equal masses. In Ref. [44] the method was applied to a three-point function111See Refs.[64, 65] for two recent publications on the method. occurring in the calculation of the MIs that are involved in the two-loop corrections to the electroweak form factor [66, 67].

The paper is structured as follows. In Section 2, we discuss the MIs entering the 6-denominator topology of Fig. 1. We focus on the two crossed MIs (${\mathcal{T}}_{9},\>{\mathcal{T}}_{10}$) for which we present the relevant second-order linear differential equation that will be solved expanding the solution by series near the singular points. Section  3 is devoted to the discussion of the solution for ${\mathcal{T}}_{9}$ in the region $x\geq 0$. We present first the series in the two singular points $x=0$ and $x=16$ and their matching. Then, we discuss the expansion at infinity and how it can be matched to the expansion at $x=16$. In Section 4, we present the solution for ${\mathcal{T}}_{9}$ in the region $x<0$ obtained via the analytic continuation in the high-energy time-like region. In Section 5, we discuss the evaluation of the second master integral. Finally, in Section 6, we present the Fortran routine.

2 The Differential Equations for the two crossed

Master Integrals

We consider a process in which two massless particles with incoming momenta $p_{1}$ and $p_{2}$, such that $p_{1}^{2}=p_{2}^{2}=0$, annihilate into a particle with momentum $p=p_{1}+p_{2}$. We define the Mandelstam invariant $S=-(p_{1}+p_{2})^{2}$ and the dimensionless parameter

[TABLE]

where $s=S/m^{2}$ and $m$ is the mass of a massive state that runs into the loops.

The 6-denominator topology we are interested in relevant for this process is shown in Fig. 1. The dimensionally regularized scalar integrals belonging to that topology can be expressed in terms of

[TABLE]

In Eq. (2), $D_{i}$, $i=1,...,7$, are the denominators to which the following routing is assigned

[TABLE]

with $k_{1}$ and $k_{2}$ the loop momenta; $a_{i}$, with $i=1,...,7$, are integer numbers, $d=4-2\epsilon$ is the dimension of the space-time, and the normalization is such that222Note that we present, in the paper and in the routine, the euclidean version of the MIs, before Wick rotation.

[TABLE]

where $\mu$ is the scale of dimensional regularization.

The reduction to the MIs of the family in Eq. (2) are performed using the computer programs FIRE [7, 10, 11] and Reduze 2 [8, 9]. There are 10 MIs in total, shown in Fig. 2. All of them are known in the literature from previous works [68, 69, 67, 45].

We focus on the evaluation of ${\mathcal{T}}_{9}$ and ${\mathcal{T}}_{10}$ using the semi-analytic approach followed in Refs. [63, 70]. We concentrate on the system of first-order linear differential equations that involves the two coupled 6-denominator MIs ${\mathcal{T}}_{9}$ and ${\mathcal{T}}_{10}$. The two MIs are finite in $\epsilon$. Moreover, in all the processes mentioned in the introduction, at the NNLO, they enter in the calculation of the finite part of the cross section so that only the ${\mathcal{O}}(\epsilon^{0})$ is needed. At the ${\mathcal{O}}(\epsilon^{0})$, we find:

[TABLE]

where $\Omega_{2}(x)$ contains the MIs of the subtopologies and, at this order in $\epsilon$, is a function that can be expressed in terms of logarithms and dilogarithms of the variable $x$.

The system is equivalent to a single second-order linear differential equation for one of the two MIs involved. Let us consider ${\mathcal{T}}_{9}$. We find:

[TABLE]

with ($\Omega(x)=(4m^{2}/x)\Omega_{2}$)

[TABLE]

where we used the notation introduced in Refs. [67, 71] for the repeated integration over square roots

[TABLE]

The function $H(-r,-r,x)$ is real when $x\geq 0$. In the Minkowski region, $x\to-s-i0^{+}$, with $s>0$, $H(-r,-r,x)$ is real if $0<s<4$. For $s>4$, it develops an imaginary part due to the branch cut of the square root. The function $H(r,0,x)$ is real if $0<x<4$, while for $x>4$ the square root of the integrand has a branch cut. The result is purely imaginary and the sign depends on the sign of the small imaginary part that we add to $x$ to chose on which part of the cut we are. The same happens for the square root $\sqrt{x(4-x)}$ in Eq. (10). It is real for $0<x<4$ and purely imaginary for $x>4$. The combination $H(r,0,x)/\sqrt{x(4-x)}$ is real on the entire $x>0$ axis. Using consistently the same prescription for $H(r,0,x)$ and for $\sqrt{x(4-x)}$, we find that the ratio is real and independent on the prescription used.

Eq. (7) belongs to the Fuchs class, i.e. it has regular singular points only, eventually including the point at infinity. In our case, the singularities on the real axis are located at $x=0\,,\>x=16$, while also the point at infinity, $x=\infty$, is singular, as can be seen replacing the variable $x$ with $y=1/x$ and studying the equation in $y=0$.

The solution of the homogeneous equation associated to Eq. (7) can be expressed in terms of complete elliptic integrals of the first kind, and the particular solution is expressed as repeated integrations over the elliptic kernel, as it was discussed in detail in Ref. [45]. However, in this paper we are going to use another approach for the solution of the second-order differential equation. We will use the differential equation to generate power series expansions around the singular points and at infinity. Each series is determined up to two arbitrary constants. We will impose the constants of the series in $x=0$, since we know the initial conditions for ${\mathcal{T}}_{9}$ in that point. Then, the series are matched two-by-two in a point which lies inside both convergence domains. In this way, we will be able to fix all the constants and have a representation by series on the whole real axis. Our ultimate goal is to be able to evaluate precisely the function ${\mathcal{T}}_{9}$ on the whole real axis. In order to achieve the required precision it can be useful to supplement the original expansion in $x=0$, $x=16$ and infinity, with additional expansions around regular points.

Once the first master integral ${\mathcal{T}}_{9}$ has been determined, we can find the expression of the second, ${\mathcal{T}}_{10}$, using Eq. (5):

[TABLE]

3 ${\mathcal{T}}_{9}$ evaluation for $x\geq 0$

In this Section we discuss the solution of Eq. (7) in the region $x\geq 0$. ${\mathcal{T}}_{9}$ is obtained through the series in the singular regular points $x=0,\>x=16$ and $x=\infty$ that are then matched to cover the entire region $x\geq 0$. In all points, we first solve the homogeneous equation and then the complete equation, obtaining all the coefficients of the series in terms of the first two unknown coefficients. These unknowns will be fixed from the behaviour of the solution in one point, with the matching procedure.

3.1 The solution around $x=0$

The point $x=0$ allows us to impose the initial conditions and, therefore, to determine the two constants of integration that come from the general solution of the second-order linear differential equation (7). For this purpose, it is sufficient to know the behaviour of the master ${\mathcal{T}}_{9}$ for $x\to 0$ that can be obtained, for example, via a large-mass asymptotic expansion of the integral,

[TABLE]

This implies that in the solution no terms with inverse powers of $x$ appear, fixing the constants of integration.

We first consider the homogeneous equation

[TABLE]

The functions $p(x)$ and $q(x)$ have the following expansion in $x=0$:

[TABLE]

Since $x=0$ is a singular regular point, we look for a power series solution of the form:

[TABLE]

where $a_{n}$ are numerical coefficients determined from the differential equation itself and from the initial conditions. Substituting the solution (18) in Eq. (15), we obtain the characteristic equation for the determination of $\alpha$:

[TABLE]

with double solution in $\alpha=-3/2$. This implies two independent prefactors of the form $1/(x\sqrt{x})$ and $\log{x}/(x\sqrt{x})$. Therefore, the general solution of the homogeneous differential equation (15) is

[TABLE]

where we have absorbed a $1/x$ factor inside the series. The series (20) converges in a circle of radius $r=16$, i.e. up to the nearest divergence point on the real axis.

Expanding now the differential equation (15) in $x=0$ and substituting the general solution (20), we can fix all the coefficients of the series in terms of the first two coefficients, $a_{-1}$ and $b_{-1}$. The first few coefficients are:

[TABLE]

The general solution for ${\mathcal{T}}_{9}^{(0)}$ is a combination of two independent solutions, that can be found imposing, for instance, $a_{-1}=0$ and $b_{-1}=1$ or $a_{-1}=1$ and $b_{-1}=0$.

Let us now consider Eq. (7). The expansion of the function $\Omega(x)$ around $x=0$ is333In order to simplify the notation from now on we set $m^{2}=1$.

[TABLE]

with first coefficients

[TABLE]

Therefore, the inhomogenous term has a double pole in $x=0$, multiplied also by a single $\log{x}$. We look for a particular solution of Eq. (7) in $x=0$ of the form:

[TABLE]

Substituting Eq. (28) in the second-order differential equation expanded around $x=0$ we obtain, as in the case of the general solution of the homogeneous equation, terms $p_{n}$ and $q_{n}$ that depend on $p_{-1}$ and $q_{-1}$. However, in this case we are looking for a particular solution, since the general solution of the homogeneous equation is already known by Eq. (20). We can then choose to set

[TABLE]

finding the following first terms of the series in Eq. (28):

[TABLE]

The general solution of the complete equation is therefore:

[TABLE]

To determine completely the solution, we have to impose the initial conditions. Since ${\mathcal{T}}_{9}(x)$ can have at most a logarithmic singularity for $x\to 0$, the coefficients of the power singularities must vanish:

[TABLE]

and, as a consequence, all the $a_{n}$ and $b_{n}$ coefficients vanish.

Therefore, the solution of the complete equation reduces to

[TABLE]

where the first few coefficients $p_{n}$ and $q_{n}$ are given in Eqs. (30–32).

The solution given in Eq. (35) is real for $x>0$. However, in the physical region, $x<0$ ($s>0$), it develops an imaginary part that can be determined using the Feynman prescription $x\to-s-i0^{+}$. This means that the logarithmic terms develop an explicit imaginary part:

[TABLE]

Then, ${\mathcal{T}}_{9}(x)$ becomes complex for $x<0$ ($s>0$) with:

[TABLE]

3.2 The solution around $x\,=\,16$

The series in $x=0$ is completely determined. The following singular regular point we have to consider is $x=16$. Since the singular point closest to $x=16$ is $x=0$, the radius of convergence of the series in $x=16$ is $r=16$.

By solving the characteristic equation as in the previous subsection, we obtain a double zero in zero, so that the homogeneous equation has a solution of the form:

[TABLE]

The coefficients are, of course, different from the ones of the previous section, although we use the same notation to avoid introducing too many symbols. The first few coefficients read:

[TABLE]

The expansion of the inhomogenous term around $x=16$ is of the following form:

[TABLE]

where the first three coefficients $q_{n}$ are:

[TABLE]

In particular, note that there is no logarithmic term in Eq. (43).

The particular solution of the non-homogeneous equation in $x=16$ reads:

[TABLE]

The coefficients $r_{i}$ and $p_{i}$ depend on the $r_{0}$ and $p_{0}$, which are undetermined. However, since we are looking for a particular solution we can set from the beginning $r_{0}=0$ and $p_{0}=0$. This, in turn, forces the entire series of the logarithmic part of Eq. (47) to vanish, $p_{n}=0\quad\mbox{for all}\quad n=1,2,...$ . Therefore, we have a simple power series, with the first three terms given by:

[TABLE]

The general solution of the differential equation is given by

[TABLE]

Note that the integral ${\mathcal{T}}_{9}(x)$ should be real in the Euclidean region. However, the logarithmic terms, that come from the homogeneous solution, are responsible of the appearance of an imaginary part that cannot be there. We have, therefore, to impose that $b_{0}=0$. This condition implies that the logarithmic part of the expansion vanishes completely. The solution in Eq. (51) becomes a simple power series and depends on a single condition, $a_{0}$, that can be fixed as explained in the following section.

3.3 Matching the series in $x=0$ and $x=16$

The series expansion around $x=0$ is completely determined by imposing the initial conditions. The series in $x=16$, instead, depends on a single undetermined constant of integration, $a_{0}$. We can compute $a_{0}$, imposing that the series in $x=0$ and the one in $x=16$ assume the same value in a given point in the intersection of the respective domains of convergence. Since both series have radius of convergence $r=16$, in principle it is sufficient to impose that both series have the same value in any point $x\in(0,16)$.

Dealing with infinite series would exactly determine the coefficient $a_{0}$. However, we can only determine an arbitrary, but finite, number of coefficients of both series. Therefore, $a_{0}$ will be determined in an approximate way.

The number of terms in the series depend on the relative precision at which we want to be able to compute ${\mathcal{T}}_{9}(x)$ in a given point of the real axis. Our goal is to provide a double precision numerical routine, using a relatively small number of terms in the series (around 50 or less).

If we want to use just the series in $x=0$ and $x=16$ and we want to be able to provide such a precision, we have to deal with a large number of terms in the series. In order to keep the number of terms of the order of 50, and relative double precision within the interval $0<x<16$, we have to add series in intermediate points.

All the points in the interval $x\in(0,16)$ are regular points for the differential equation and they will result in simple power series (without the logarithmic part). In particular, we added series in $x=2$, 4, and 8. The procedure of matching is, therefore, performed as follows. We match the series in $x=0$ with the one in $x=2$. As a matching point we choose $x=1.5$. Then, in the point $x=3.25$ the series in $x=2$ is matched with the one in $x=4$, while in $x=6$, the series in $x=4$ is matched with the one in $x=8$ . Finally, the series in $x=8$ is matched with the one in $x=16$, in the point $x=12$.

The actual point in which we match two series is of course arbitrary. Nevertheless, a bad choice would lower the precision of the matched series. This would, in turn, lower the precision for all $x$ above the matching point. A possible approach for a good choice is the following. We first start with a point that assures a good convergence for both the series and we determine the unknown constants. Then, we vary a bit the point of the matching and we look at the corresponding variation of the significant digits of the constants. A good matching point maximises the number of stable digits in the result for the unknown constants.

3.4 The solution around $x\,=\,\infty$

We consider now the expansion of ${\mathcal{T}}_{9}$ around $x\,=\,\infty$. Since the closest singularity to $x=\infty$ is at $x=16$, we expect the expansion around infinity to be convergent outside the circle of radius 16, i.e. for $|x|>16$.

The expansion at infinity can be studied systematically by performing the following change of variable $x=1/y$ and, then, considering the limit $y\to 0$.

The homogeneous equation in $y\to 0$ limit reads

[TABLE]

We look for a solution of the form

[TABLE]

The characteristic equation gives $(\beta-2)^{2}=0$, with a double zero in $\beta=2$. Therefore, the solution of the homogeneous equation, in the original variable $x=1/y$ is

[TABLE]

with the coefficients $a_{n}$ and $b_{n}$ expressed in terms of the lowest-order ones, $a_{2}$ and $b_{2}$ as shown for the first few terms:

[TABLE]

The expansion of the non-homogeneous term $\Omega(x)$ around $x=\infty$ is of the form:

[TABLE]

where the lowest-order coefficients read:

[TABLE]

The differential equation involves a second derivative and the non-homogeneous term has double logarithmic terms. Therefore, the particular solution must contain up to four powers of the logarithm:

[TABLE]

Substituting Eq. (62) into the non-homogeneous equation, we obtain the following first few coefficients:

[TABLE]

Finally, the general solution is given by:

[TABLE]

where we set:

[TABLE]

The coefficients of the power series and of the single logarithm depend upon the two constants of integration, while the coefficients of the double, triple and fourth logarithm are uniquely determined. As in the case of the series in $x=16$, the two constants have to be determined matching the solution in $x=\infty$ with the one in $x=16$, in an intermediate point chosen in the range $16<x<32$ (the series in $x=16$ has radius of convergence $r=16$). However, in order to improve the precision in the determination of the integral ${\mathcal{T}}_{9}(x)$, without adding too many terms in the series expansions, it is better to add the expansions in three additional regular points: $x=32$, $x=64$ and $x=128$, before the matching with $x=\infty$.

4 ${\mathcal{T}}_{9}$ evaluation for $x<0$ ($s>0$)

The solution for ${\mathcal{T}}_{9}$ in the region $x<0$ can be constructed starting from the expansion of the amplitude ${\mathcal{T}}_{9}(x)$ for large time-like momenta, namely for $x\to-\infty$ ($s\to\infty$), that can be found from the asymptotic expansion in the space-like region ($x\to\infty$) by analytic continuation. With the Feynman prescription

[TABLE]

we have to consider that the logarithm develops an immaginary part as in Eq. (36):

[TABLE]

Then ${\mathcal{T}}_{9}(s)$ becomes complex and its real and immaginary parts are given by:

[TABLE]

The series in $x=0$ has a convergence radius $r=16$ and the series at infinity converges in $|x|>16$. In order to determine ${\mathcal{T}}_{9}$ in all points of the time-like region with the required precision, we need additional expansion points to sew the series at infinity with the one in $x=0$. Since in the region $-\infty<x<0$ ($0<s<\infty$) there are no singular points, the points to be added will be regular points, and the corresponding series will be simple power series.

We added the following points: $s=4$, $s=8$, $s=16$, $s=32$, $s=64$ and finally $s=128$. We will discuss extensively just $s=16$.

4.1 The solution around $s=16$

The point $s=16$ is a regular point. Therefore, the expansion of the homogeneous solution is a power series

[TABLE]

with the first few coefficients given in terms of $a_{0}$ and $a_{1}$ by

[TABLE]

The expansion of the inhomogneous term around $s=16$ is of the following form

[TABLE]

where

[TABLE]

Therefore, the particular solution of the differential equation in $s=16$ is, again, a power series in which the coefficients $p_{n}$, $n\geq 2$, depend upon the first two coefficients, $p_{0}$ and $p_{1}$. Since we are now looking for a particular solution, we can set $p_{0}=0$ and $p_{1}=0$, finding

[TABLE]

with the first few coefficients that read

[TABLE]

Finally, the general solution is given by

[TABLE]

where

[TABLE]

5 Expansions for the Master Integral ${\mathcal{T}}_{10}$

The second MI is directly determined from the first one by means of Eq. (13):

[TABLE]

Knowing the series expressions for ${\mathcal{T}}_{9}$, Eq. (88) allows to determine ${\mathcal{T}}_{10}$ performing a simple derivative.

The matching conditions that we imposed for the series expansions in the various points of the real axis, for the complete determination of ${\mathcal{T}}_{9}$, are still valid for ${\mathcal{T}}_{10}$. In principle, one simply has to take the derivative of each of these series, the series itself, and combine them in order to fulfill Eq. (88). In the case of infinite series, there would be no difference in the determination and precision of ${\mathcal{T}}_{10}$ with respect to what we found for ${\mathcal{T}}_{9}$. However, we deal with truncated series and this means that the optimal choice for a matching point of two series for ${\mathcal{T}}_{9}$ can be less optimal for the corresponding series of ${\mathcal{T}}_{10}$. Therefore, we decided to determine the matching points independently for the series of ${\mathcal{T}}_{9}$ and ${\mathcal{T}}_{10}$.

We used the criterium of maximization of the number of stable digits in the determination of the unknown constants, under the variation of the matching point. In so doing, we found that the matching points for corresponding pair of series of ${\mathcal{T}}_{9}$ and ${\mathcal{T}}_{10}$ give rise to slightly different matching constants. We used the difference between the values of the matching constants as an indicator for the precision at which we can claim the series reproduce the numerical value of the masters. In all the matching points, we found corresponding matching constants that agree with double precision.

6 The Fortran Routine

In this section we give details on the numerical routines that accompany the paper.

The routine implements the series in the various points of the real domain discussed in the previous sections. In some points (in particular in $x=2,4,8,16,32,64,128$ and $s=4,8,12,16,32,64,128$), in order to improve the convergence of the series, we performed the Bernoulli variable transformation [72, 73], which is defined as

[TABLE]

for a series expansion around $x_{0}$, with nearest singular points $a$ and $b$, such that $a<x_{0}<b$. This change of variable usually increases the convergence of the series near the point of expansion (see for instance Refs. [63, 44]). Althought in the points indicated above we found a considerable increase in such convergence, resulting in an increase of the number of reliable digits of the final result, in $x=0$ and $x=\pm\infty$ the original power series worked at the same level of accuracy (or sometimes even better) or had better numerical behaviour. Therefore, the routines are written using the original series in $x=0$ and $x=\pm\infty$, and the series in the Bernoulli variable in all the regular points and $x=16$.

The numerical program consists of the header file main_elliptic.f and the two main files MI1.f and MI2.f, which compute the master integrals $\mathcal{T}_{9}$ and $\mathcal{T}_{10}$ respectively. The program is written in FORTRAN. Several files contain the lengthy formulae of the expansions around the various points.

The program can be used in the two following ways:

• •

Way 1: As a whole with output onto the screen and into an outputfile. After unzipping the files

[TABLE]

the program can then be compiled with the provided makefile, meaning by typing

[TABLE]

and run by the command

[TABLE]

If no input value for $x=-s$ is given, the program interrupts and asks to input a value. If instead no input for the name of the output file is given, the output is written into a default file named output_MI.dat.

• •

Way 2: Inside another program. In this case only the files MI1.f for $\mathcal{T}_{9}$ and MI2.f for $\mathcal{T}_{10}$ are needed as well as all files in the folder seriesexpansions. The makefile of the other program must then be adjusted by adding MI1.o or MI2.o to the files to be compiled. The function complex16 MI1(double precision x) for $\mathcal{T}_{9}$ or complex16 MI2(double precision x) for $\mathcal{T}_{10}$ can then be called directly within any other FORTRAN program.

In the following we list the various files and explain them in more details.

• •

main_elliptic.f: The main file calls the functions MI1(x) and MI2(x) for the value $x$ given by the user as an argument when running the program and writes the output. This file is not needed if the user wants to call the integrals from within his/her own program.

• •

MI1.f: Computes the master integral $\mathcal{T}_{9}$ and can be called by the user directly if he/she decides to call the integral from within his/her own program. MI1.f decides which series expansion is needed for the given value of $x$ and returns the corresponding value. It needs the help files that are provided in the folder seriesexpansions.

• •

MI2.f: Same as MI1.f but for $\mathcal{T}_{10}$.

• •

seriesexpansions/MI1_in_x_#.f: Help files that contain the lengthy expressions for the series expansions of $\mathcal{T}_{9}$ around $x=0,2,4,8,16,32,64,128$ and $\infty$, where # stands for the respective value of $x$.

• •

seriesexpansions/MI1_in_s_#.f: Help file that contains the lengthy expressions for the series expansions of $\mathcal{T}_{9}$ around $s=4,8,12,16,32,64,128$ and $\infty$.

• •

seriesexpansions/MI2_in_x_#.f and seriesexpansions/MI2_in_s_#.f: Same as seriesexpansions/MI1_in_x_#.f and seriesexpansions/MI1_in_s_#.f respectively but for $\mathcal{T}_{10}$.

6.1 Numerical Checks

We performed several numerical checks both to ensure the correctness of our results and to check the numerical accuracy. We will describe them in the following.

• •

As outlined in Section 3.3, we determined the constants of integration by matching the series in points within the radius of convergence of two series, starting from $x=0$ which we have determined completely, going to $x=\infty$ and $x=-\infty$. We did this procedure for both $\mathcal{T}_{9}$ and $\mathcal{T}_{10}$ separately. If we would be able to expand the series to arbitrary high order, the constants of integration of $\mathcal{T}_{9}$ and $\mathcal{T}_{10}$ would be the same. However, we work with truncated series, and the determination of the constants depend upon the details of the series used, as for instance the form of the coefficients and the number of terms. As a consequence, the integration constants are not exactly the same. This allows us to use the comparison of the matching constants between $\mathcal{T}_{9}$ and $\mathcal{T}_{10}$ for each series to determine internally the numerical accuracy of the procedure. Doing so we find agreement in all the series to double precision accuracy.

• •

As an internal check and as a determination of the accuracy of the result, we adopted the following strategy. The series in $x=0$ is completely determined, since we impose the initial conditions in that point. Starting from $x=0$, we match the undetermined constants of the series as described in the paper up to the series in $x=\infty$. Independently, we perform the same procedure in the Minkowski region, starting from $s=0$ up to the series in $s=\infty$. Now with the series in $x=\infty$, we perform an analytic continuation to the Minkowski region, $s>0$, and we numerically evaluate the series in $s=1000$, comparing the result with the numerical evaluation, in the same point, of the series obtained with the matchings in the Minkowski region. We find that the two numbers agree with double precision.

• •

We cross-checked our numerical routines against PySecDec in several points of the entire domain, in both the Euclidean and Minkowski regions. We found complete agreement within the numerical accuracy of PySecDec, which is limited to 5-6 digits.

• •

The most stringent test was the one done with the numbers coming from the exact solution of Ref. [45]. We could check our routines against the numbers provided by the authors of Ref. [45] in $x=3,13,50$ and $s=3,5,18,50$, finding an agreement to double precision accuracy444In the comparison between the numbers of our routines and the numbers of Ref. [45] it must be remembered that the normalization of the integrals are different in the two works, as can be seen from our Eq. (4) and Eq. (2.2) of Ref.[45]. In particular, in order to match the numbers of Ref. [45], our $\mathcal{T}_{9}$ has to be multiplied by $16$, while $\mathcal{T}_{10}$ for $-16$..

7 Conclusions

In this paper we presented a semi-analytical evaluation for the two MIs of the crossed vertex topology with a closed massive loop, implemented in a Fortran numerical routine.

The two MIs can be expressed in power series of the dimensional parameter $\epsilon=(4-d)/2$. Each order in $\epsilon$ fulfills a system of two coupled first-order linear differential equations, that admits solutions in terms of one-fold integrals of complete elliptic integrals of the first and second kind times polylogarithmic terms (see Ref. [45]). In the present paper we focus on the solution of the differential equations for the ${\mathcal{O}}(\epsilon^{0})$, which is relevant for phenomenological applications at the NNLO.

In order to implement the solutions in a Fortran numerical routine, for the precise evaluation of the two MIs, we followed a standard approach that was used in the past for the study of the equal-mass sunrise and the three-point function with two massive exchanges, namely the solution by series of the equivalent second-order linear differential equation for one of the masters. The other master is then calculated by a simple derivative, once the first master is known.

Expanding the master in series in the singular points of the differential equation we were able to directly construct a solution that covers the entire range $-\infty\leq x\leq\infty$ which is suitable for precise numerical evaluations.

The Fortran routine presented in this work returns the numerical value of the two MIs for every real value of the dimensionless parameter they depend on, with double precision accuracy.

Acknowledgments

We wish to thank Andreas von Manteuffel and Lorenzo Tancredi for having provided us with the numeric results of the two master integrals published in Ref. [45] in additional points of the phase spase, for cross-checks. P.P.G. acknowledges the support of the Spanish Agencia Estatal de Investigacion through the grant “IFT Centro de Excelencia Severo Ochoa SEV-2016-0597”. R.G. is partially supported by a European Union COFUND/Durham Junior Research Fellowship under the EU grant number 609412. R.G. thanks the Galileo Galilei Institute for Theoretical Physics for the hospitality and the INFN for partial support during the completion of this work