3-loop Feynman Integral Extrapolations for the Baseball Diagram

Higher order corrections are required for accurate theoretical prediction improvements in the technology of high energy physics experiments. Feynman loop integrals arise in the calculations of the Feynman diagrammatic approach, which is commonly used to address higher order corrections. Loop integrals may suffer from integrand singularities or irregularities at the boundaries and/or in the interior of the integration domain.

We recently explored methods for higher-order corrections to 2-loop Feynman integrals in the Euclidean or physical kinematical region [1], using numerical extrapolation and adaptive iterated integration. Our current goal is to address a 3-loop two-point integral for the “baseball” diagram of Figure 1 with six internal lines.

A representation of an $L$-loop Feynman integral with $N$ internal lines is given by (e.g., [2, 3, 4, 5])

		$\displaystyle{\mathcal{F}}={\Gamma(N-\frac{\nu L}{2})}\,(-1)^{N}\int_{{\mathcal{C}}_{N}}\prod_{r=1}^{N}dx_{r}\,\delta(1-\sum x_{r})\,U^{-\nu/2}(V-i\varrho)^{\nu L/2-N}$		(1)
		$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{or,\leavevmode\nobreak\ for\leavevmode\nobreak\ }L=3,N=6:$
		$\displaystyle{\mathcal{F}}={\Gamma(3\varepsilon)}\,\int_{{\mathcal{C}}_{6}}\prod_{r=1}^{6}dx_{r}\,\delta(1-\sum x_{r})\,U^{\leavevmode\nobreak\ \varepsilon-2}(V-i\varrho)^{-3\varepsilon}$		(2)

where $V=M^{2}-W/U,\leavevmode\nobreak\ \leavevmode\nobreak\ M^{2}=\sum_{r}m_{r}^{2}x_{r}$;   $U$ and $W$ are polynomials determined by the topology of the corresponding diagram and physical parameters;  $\nu=4-2\varepsilon$ is the space-time dimension; $\varrho$ is a regularization parameter and we can set $\varrho=0$ unless $V$ vanishes in the domain;  ${\mathcal{C}}_{N}$ is the $N$-dimensional unit hypercube. We consider the integral without the prefactor $\Gamma(3\varepsilon).$

We use automatic integration, which is a black-box approach for producing (as outputs) an approximation $Qf$ to an integral $\mathcal{I}=\int_{\mathcal{D}}f(\vec{x})\leavevmode\nobreak\ d\vec{x}$ and an estimate ${\mathcal{E}}f$ of the actual error $Ef=|Qf-\mathcal{I}|,$ with the goal of satisfying an accuracy requirement of the form $|\,Qf-\mathcal{I}\,|\leavevmode\nobreak\ \leq\leavevmode\nobreak\ {\mathcal{E}}f\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \max\,\{\,t_{a}\,,\,t_{r}\,|\mathcal{I}|\,\},$ where the integrand function $f,$ $d$-dimensional region $\mathcal{D}$ and (absolute/relative) error tolerances $t_{a}$ and $t_{r},$ respectively, are specified as part of the input.

For ultra-violet (UV) singularities at the boundaries of the domain, we make use of lattice rules for Quasi-Monte Carlo (QMC) integration [6] with a boundary transformation that is capable of smoothing singularities [7, 8]. Another non-adaptive approach is based on a double-exponential (DE) formula [9, 10, 11, 12]. Adaptive integration may be useful in fairly low dimensions to treat interior singularities by intensive region partitioning around hot spots [13, 14, 15]. We give results obtained with 1D iterated adaptive integration by the program Dqagse from the Quadpack package [13, 16]. In [1] we applied adaptive integration with an extrapolation on the parameter $\varrho$ that adjusts the factor in $V$ in the integrand denominator of Eq. (1). We termed the combined extrapolations in $\varepsilon$ and $\varrho$ as a “double extrapolation.”

Subsequently, Section 2 defines the baseball integral; linear extrapolation and double extrapolation methods are covered in Section 3; results are given in Section 4 and conclusions in Section 5.

A set of 3-loop diagrams with six internal lines, and their associated integrals are calculated by the authors after separating the ultra-violet divergence  [17]. It is our goal to demonstrate our techniques for numerical integration and extrapolation using another 3-loop diagram in the massive case. Hereafter, we assume that all masses of the internal lines equal 1. The derivation is based on a sector decomposition in 5D parameter space; thus the integrals are labeled by permutations of the 5-tuple $(1,2,3,4,5).$ There are $5!=120$ integrals, whereas many of the integrals may coincide through symmetries and the number of integrals to be computed is reduced. We deal with the baseball diagram shown in Figure 1. Then after assuming equal masses, 15 integrals remain. In this paper, we will study the ${\mathcal{I}}_{51234}$ integral of the baseball diagram, which is the integral that has the most singular UV divergence among the 15, with an ${\mathcal{O}}(\frac{1}{\varepsilon^{2}})$ term, and the rest start at the ${\mathcal{O}}(\frac{1}{\varepsilon})$ or ${\mathcal{O}}(1)$ order.

