The Method of Arbitrarily Large Moments to Calculate Single Scale Processes in Quantum Field Theory

TL;DR

This paper introduces a novel method leveraging differential and difference equations to efficiently compute Mellin moments of single scale quantities in quantum field theory, enabling precise numerical and potential analytic reconstructions.

Contribution

The paper presents a new approach to calculate large sets of Mellin moments using differential/difference equations derived from integration-by-parts identities, simplifying complex quantum calculations.

Findings

Allows calculation of many Mellin moments efficiently

Enables analytic reconstruction of physical quantities

Provides highly precise numerical representations

Abstract

We device a new method to calculate a large number of Mellin moments of single scale quantities using the systems of differential and/or difference equations obtained by integration-by-parts identities between the corresponding Feynman integrals of loop corrections to physical quantities. These scalar quantities have a much simpler mathematical structure than the complete quantity. A sufficiently large set of moments may even allow the analytic reconstruction of the whole quantity considered, holding in case of first order factorizing systems. In any case, one may derive highly precise numerical representations in general using this method, which is otherwise completely analytic.

Tables

Table 1. Table 1: The number of moments needed to find the associated difference equations.

	${\epsilon }^{-3}$	${\epsilon }^{-2}$	${\epsilon }^{-1}$	${\zeta }_2$	${\epsilon }^0$	${\zeta }_2$	${\zeta }_3$
$J_1$	8	24	48	8	89	24	8
$J_2$	8	23	42	8	81	21	8
$J_3$	-	8	19	-	63	-	5

Table 2. Table 2: The order and degree of the associated difference equations.

	${\epsilon }^{-3}$		${\epsilon }^{-2}$		${\epsilon }^{-1}$		${\zeta }_2$		${\epsilon }^0$		${\zeta }_2$		${\zeta }_3$
	$O$	$R$	$O$	$R$	$O$	$R$	$O$	$R$	$O$	$R$	$O$	$R$	$O$	$R$
$J_1$	1	2	2	5	3	11	1	2	4	20	2	5	1	2
$J_2$	1	2	2	5	3	10	1	2	4	18	2	5	1	2
$J_3$	-	-	1	2	2	4	-	-	4	12	-	-	1	1

Equations

$${\hat{I}}_i(x)=\sum _{n=0}^{\infty }{I}_i(n)x^n=\sum _{n=0}^{\infty }\left(\sum _{k=o}^{\infty }{I}_i^{(k)}(n){\epsilon }^k\right)x^n$$

$$I_i^{(k)}(0),I_i^{(k)}(1),I_i^{(k)}(2),\dots ,I_i^{(k)}(s)\in \mathbb{R}$$

$$D_x\left(\begin{array}{c}{\hat{I}}_1(x)\\ {\hat{I}}_2(x)\\ \dots \\ {\hat{I}}_m(x)\end{array}\right)=A\left(\begin{array}{c}{\hat{I}}_1(x)\\ {\hat{I}}_2(x)\\ \dots \\ {\hat{I}}_m(x)\end{array}\right)+\left(\begin{array}{c}{\hat{r}}_1(x)\\ {\hat{r}}_2(x)\\ \dots \\ {\hat{r}}_m(x)\end{array}\right)$$

$${\hat{r}}_i(x)=\sum _{n=0}^{\infty }{r}_i(n)x^n=\sum _{n=0}^{\infty }\left(\sum _{k=o}^{\infty }{r}_i^{(k)}(n){\epsilon }^k\right)x^n$$

$$\sum _{k=0}^m{a}_k(x,\epsilon )D_x^{(k)}{\hat{I}}_1(x)=\sum _{i=1}^m{d}_i(x,\epsilon ){\hat{r}}_i(x),$$

$${\hat{I}}_i(x)$$

$$\sum _{k=0}^d{b}_k(n,\epsilon )I_1(n+k)=\rho (n,\epsilon ),$$

$$\rho (n,\epsilon )=\sum _{k=o}^{\infty }{\rho }^{(k)}(n){\epsilon }^k,$$

$$I_1(n)=\sum _{k=o}^{\infty }{I}_1^{(k)}(n){\epsilon }^k$$

$$\sum _{k=0}^{d^{\prime }}{b}_k(n,0)I_1^{(o)}(n+k)={\rho }^{(o)}(n)$$

$$I_1^{(o)}(n)\leftarrow \frac{{\rho }^{(o)}(n^{\prime })-\sum _{k=0}^{d^{\prime }-1}{b}_k(n^{\prime },0)I_1^{(o)}(n^{\prime }+k)}{b_{d^{\prime }}(n^{\prime },0)}$$

$$\sum _{k=0}^{d^{\prime }}{b}_k(n,\epsilon )I_1^{\prime }(n+k)=\sum _{k=o+1}^{\infty }{\rho }^{\prime (k)}(n){\epsilon }^k$$

$$\sum _{k=0}^O{a}_k(n)F(n+k)=0$$

$$\frac{1}{4^n}\left(\genfrac{}{}{0ex}{}{2n}{n}\right)\left\{\sum _{k=1}^n\frac{4^k{S}_1(k-1)}{\left(\genfrac{}{}{0ex}{}{2k}{k}\right)k^2}-7{\zeta }_3\right\},$$

Full text

††preprint: DESY 17-009, DO-TH 17/01

I INTRODUCTION

Single scale higher order QED and QCD calculations in the massless Vermaseren:2005qc ; Ruijl:2016pkm and massive Bierenbaum:2009mv cases at fixed Mellin moment $n$ are given by polynomials of rational numbers and a series of a few special constants, the multiple $\zeta$-values, and possible generalizations thereof MZVX . This is irrespectively the case, whether or not the functional representation for general values of $n$ obeys an equation factorizing in first order, Vermaseren:2005qc , or not ELL . If one has access foot1 to an algorithm, through which a large number of moments, e.g. $N=2000$ or larger, can be calculated, the method of guessing, cf. GUESSING , for holonomic problems, which often appear in physics applications, allows one to gain one large difference equation describing the corresponding problem Blumlein:2009tj . If this equation is solvable in difference ring theory Summation one finds the solution for general values of $n$ for this quantity without any further assumptions, e.g. made in Ruijl:2016pkm ; Velizhanin:2014fua . In any case, significantly more moments will allow one to constrain the considered quantity much better numerically using approximation methods, e.g. through Chebyshev-polynomials or other interpolation methods.

In the following we describe an algorithm through which the system of differential equations, or associated to it, that of difference equations, available by the integration-by-parts relations Chetyrkin:1981qh , can be used to compute a large number of Mellin moments for the master integrals, and through them for the whole problem. The solution of the associated difference equations will need a relatively low number of initial values which have to be provided. The corresponding sequences of rational numbers mentioned above actually form the problematic part in gaining the general $n$ result from the moments, since very involved, and in some cases yet unknown, functions span the corresponding sequences. In any case, the method allows either to find the one-dimensional distribution from a large but finite amount of moments, or at least to constrain it numerically at high accuracy.

We will give a brief illustration of the algorithm in case of a system of massive 3-loop master integrals, for which first a larger number of moments is generated, corresponding difference equations are found and solved for general values of the Mellin variable $n$.

II The Algorithm

Single scale master integrals can be represented as analytic functions $\hat{I}_{i}(x)=\sum_{n=0}^{\infty}I_{i}(n)x^{n}$ with $1\leq$ $i\leq m$. We aim at computing a large number of coefficients $I_{i}(0),I_{i}(1),\dots,I_{i}(s)$. Usually the coefficients $I_{i}(n)$ depend on the dimensional parameter $\varepsilon$ which itself can be expanded in a Laurent series in $\varepsilon$ of a certain order $o\in\mathbb{Z}$. We are interested in calculating the coefficients $I_{i}^{(k)}(n)\in\mathbb{R}$ of the expansions

[TABLE]

up to a certain degree in $\varepsilon$, $t_{i}\in\mathbb{Z}$. More precisely, we want to compute for $1\leq i\leq m$ the initial values

[TABLE]

for the $\varepsilon$-orders $o\leq k\leq t_{i}$.