The integral is given by

	$\displaystyle{\mathcal{I}_{51234}}$	$\displaystyle=$	$\displaystyle\int d\Gamma_{X}\leavevmode\nobreak\ t^{-1+2\varepsilon}v^{-1+\varepsilon}u^{\varepsilon}\,G^{-3\varepsilon}f^{-2+\varepsilon}=I_{00}+I_{01}+I_{10}+I_{11}$		(3)
	$\displaystyle I_{00}$	$\displaystyle=$	$\displaystyle\frac{1}{\varepsilon^{2}}\frac{1}{2}\int d\Gamma_{0}\leavevmode\nobreak\ u^{\varepsilon}G_{0}^{-3\varepsilon}f_{0}^{-2+\varepsilon}$		(4)
	$\displaystyle I_{01}$	$\displaystyle=$	$\displaystyle\frac{1}{\varepsilon}\frac{1}{2}\int d\Gamma_{a}\leavevmode\nobreak\ v^{-1+\varepsilon}u^{\varepsilon}\left(G_{a}^{-3\varepsilon}f_{a}^{-2+\varepsilon}-G_{0}^{-3\varepsilon}f_{0}^{-2+\varepsilon}\right)$		(5)
	$\displaystyle I_{10}$	$\displaystyle=$	$\displaystyle\frac{1}{\varepsilon}\int d\Gamma_{b}\leavevmode\nobreak\ t^{-1+2\varepsilon}u^{\varepsilon}\left(G_{b}^{-3\varepsilon}f_{b}^{-2+\varepsilon}-G_{0}^{-3\varepsilon}f_{0}^{-2+\varepsilon}\right)$		(6)
	$\displaystyle I_{11}$	$\displaystyle=$	$\displaystyle\int d\Gamma_{X}\leavevmode\nobreak\ t^{-1+2\varepsilon}v^{-1+\varepsilon}u^{\varepsilon}\times$		(7)
			$\displaystyle\left(G^{-3\varepsilon}f^{-2+\varepsilon}-G_{a}^{-3\varepsilon}f_{a}^{-2+\varepsilon}-G_{b}^{-3\varepsilon}f_{b}^{-2+\varepsilon}+G_{0}^{-3\varepsilon}f_{0}^{-2+\varepsilon}\right)$
		$\displaystyle=$	$\displaystyle\int d\Gamma_{X}\leavevmode\nobreak\ t^{-1+2\varepsilon}v^{-1+\varepsilon}u^{\varepsilon}\left(G^{-3\varepsilon}f^{-2+\varepsilon}-G_{b}^{-3\varepsilon}f_{b}^{-2+\varepsilon}\right)$

where

	$\displaystyle f$	$\displaystyle=$	$\displaystyle 1+w+u+uw+tu+tuw+tuvw+tu^{2}vw$		(8)
	$\displaystyle f_{0}$	$\displaystyle=$	$\displaystyle f|_{t=0,v=0},\leavevmode\nobreak\ \leavevmode\nobreak\ f_{a}=f|_{t=0},\leavevmode\nobreak\ \leavevmode\nobreak\ f_{b}=f|_{v=0}$		(9)
	$\displaystyle q$	$\displaystyle=$	$\displaystyle z(1-z)+wz(1-z)+uz(1-z)+uwz(1-z)+tuz+tuwz+tuvw(1-z)$
			$\displaystyle+tu^{2}vw(1-z)+t^{2}u^{2}vw$

	$\displaystyle q_{0}$	$\displaystyle=$	$\displaystyle q|_{t=0,v=0},\leavevmode\nobreak\ \leavevmode\nobreak\ q_{a}=q|_{t=0},\leavevmode\nobreak\ \leavevmode\nobreak\ q_{b}=q|_{v=0}$		(10)
	$\displaystyle G$	$\displaystyle=$	$\displaystyle(1+t(1+u(1+v(1+w))))-s\frac{q}{f}$		(11)
	$\displaystyle G_{0}$	$\displaystyle=$	$\displaystyle G|_{t=0,v=0},\leavevmode\nobreak\ \leavevmode\nobreak\ G_{a}=G|_{t=0},\leavevmode\nobreak\ \leavevmode\nobreak\ G_{b}=G|_{v=0}$		(12)

and

	$\displaystyle\int d\Gamma_{X}=\int_{0}^{1}dz\int_{0}^{1}dt\int_{0}^{1}du\int_{0}^{1}dv\int_{0}^{1}dw,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \int d\Gamma_{0}=\int_{0}^{1}dz\int_{0}^{1}du\int_{0}^{1}dw$
	$\displaystyle\int d\Gamma_{a}=\int_{0}^{1}dz\int_{0}^{1}du\int_{0}^{1}dv\int_{0}^{1}dw,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \int d\Gamma_{b}=\int_{0}^{1}dz\int_{0}^{1}dt\int_{0}^{1}du\int_{0}^{1}dw$

Note that $f_{0}=f_{a},$ $q_{0}=q_{a}$ and $G_{0}=G_{a};$ thus the integrand of Eq. (5) is zero, and $I_{11}$ simplifies to Eq. (7). Here, $s$ equals $s=P^{2}$, where $P$ is the external momentum; $s$ is measured in units of $m^{2}$.