In our approach we rely on the property that these unknown functions $\hat{I}_{i}(x)$ are usually described by a coupled system of first oder linear differential equations

[TABLE]

with $D_{x}=d/dx$ and where $A$ is an $m\times m$ matrix with entries consisting of polynomials in $\varepsilon$ and $x$ and the entries $\hat{r}_{i}$ can be given in form of the expansions

[TABLE]

with $r_{i}^{(k)}(n)\in\mathbb{R}$ for $n\in\mathbb{N}$. Another important assumption is that the coefficients $r_{i}^{(k)}(n)$ can be determined efficiently by using either the method under consideration in a recursive fashion or by using, e.g., symbolic summation and integration techniques Summation ; VLadders . Further, we assume that for reasonable small numbers $s^{\prime}_{i}$ the coefficients $I_{i}^{(k)}(n)$ with $0\leq n\leq s^{\prime}_{i}$ can be computed up to certain $\varepsilon$-degrees $t^{\prime}_{i}$, i.e., $o\leq k\leq t^{\prime}_{i}$ as a preprocessing step. Given this input we introduce the following efficient algorithm that computes the coefficients (2) up to large values of $n$.

1. (1)

Using decoupling algorithms UNCOUPL ; OreSys transform the system (3) symbolically to one scalar linear differential equation

[TABLE]

where the $a_{i}(x,\varepsilon)$ and $r_{i}(x)$ are polynomials in $x$ and $\varepsilon$; in addition, one obtains identities of the form

[TABLE]

for $2\leq i\leq n$ where the $e_{i,j,k}(x,\varepsilon)$ and $f_{i,j,k}(x,\varepsilon)$ are rational functions in $x$ and $\varepsilon$.

The algorithm proceeds now as follows: Compute in steps (2)-(4) the initial values for $\hat{I}_{1}(x)$ by using the scalar differential equation (5), and compute afterwards in step (5) the initial values for the remaining integrals with the given formulas (6).

1. (2)

Plug $\hat{I}_{1}(x)=\sum_{n=0}^{\infty}I_{1}(n)x^{n}$ and $\hat{r}_{i}(x)=\sum_{n=0}^{\infty}r_{i}(n)x^{n}$ into (5) and eliminate $D_{x}$ by using the property $D_{x}\sum_{n=0}^{\infty}h(n)x^{n}=\sum_{n=1}^{\infty}n\,h(n)x^{n-1}$ for a power series $\sum_{n=0}^{\infty}h(n)x^{n}$. Then by coefficient comparison w.r.t. $x^{n}$ and using appropriate shifts one gets a linear recurrence of the form

[TABLE]

with $\rho(n,\varepsilon)=\sum_{j=1}^{m}\sum_{k=0}^{l}g_{j,k}(n,\varepsilon)r_{j}(n+k)$ for some $l\in\mathbb{N}$ where the $b_{k}(n,\varepsilon)$ and $g_{j,k}(n,\varepsilon)$ are polynomials in $n$ and $\varepsilon$. Finally, divide the equation by a factor $\varepsilon^{u}$ for some $u\geq 0$ in order to obtain updated polynomials $b_{k}(n,\varepsilon)$ in $n,\varepsilon$ where not all $b_{k}(n,0)$ with $0\leq k\leq d$ are zero; the $g_{j,k}(n,\varepsilon)$ are now polynomials in $n$ and Laurent polynomials in $\varepsilon$. 2. (3)

We write the right hand side of (7) in the expanded representation

[TABLE]

with $\rho^{(k)}(n)\in\mathbb{R}$. Since the $r_{i}^{(k)}(n)\in\mathbb{R}$ in (4) can be computed efficiently by assumption, the coefficients $\rho^{(k)}(n)$ for sufficiently large $k\in\mathbb{Z}$ and $n\in\{0,1,\dots,s\}$ can be obtained explicitly without any cost. 3. (4)

We proceed as follows; compare BKSS:12 , Lemma 1. Plugging

[TABLE]

into (7) and doing coefficient comparison w.r.t. $\varepsilon^{o}$ yield the constraint

[TABLE]

for some $d^{\prime}\leq d$ with $b_{d^{\prime}}(n)\neq 0$. Choose $\delta\in\mathbb{N}$ such that $b_{d^{\prime}}(n)\neq 0$ for all $n\in\mathbb{N}$ with $n\geq\delta$. Then we can compute with the first values $I_{1}^{(o)}(0),\dots,I_{1}^{(o)}(n+d^{\prime}+\delta-1)$ and the rule

[TABLE]

for $n^{\prime}=n-d^{\prime},n\geq d^{\prime}+\delta$ all the other values in linear time. Now insert (9) with the explicitly computed values $I_{1}^{(o)}(n)$ with $0\leq n\leq s$ into (7) and move these given values to the right hand side. This yields

[TABLE]

for $I^{\prime}_{1}(n)=\sum_{k=o+1}^{\infty}I_{1}^{(k)}(n)$ where the $\rho^{\prime(k)}(n)$ for $0\leq n\leq s$ and sufficiently large $k$ are given explicitly. Now we are in the position to repeat this tactic iteratively to compute the remaining coefficients: next $I_{1}^{o+1}(n)$ with $0\leq n\leq s$, afterwards $I_{1}^{o+2}(n)$ with $0\leq n\leq s$, and eventually $I_{1}^{t_{i}}(n)$ with $0\leq n\leq s$. 4. (5)

Now expand the rational functions $f_{i,j,k}(x,\varepsilon)$ and $g_{i,j,k}(x,\varepsilon)$ in (6) in a Laurent-series expansion w.r.t. $\varepsilon$ and compute their coefficients as power series in $x$ up to the necessary orders. Finally, combine all the explicitly given (finite) expansions in (6) using component-wise addition and Cauchy-products which yield the coefficients $I_{i}^{(k)}(n)$ for $2\leq i\leq m$, $o\leq k\leq t_{i}$ and $0\leq n\leq s$.

The following remarks are in order; related considerations have been applied also to our algorithms to solve coupled systems in terms of nested sums over hypergeometric products SymbolicUncouplingMethod1 ; VLadders ; SymbolicUncouplingMethod2 ; SymbolicUncouplingMethod3 .

1. (i)

Our algorithm requires that sufficiently many initial values/coefficients of $\hat{r}_{i}(x)$ and $\hat{I}_{1}(x)$ up to the right $\varepsilon$-order are computed as a preprocessing step. The necessary bounds for these numbers can be determined by analyzing the formulas (7) and (6) accordingly. 2. (ii)

For simplicity, we assumed that only one scalar differential equation arises, as it happens in most examples. In general, several scalar differential equations might arise, i.e., steps (2)–(4) have to be executed several times. 3. (iii)

One can choose any $\lambda$ with $1\leq\lambda\leq m$ to determine a scalar differential equation in $\hat{I}_{\lambda}(x)$. Different choices of $\lambda$ might lead to different recurrences (7) with different orders $d$ and formulae (6) of different size. In our calculations we analyze for each $\lambda$ ($1\leq\lambda\leq m$) the obtained symbolic formulae and choose that $\lambda$ for the calculation steps (3)-(5) that serves us best, e.g., to minimize the required $\varepsilon$-orders for the $\hat{r}_{i}(x)$ or minimize the the recurrence order $d$. 4. (iv)

Often it is a challenge to determine the first initial values of $I_{i}^{k}(n)$ to activate the recurrence formula (10). Thus it is highly desirable to keep the order $d$ in (7) as small as possible. Since $d$ gets smaller if the degrees of the $a_{i}(x,\varepsilon)$ w.r.t. $x$ in (5) can be made smaller, the following refinement can be applied. Divide (5) by the polynomial $h=\gcd_{x}(a_{0},\dots,a_{n})$ in $x$. In most applications this reduces the degrees of the polynomials $a_{i}(x,\varepsilon)$ w.r.t. $x$ heavily and thus produces recurrences with much lower order, e.g., $d=4$ instead of $d=20$ in (7). The price to be paid is to determine the coefficients $\rho^{(k)}(n)$ in (8) by extracting the $n$th coefficients from the right hand side of (5). Here the more involved calculation steps similarly as given in step (5) must be carried out.