Results for $s=1$ and $s=4$ are summarized in Table 1, showing leading coefficients of the Laurent series expansion of the integrals $I_{00},I_{!0}$ and $I_{11}.$ The integration uses a lattice rule with $\leavevmode\nobreak\ 10M$ points, with two copies in each coordinate direction, and the Sidi $\Psi_{6}$ transformation [8]. The expansion with respect to $\varepsilon$ is performed with a linear (single) extrapolation. “NI” refers to our “Numerical Integration.”

The “exact” values listed for comparison were obtained analytically and with the Mathematica functions Integrate and NIntegrate. Analytic exact values are given for $C_{-2},$ $C_{-1}$ and $C_{0},$ and numerical values for $C_{1}.$ The $C_{-2}$ term appears only in $I_{00}$ and is $s$-independent. Similarly, $C_{-1}$ of $I_{10}$ is $s$-independent. When the Mathematica Integrate function fails at the evaluation of $C_{0}$ and $C_{1}$ for $I_{10}$ and $I_{11}$ for $s=1$ and 4, NIntegrate is used. The lattice rule integrations were carried out on an Nvidia Quadro GV100 GPU, each taking time of a fraction of a second for about 20 integrations per problem, although less were used.

Sample results using double extrapolations are further given in Table 2 for $I_{00}$ with $s=4.5,\,5,\,7$ and 10. The $\varrho$ values for the extrapolation follow a geometric sequence $\varrho_{j}=2^{-j}$ for $j=1,2,\ldots,11,$ whereas the $\varepsilon$ values form a Bulirsch type sequence [26], $\varepsilon_{\ell}=1/4,1/6,1/8,1/12,1/16,1/24,\ldots.$ The $\varepsilon$ results listed in Table 2 are incurred when the estimated error (based on successive differences) no longer decreases. $\varepsilon$-ext denotes the number of $\varepsilon$ extrapolations to obtain the result for this value of $s.$ T[s] denotes the corresponding time (for integration) incurred through this number of $\varepsilon$-ext (not including the iteration where the accuracy stops improving). The integration time dominates the overall time. The larger values of $s$ require more extrapolations for convergence, hence utilize more time. These computations were performed (sequentially) on an x86_64 machine under GNU/ Linux. The method is also effective for $I_{10}$ (4D), but needs parallelization for $I_{11}$ (5D).

As part of our discussion, let’s compare the performance of single and double extrapolations for the same problem. To conclude initially, we observed that double extrapolation may be justified to achieve better accuracy. This is illustrated for $I_{00}$ with $s=7.$ The integrations are done by iterated adaptive integration using the routine Dqagse from Quadpack [13, 16] with a maximum of 50 subdivisions in each coordinate direction, relative error tolerance $10^{-9}$ in the outer coordinate and $5\times 10^{-10}$ in subsequent directions.

A double extrapolation with $\varepsilon=1/1.15^{\ell},\leavevmode\nobreak\ \varrho=1/1.2^{k},25\leq\ell\leq 30,\leavevmode\nobreak\ 10\leq k\leq 20$ yields for $C_{-2},C_{-1},C_{0}$ and $C_{1}$ at the $5^{th}$ linear system: 0.1250000066$\pm$7.5e-08, 0.2687828385$\pm$1.7e-05,
-2.7440527816$\pm$1.5e-03, -4.5561503884. The time for all six $\varepsilon$-cycles is 18.1 seconds. For the same problem, a single extrapolation with  $\varepsilon=1/1.15^{\ell},\ell\geq 25,\leavevmode\nobreak\ \leavevmode\nobreak\ \varrho=0$   gives as the best result observed (at system 4): 0.1249996422$\pm$4.7e-06, 0.2688439263$\pm$7.8e-04, -2.7476802260$\pm$4.8e-02, -4.4743466989. In order to compare the time with the double extrapolation, the result for $C_{-2},C_{-1},C_{0},C_{1}$ at system 5 is: 0.1250017027$\pm$2.1e-06, 0.2683866114$\pm$4.6e-04, -2.7074720045$\pm$1.5e-02, -6.2248242790, with a total time (for integrations 1 to 6) of 4.9 seconds.

Finally, graphs are displayed in Figure 2 for the real and imaginary part of the total integral ${\mathcal{I}}_{51234}$ of Eq. (2), where the data were obtained using DE integration and single extrapolation with regard to $\rho$ after an expansion in $\varepsilon$  [17].   

Whereas symbolic or symbolic/numerical calculations are performed for some challenging problems using existing software packages, we focus on the development of fully numerical methods for the evaluation of Feynman loop integrals. The integration strategies adhere to black-box approaches for generating an approximation, assuming little or no knowledge of the problem, apart from the specification of the integrand function. We demonstrated efficient strategies based on lattice rules evaluated on GPUs and double-exponential integration, as well as iterated adaptive integration.

\ack

We acknowledge the support by JSPS KAKENHI Grant Number JP20K11858, JP20K03941 and JP21K03541, and the National Science Foundation Award Number 1126438 that funded work on multivariate integration.