All these aspects have been implemented within our package SolveCoupledSystem VLadders ; SymbolicUncouplingMethod2 ; SymbolicUncouplingMethod3 which uses subroutines of the summation package Sigma Summation and the uncoupling package OreSys OreSys .

III An Example

Let us consider a set of three 3-loop master integrals $\{\hat{J}_{1}(x,\varepsilon),\hat{J}_{2}(x,\varepsilon),\hat{J}_{3}(x,\varepsilon)\}$ which obey a $3\times 3$ inhomogeneous linear coupled system, cf. (3), as an illustrative example. They contribute to the massive operator matrix element $A_{gg}^{(3)}$ Ablinger:2014uka , and are given in Eqs. (3.62–3.64). In a first step, we will use our above method to compute a large series of fixed moments for $J_{1,2,3}(n,\varepsilon)$, and in a second step we will use this information to guess recurrences for these functions. In our example the integrals $J_{i}$ have poles up to $\varepsilon^{-3}$ and we would like to compute the moments up to the constant term $O(\varepsilon^{0})$. Since the number of moments $n_{0}$ needed to guess a recurrence is not known a priori, the moments are computed up to a reasonable number. If this is not sufficient, further moments can be computed reusing the previous calculations. E.g., using our algorithm the time to generate $n_{0}=500$ moments amounts to 48 sec and $n_{0}=2000$ moments require 569 sec in the present example. In a nutshell, higher moments can therefore easily be produced.

The guessing method to find the minimal difference equation requires the numbers $n_{0}$, depending on $\varepsilon$ and whether the respective set corresponds to the purely rational term or a respective contribution $\propto\zeta_{i}$. The values are summarized in Table I. Here the highest number turns out to be $n_{0}=89$. The guessing method GUESSING yields a difference equation

[TABLE]

for each power in $\varepsilon$. Here $O$ and $R$ denote the order and degree of the difference equation, i.e. the number of shift-operators and the maximal degree of the polynomials $a_{k}(n)$. The orders and degrees of the different difference equations are listed in Table II. The number of moments always includes a sufficiently large safety margin to validate the corresponding difference equations, which may even be enlarged. We generate additional moments according to their order to provide the needed initial values to solve the difference equations.

In the present case, using the algorithm of Ref. VLadders , all difference equations can be solved completely using difference ring theory Summation . Besides rational terms in $n$ the result for $J_{1}(n)$ and $J_{2}(n)$ is spanned by the harmonic sums HSUM up to $S_{1,2}(n)$. $J_{3}(n)$, furthermore, contains the finite binomial sum

[TABLE]

multiplied by a rational term in $n$, see also Ablinger:2014bra . All appearing letters building the nested sums forming the master integrals are found automatically. The above example is a simpler one. Let us mention that, for comparison, the computational time for the complete pole terms of more involved massive 3-loop examples requires the knowledge of about $\sim 1000$ even moments Blumlein:2009tj and is estimated to amount to several weeks PREP .

In case the difference equations would turn out not to be first order factorizable using the algorithm of Ref. VLadders , at least the first order factorizable terms would be gained analytically, leaving the non-factorizable part behind for further analysis using different methods, like those for elliptic integrals and their potential generalization.

Having the above number of moments, Eq. (1) yields a first numerical approximation of $\hat{J}_{k}(x)$ also up to the required power in $\varepsilon$.

The present method is suitable to obtain precise numerical representations in case of various current massless and massive calculations at 3- and 4-loop order. In case one may solve the associated difference equations algebraically, one will even obtain the complete analytic result without any prejudice, i.e. a sufficiently large set of scalar moments allows one to fully reconstruct a one-dimensional distribution. This automatic method is therefore suited to carry out many more higher loop calculations in contemporary elementary particle physics in a completely or at least widely analytic form.

We would like to thank J. Ablinger, A. Behring, A. De Freitas, M. Kauers, and P. Marquard for discussions. This work was supported in part by the Austrian Science Fund (FWF) grant SFB F50 (F5009-N15), the European Commission through contract PITN-GA-2012-316704 (HIGGSTOOLS).