Large scale analytic calculations in quantum field theories

TL;DR

This paper surveys the mathematical structures and calculation methods used in high-order perturbative quantum field theory, focusing on zero- and single-scale quantities and their associated function spaces.

Contribution

It provides a comprehensive overview of the mathematical frameworks and techniques for perturbative calculations in relativistic quantum field theories.

Findings

Detailed classification of function spaces used in calculations

Analysis of mathematical structures underlying perturbative methods

Overview of calculation techniques for higher-order corrections

Abstract

We present a survey on the mathematical structure of zero- and single scale quantities and the associated calculation methods and function spaces in higher order perturbative calculations in relativistic renormalizable quantum field theories.

Equations

$$\left\{\ln (2),{\zeta }_2,{\zeta }_3,{\mathrm{Li}}_4\left(\frac{1}{2}\right)\right\}.$$

$$0.3395454690873598695906678484608602061387815339795751791304750$$

$$222490137419723806082682624196443182167020255697096551752247012$$

$$11749559277$$

$$H_{-1,0,0,1}(1)=-\frac{1}{12}{\ln }^4(2)+\frac{1}{2}{\ln }^2(2){\zeta }_2+\frac{3}{5}{\zeta }_2^2-\frac{3}{4}\ln (2){\zeta }_3-2{\mathrm{Li}}_4\left(\frac{1}{2}\right).$$

$$B(a_1,a_2)$$

$${}_3{F}_2(a_1,a_2,a_3;b_1,b_2;x)$$

$$\frac{1}{(a+b{)}^{\alpha }}=\frac{1}{\Gamma (\alpha )}\frac{1}{2\pi i}{\int }_{-i\infty }^{i\infty }dz\Gamma (\alpha +z)\Gamma (-z)\frac{b^z}{a^{\alpha +z}},\alpha \in \mathbb{R},\alpha >0,$$

$$\displaystyle\frac{d}{dx}\left(\begin{array}[]{c}f_{1}\\ \vdots\\ f_{n}\end{array}\right)=\left(\begin{array}[]{ccc}A_{11}&\ldots&A_{1,n}\\ \vdots&&\vdots\\ A_{n1}&\ldots&A_{n,n}\end{array}\right)\left(\begin{array}[]{c}f_{1}\\ \vdots\\ f_{n}\end{array}\right)+\left(\begin{array}[]{c}g_{1}\\ \vdots\\ g_{n}\end{array}\right),$$

$$\sum _{k=0}^n{p}_{n-k}(x)\frac{d^{n-k}}{d{x}^{n-k}}{f}_1(x)=\overline{g}(x),$$

$$f_1(x)$$

$$g_k(x)$$

$${\partial }_{m}f(\epsilon ,x_n)=A_m(\epsilon ,x_n)f(\epsilon ,x_n)$$

$$A_m^{\prime }=B^{-1}{A}_{m}B-B^{-1}({\partial }_{m}B),$$

$${\partial }_{m}f(\epsilon ,x_n)=\epsilon A_m(x_n)f(\epsilon ,x_n),$$

$$G(N)=\mathbf{M}[f(x)](N)={\int }_0^{1}dx{x}^{N-1}f(x).$$

$$\sum _{l=0}^m{P}_l(x,\epsilon )\frac{d^l}{d{x}^l}I(x,\epsilon )$$

$$\sum _{l=0}^m{R}_l(N,\epsilon )I(N+l,\epsilon )$$

$$\psi (x)={\psi }_1^{(0)}(x)\left[C_1-\int dx{\psi }_2^{(0)}(x)\frac{N(x)}{W(x)}\right]+\{1\to 2\},$$

$${\mathbb{H}}_{a_1,...,a_{m-1};a_m,F_m(r(y_m)),a_{m+1},...a_q}(x)=$$

$${\int }_0^{x}d{y}_1{f}_{a_1}(y_1){\int }_0^{y_1}d{y}_2...{\int }_0^{y_{m-1}}d{y}_m{f}_{a_m}(y_m)F_m[r(y_m)]H_{a_{m+1},...,a_q}(y_m),$$

$$k^2=r(x)=\frac{{\vartheta }_2^4(q)}{{\vartheta }_3^4(q)},$$

$$\eta (\tau )=q^{\frac{1}{12}}\prod _{k=1}^{\infty }(1-q^{2k}).$$

$${\mathrm{ELi}}_{n,m}(x,y)=\sum _{k=1}^{\infty }\sum _{l=1}^{\infty }\frac{x^k}{k^n}\frac{y^l}{l^m}{q}^{kl}$$

$$\frac{1}{1-x(1-\eta )},\frac{\sqrt{x(1-x)}}{\eta +x(1-\eta )},\sqrt{x(1-\eta (1-x)},\eta \in [0,1].$$

$$\frac{x}{\sqrt{1-x^2}\sqrt{1-k^2{x}^2}},\frac{x}{\sqrt{1-x^2}\sqrt{1-k^2{x}^2}(k^2(1-x^2(1-z^2))-z^2)},$$

$$S_{b,\vec{a}}(N)=\sum _{k=1}^N{g}_b(k)S_{\vec{a}}(n),S_{\varnothing }=1,g_c\in \bar{\mathfrak{A}},$$

$$\sum _{i=1}^N\left(\genfrac{}{}{0ex}{}{2i}{i}\right)(-2{)}^i\sum _{j=1}^i\frac{1}{j\left(\genfrac{}{}{0ex}{}{2j}{j}\right)}{S}_{1,2}\left(\frac{1}{2},1\right)(j)={\int }_0^{1}dx\frac{(-x{)}^N-1}{x+1}\sqrt{\frac{x}{8-x}}$$

$$\displaystyle\times\Bigl{[}{\rm H}_{\sf w_{12},1,0}(x)-2{\rm H}_{\sf w_{13},1,0}(x)-\zeta_{2}\left({\rm H}_{\sf w_{12}}(x)-2{\rm H}_{\sf w_{13}}(x)\right)\Bigr{]}$$

$$-\frac{5{\zeta }_3}{8\sqrt{3}}{\int }_0^{1}dx\frac{(-2x{)}^N-1}{x+\frac{1}{2}}\sqrt{\frac{x}{4-x}}+c_1{\int }_0^{1}dx\frac{(-8x{)}^N-1}{x+\frac{1}{8}}\sqrt{\frac{x}{1-x}},$$

$${\mathrm{H}}_{b,\vec{a}}(x)={\int }_0^{x}dy{f}_b(y){\mathrm{H}}_{\vec{a}}(y),{\mathrm{H}}_{\varnothing }=1,f_c\in \mathfrak{A},$$

$$P(x)=\prod _{k=1}^n(x-a_k)\prod _{l=1}^m(x^2+b_{l}x+c_l),a_k,b_l,c_l\in \mathbb{R}$$

Full text

1 Introduction

Precise theoretical predictions within the Standard Model of elementary particles are indispensable for the concise understanding of the fundamental parameters of this physical theory and the discovery of its potential extensions. At the experimental side highly precise measurements exist at $e^{+}e^{-}$, $ep$ and $pp$–colliders as at LEP, HERA, and the LHC. In the near future the high luminosity phase of the LHC will even provide much more precise data. Other facilities, like the ILC ILC and a possible FCC FCC , are currently planned. During the last three decades enormous efforts have been made to calculate key observables measured at these colliders at higher and higher accuracy, to meet the challenge provided by the accuracy of the experiments.

For zero–scale quantities currently analytic massless calculations can be performed at the five–loop and for massive calculations at the four-loop level. Single scale calculations are performed in both cases at the three–loop level. To perform these large scale calculations very demanding efforts are needed at the side of their automation, computer–algebraic implementation, and the use of highly efficient mathematical technologies. Therefore, the present problems can only be solved within a very close interdisciplinary cooperation between experts in all these different fields and it cannot be the sole tasks for theoretical physicists anymore.

While at one–loop order the mathematical solution for many scattering processes has been known early, cf. tHooft:1978jhc ; tHooft:1973wag ; Veltman:1994wz , systematic representations at higher loop order turned out to be more difficult. The core problem concerns the analytic integration of Feynman parameter integrals. Here integration is understood as anti–differentiation. An essential question is to determine the final solution space to which the respective integrals do belong and its mathematical structure, and to find the irreducible objects through which the corresponding integrals are represented. Furthermore, one needs efficient mathematical and computer–algebraic technologies to map the given Feynman parameter integrals into the latter quantities.

In this paper we give a survey on the main technological steps to calculate higher loop zero– and single–scale quantities in renormalizable quantum field theories, with the focus on analytic integration techniques and the occurring function spaces. The systematic theory of integration in this field is vastly developing and many more new structures are expected to be revealed in the future at higher loop levels and by considering the production of more particles in the final state of the respective scattering processes. These calculations are needed to obtain stable theoretical predictions for the experimental precision measurements at the present and future colliders, which operate at high luminosity.

The paper is organized as follows. In Section 2 we summarize the main steps in multi–loop perturbative calculations. Different methods used in symbolic calculations of zero– and single–scale Feynman parameter integrals are described in Section 3. In Section 4 a hierarchy of function spaces, mainly for single–scale integrals, is discussed which emerge in present multi–loop calculations. Here we consider as well the representations in Mellin–$N$ and $x$–space. Section 5 contains the conclusions.111 For other recent surveys on integration methods for Feynman integrals see Weinzierl:2010ps ; Ablinger:2013eba ; Ablinger:2013jta ; Weinzierl:2013yn ; Duhr:2014woa ; Blumlein:2018cms .

2 Main Steps in Multi–Loop Perturbative Calculations

In most of the large projects, which are currently dealt with, the Feynman diagrams are generated using packages like QGRAF Nogueira:1991ex and performing the color algebra for the gauge groups using Color vanRitbergen:1998pn . Standardized algorithms to obtain Feynman parameterizations exist, cf. e.g. NAKANISHI ; LEFSCHETZ ; Bogner:2010kv . At growing complexity, to perform the Dirac- and spin-algebra will be a challenge even to FORM Vermaseren:2000nd ; Tentyukov:2004hz ; Tentyukov:2007mu ; Ruijl:2017dtg . One further maps the set of the contributing Feynman integrals to master integrals using the integration–by–parts (IBP) technique Chetyrkin:1981qh based on Laporta’s algorithm Laporta:2001dd , of which several implementations exist, cf. e.g. Smirnov:2008iw ; Smirnov:2019qkx ; Studerus:2009ye ; vonManteuffel:2012np ; MARSEID and others. The remaining main step is then the integration of the master integrals. One possibility to inspect the problem on hand, is to analyze the associated system of first order differential equations for the master integrals. Sometimes it is also useful to consider, in addition, the related system of linear difference equations. One may decouple these systems using the algorithms implemented in the packages ORESYS ; BCP13 , as e.g. Zürcher’s algorithm Zuercher:94 . This leads to a single differential equation or difference equation of large order and degree and associated determining equations for the remaining master integrals. If the former equations can be factored at first order, it is known that the master integrals can be obtained in terms of indefinitely nested sums or iterated integrals over certain alphabets, which are revealed in the solution process, e.g. using difference field and ring theory Karr:1981 ; Bron:00 ; Schneider:01 ; Schneider:04a ; Schneider:05a ; Schneider:05b ; Schneider:07d ; Schneider:10b ; Schneider:10c ; Schneider:15a ; Schneider:08c ; Schneider:08d ; Schneider:08e , algorithmically implemented in the package Sigma Schneider:2007a ; Schneider:2013a . This applies to a wide class of physical cases. Most of the integration and summation methods described in Section 3 apply to them and allow to obtain the integrals analytically in terms of the mathematical functions described in Section 4. Finally, efficient numerical representations of these functions have to be provided to obtain numerical predictions of the different observables for the experiments.

3 Symbolic Integration of Feynman Parameter Integrals

In the following we summarize main aspects of the analytic integration of multi–loop Feynman parameter integrals. Of course these integrals can also be evaluated numerically, without observing their particular analytic structure, to some accuracy and methods exist to separate the different pole contributions in $\varepsilon$, cf. e.g. Binoth:2000ps ; Nagy:2006xy ; Anastasiou:2007qb ; Smirnov:2008py ; Carter:2010hi ; Smirnov:2009pb ; Becker:2010ng ; Becker:2012aq ; Becker:2011vg ; Smirnov:2015mct ; Borowka:2018dsa , which we will not discuss in the following. These methods play a role, however, also for testing analytic results. In calculating all the integrals required to solve a large scale problem, it is usually necessary to combine different analytic methods, at least for the sake of efficiency. This requirement finally led to the creation of these methods. In the future even more and further refined technologies will be needed to solve more enhanced problems. Finally, one ends up with sets of irreducible functions which span the solutions, see Section 4. The numerical representation of these functions is necessary and will be discussed in Section 4.3.

Non of the different techniques described in the following are of universal character. In particular the solution of the most advanced problems will need a combined and sensible use of various of them. All of them have to be handled with care to achieve a steady stepwise reduction of the problem on hand and to avoid to enlarge the complexity, given the limited time and memory resources for the corresponding computer algebraic calculations. This will also apply to future developments, since more complex calculations will require further new and advanced technologies.

Many of the formalisms described below lead to summation problems. Their solution requires dedicated and efficient algorithms in difference field theory as implemented in the packages Sigma Schneider:2007a ; Schneider:2013a , EvaluateMultiSums and SumProduction Ablinger:2010pb ; Blumlein:2012hg ; Schneider:2013zna , see also Schneider:19 .

3.1 The PSLQ Method

The PSLQ method applies to the solution of zero dimensional quantities, i.e. physical quantities given by pure numbers. If the pool of constants is known or can be guessed over which the corresponding quantity has a polynomial representation over $\mathbb{Q}$, a highly precise numerical representation of the quantity and the individual monomials allows to determine the corresponding rational coefficients, cf. PSLQ . This method has been applied recently in a massive calculation of the five–loop QCD $\beta$-function Luthe:2017ttg . Here the individual master integrals certainly contain also constants of elliptic nature and probably beyond. However, they all cancel in the final result, which is spanned by multiple zeta values (MZVs) Borwein:1999js ; Blumlein:2009cf , more precisely by $\{\zeta_{2},\zeta_{3},\zeta_{5},\zeta_{7}\}$, beyond pure rational terms. Let us illustrate the method by an example. We would like to determine the harmonic polylogarithm $H_{-1,0,0,1}(1)$, cf. Section 4.2, which is given by a polynomial of MZVs up to weight w=4. I.e. we have to apply the PSLQ method over all monomials up to w=4

[TABLE]

An approximate numerical value of $H_{-1,0,0,1}(1)$ is

[TABLE]

and PSLQ yields

[TABLE]

In particular, monomials like $\ln(2),\ln^{2}(2),\ln^{3}(2),\zeta_{2},\zeta_{3}$ do not contribute here.

3.2 Hypergeometric Functions and their Generalizations

Simpler Feynman–parameter integrals have representations in terms of generalized hypergeometric functions HYPKLEIN ; HYPBAILEY ; SLATER1 and their generalizations such as Appell-, Kampe-De-Feriet- and related functions APPEL1 ; APPEL2 ; KAMPE1 ; EXTON1 ; EXTON2 ; SCHLOSSER ; Anastasiou:1999ui ; Anastasiou:1999cx ; SRIKARL ; Lauricella:1893 ; Saran:1954 ; Saran:1955 . This is due to the hyperexponential nature of the Feynman–parameter integrals, implying real exponents due to the dimensional parameter $\varepsilon$. These representations map multiple integrals to single series (for generalized hypergeometric functions) and double infinite series (e.g. for Appell series), which finally have to be solved by applying summation theory. The simplest function is Euler’s Beta-function implying the series of ${}_{p+1}F_{p}$ functions

[TABLE]

Representations of this kind are usually sufficient for massless and massive single–scale two–loop problems Hamberg:1990np ; HAMBERG ; Buza:1995ie ; Bierenbaum:2007qe . In the case of three–loop ladder graphs Appell-functions are appearing Ablinger:2012qm ; Ablinger:2015tua . There are some more classes of higher transcendental functions of this kind, which have been studied in the mathematical literature EXTON1 ; EXTON2 ; SRIKARL . The corresponding representations allow the expansion in the dimensional parameter $\varepsilon$. At a given level in the calculation of Feynman diagrams one will not find corresponding known function representations and one has to invoke other methods of integration. One way to derive analytic infinite sum representations are Mellin–Barnes integrals to which we turn now.

3.3 Analytic Solutions using Mellin–Barnes Integrals

The higher transcendental functions discussed in Section 3.2 have representations in terms of Pochhammer–Umlauf integrals POCHHAMMER ; HYPKLEIN ; KF and related to it, by Mellin–Barnes integrals BARNES1 ; MELLIN1 . They are defined by

[TABLE]

cf. e.g. SMIRNOV . Here the contour integral is understood to be either being closed to the left or the right surrounding the corresponding singularities. The Mellin–Barnes decomposition is analogous to the binomial (series) expansion for $\alpha<0$. After its application, various more Feynman parameters can be integrated using the technique described in Section 3.2. In every application the decomposition introduces a number of infinite sums of depth one according to the residue theorem. There exist some packages for Mellin–Barnes integrals Czakon:2005rk ; Smirnov:2009up ; Gluza:2007rt ; Gluza:2010rn , allowing also for numerical checks. Finally all the produced sums have to be solved using multi–summation methods. Therefore one is advised to apply this method very carefully. Not all expressions generated by this method can be analytically summed using the presently know technologies, cf. Schneider:2007a ; Schneider:2013a . Sometimes Mellin–$N$ space techniques may lead to elliptic structures, while $x$–space techniques do not, cf. Ablinger:2017xml , and sum–representations have to be cast back into definite integral representations first.

3.4 Hyperlogarithms

In a wide class of cases Feynman integrals can be represented by combinations of Kummer–Poincaré integrals KUMMER ; POINCARE ; LADAN ; CHEN ; GONCHAROV for (a part) of their expansion coefficients in $\varepsilon$. Let us assume one can isolate these terms, see vonManteuffel:2014qoa , and forms a corresponding finite multi–integral. The method of hyperlogarithms Brown:2008um has originally intended to reorganize these integrals such that one can find a sequence of integrations being linear in the Feynman parameter on hand. If this is the case the result is given in terms of Kummer–Poincaré integrals. For a corresponding implementation see Panzer:2014caa . The method has first been applied to the usual massless Feynman integrals. A generalization for massive integrals also containing local operator insertions has been given in Ablinger:2014yaa , with an implementation in Wissbrock:2015faa . Here also certain non–linear Feynman parameter structures, breaking multi–linearity, could be integrated.

3.5 The Method of Differential Equations

In single–scale processes systems of ordinary differential equations for the master integrals are naturally obtained by the IBP–relations differentiating for a parameter $x$.222Correspondingly, in the case of more parameters, partial differential equation systems are obtained. The master integrals may then be calculated by solving these systems under given physical boundary conditions, Kotikov:1990kg ; Bern:1992em ; Remiddi:1997ny ; Gehrmann:1999as . One considers the system

[TABLE]

which may also be transformed into the scalar differential equation

[TABLE]

with $p_{n}\neq 0$, and $(n-1)$ equations for the remaining solutions, which are fully determined by the solution $f_{1}(x)$. In setting up these systems one has to perform the expansion in $\varepsilon$ in parallel in the decoupling.

An important class of differential equations is formed by the first order factorizing systems, after applying the decoupling methods Zuercher:94 ; Ablinger:2013jta encoded in Oresys ORESYS , which appear as the simplest case. Eq. (19) may be transformed into Mellin space, decoupled there and solved using the efficient methods of the package Sigma, cf. Ref. Ablinger:2015tua .

The decoupled differential operator of (20) can be written in form of a combination of iterative integrals, cf. Section 4.2,

[TABLE]

with $q_{k}(x)=0$ for $1\leq k\leq m$. Further, $\gamma_{m+1}=0$ if $\bar{g}(x)=0$ in (20), and $\gamma_{m+1}=1$ and $q_{m+1}(x)$ being a mild variation of $\bar{g}(x)$ if $\bar{g}(x)\neq 0$. These solutions are d’Alembertian Abramov:94 since the master integrals appearing in quantum field theories obey differential equations with rational coefficients, the letters $h_{i}$, which constitute the iterative integrals, have to be algebraic. The solution can be computed using the package HarmonicSums Ablinger:2017Mellin . More generally, also Liouvillian solutions Singer:81 can be calculated with HarmonicSums utilizing Kovacic’s algorithm Kovacic:86 . This algorithm has been applied in many massive three–loop calculations so far, see also Ablinger:2015tua ; Ablinger:2017ptf ; Ablinger:2017hst ; Ablinger:2018zwz .

If being transformed to the associated system of difference equations, the same holds, if this system is also first order factorizing. The solution of the remaining equations are directly obtained by the first solution.

In the multi–variate case, the $\varepsilon$–representation of a linear system of partial differential equations

[TABLE]

is important, as has been recognized in Refs. Kotikov:2010gf ; Henn:2013pwa , see also Henn:2014qga . The matrices $A_{n}$ can now be transformed in the non–Abelian case by

[TABLE]

see also NOVIKOV:80 ; Sakovich:1995 , and one now intends to find a matrix $B$ to transform (23) into the form

[TABLE]

if possible. This then allows solutions in terms of iterative integrals. A formalism for the basis change to the $\varepsilon$–basis has been proposed in Lee:2014ioa and implemented in the single–variate case in Prausa:2017ltv ; Gituliar:2017vzm and in the multi–variate case in Meyer:2017joq .

3.6 The Method of Arbitrary Large Moments

In the case of single–scale problems the corresponding class of Feynman integrals depends on a real parameter $x\in[0,1]$, which is given e.g. as the ratio of two Lorentz invariants. For any power in $\varepsilon$ one would like to find the corresponding function in $x$ analytically. In a series of cases, cf. e.g. Ablinger:2014vwa ; Ablinger:2014nga ; Henn:2016tyf ; Ablinger:2018zwz , one may represent the solution in terms of a formal Taylor series in the variable $x$. The differential equations implied by the integration-by-parts method Chetyrkin:1981qh ; Laporta:2001dd ; Studerus:2009ye ; vonManteuffel:2012np ; MARSEID can now be turned into recurrences using the Taylor series ansatz. In solving the corresponding system one may generate a large number of Mellin moments for the different projections on the individual color factors and multiple zeta values Blumlein:2009cf . This is the case independently of the fact that the corresponding $x$– or $N$–space solution is given by iterative integrals or iterative–noniterative integrals. The corresponding method has been described in Ref. Blumlein:2017dxp . These moments can then be used as an input to the method described in Section 3.7 to find the associated difference equations. In some applications for single scale massive three–loop integrals Ablinger:2017ptf 8000 moments could be calculated. This is by far more than possible using standard methods like Mincer Gorishnii:1989gt , MATAD Steinhauser:2000ry or Q2E Harlander:1997zb ; Seidensticker:1999bb . Based on this number of moments, the formal power series may be used as highly precise semi–analytic numeric representations, in case the corresponding series expansion has been performed for the physical quantity to be evaluated. If analytic continuations are still necessary, the method cannot be applied directly.

3.7 Guessing One-dimensional Integrals

As has been described in Section 3.6 single–variate multiple Feynman parameter integrals can be either expanded into formal Taylor series or can be Mellin–transformed

[TABLE]

In both cases one tries now to find the associated difference equation NOERLUND to the set of moments, e.g. $\{G(2),G(4),....,G(2m)\},m\in\mathbb{N}$ Larin:1993vu ; Larin:1996wd ; Retey:2000nq ; Blumlein:2004xt . Indeed such an equation exists in many cases, as e.g. for (massive) operator matrix elements Bierenbaum:2009mv , but also for single–scale Wilson coefficients, Ref. Vermaseren:2005qc . If a suitably large number of moments has been calculated analytically the associated series of rational numbers can now be used as input for the guessing algorithm GSAGE , which is also available in Sage SAGE , exploiting the fast integer algorithms available there. The method finally returns the wanted difference equation, and tests it by a larger series of further moments. This method has been applied in Ref. Blumlein:2009tj to obtain from more than 5000 moments the massless unpolarized three–loop anomalous dimensions and Wilson coefficients in deep-inelastic scattering Vogt:2004mw ; Moch:2004pa ; Vermaseren:2005qc . Recently, the method has been applied ab initio in the calculation of three–loop splitting functions Ablinger:2017tan and the massive two– and three–loop form factor Ablinger:2018zwz ; Ablinger:2018yae . In the case of a massive operator matrix element 8000 moments Ablinger:2017ptf could be calculated and difference equations were derived for all contributing color and $\zeta$-value structures. Analytic solutions can be found using the package Sigma Schneider:2007a ; Schneider:2013a , provided the problem is solvable in difference field theory. In other cases at least the first order factorizing parts can be factored off. Other techniques are then needed to determine the remainder part of the solution.

3.8 The Almkvist–Zeilberger Algorithm

Since Feynman parameter integrals, depending on an additional parameter $x$, can be given as integrals over $\{x_{i}|_{1=1}^{n}\}\in[0,1]^{n}$, they form the multi–integral $I(x)$, depending also on $\varepsilon$. The dependence on the real parameter $x$ may be transformed into one on an integer parameter $N$, see Section 3.6. The Almkvist–Zeilberger algorithm AZ1 ; AZ2 is providing a method to find either an associated differential equation for $I(x)$ or a difference equation for $I(N)$, the coefficients of which are either polynomials in $\{x,\varepsilon\}$ or $\{N,\varepsilon\}$,

[TABLE]

[TABLE]

Both equations may be inhomogeneous, where the inhomogeneities emerge as known functions from lower order problems. An optimized and improved algorithm for the input class of Feynman integrals has been implemented in the MultiIntegrate package Ablinger:2013hcp ; Ablinger:2015tua . It can either produce homogeneous equations of the form (27,28) or equations with an inhomogeneity formed out of already known functions.

3.9 Iterative-Noniterative Integrals and Elliptic Solutions

Non–first order factorizing systems of differential or difference equations for the master integrals, cf. Section 3.5, occur at a certain order in massive Feynman diagram calculations. Well–known examples for this are the sun–rise integral, cf. e.g. Broadhurst:1993mw ; Laporta:2004rb ; Bloch:2013tra ; Adams:2013kgc ; Adams:2014vja ; Adams:2015gva ; Adams:2015ydq , the kite integral SABRY ; Remiddi:2016gno ; Adams:2016xah , the three–loop QCD–corrections to the $\rho$–parameter Ablinger:2017bjx ; Grigo:2012ji ; Blumlein:2018aeq , and the three–loop QCD corrections to the massive operator matrix element $A_{Qg}$ Ablinger:2017ptf . After separating the first–order factorizing factors a Heun differential equation HEUN remains in the case of the $\rho$–parameter. One may write the corresponding solution also using ${}_{2}F_{1}$–functions with rational argument IVH ; Ablinger:2017bjx and rational parameters. It is now interesting to see whether these solutions can be expressed in terms of complete elliptic integrals, which can be checked algorithmically using the triangle group TAKEUCHI .

In the examples mentioned one can find representations in terms of complete elliptic integrals of the first and second kind, K and E, cf. TRICOMI ; WITWAT , and the question arises whether an argument translation allows for a representation through only K. Criteria for this have been given in Herfurtner1 ; Movasati1 . In the case of the three–loop QCD-corrections to the $\rho$–parameter, however, this is not possible.

The homogeneous solution of the Heun equations are given by ${}_{2}F_{1}$–solutions $\psi_{k}^{(0)}(x),k=1,2$, at a specific rational argument. These integrals cannot be represented such that the variable $x$ just appears in the boundaries of the integral. The inhomogeneous solution reads

[TABLE]

with $N(x)$ and $W(x)$ the inhomogeneity and the Wronskian. $C_{1,2}$ are the integration constants. Through partial integration the ratio $N(x)/W(x)$ can be transformed into an iterative integral. Since $\psi_{k}^{(0)}(x)$ cannot be written as iterative integrals, $\psi(x)$ is obtained as an iterative non–iterative integral Blumlein:2016a ; Ablinger:2017bjx of the type

[TABLE]

with $r(x)$ a rational function and $F_{m}$ a non–iterative integral. Usually more than one non–iterative integral will appear in (3.9). $F_{m}$ denotes any non–iterative integral, implying a very general representation, cf. Ablinger:2017bjx .333This representation has been used in a more special form also in Remiddi:2017har later. In Ref. Adams:2018yfj an $\varepsilon$–form for the Feynman diagrams of elliptic cases has been found recently. However, transcendental letters contribute here. This is in accordance with our earlier finding, Eq. (3.9), which, as well is an iterative integral over all objects between the individual iterations and to which now also the non–iterative higher transcendental functions $F_{m}[r(y_{m})]$ contribute. One may obtain fast convergent representations of $\mathbb{H}(x)$ by overlapping series expansions around $x=x_{0}$ outside possible singularities, see Ref. Ablinger:2017bjx for details.

Let us return to the elliptic case now. Here one one may transform the kinematic variable $x$ occurring as ${\rm\bf K}(k^{2})={\rm\bf K}(r(x))$ into the variable $q=\exp[i\pi\tau]$ analytically with

[TABLE]

by applying a 3rd order Legendre–Jacobi transformation, where $\vartheta_{l},l=1,...,4$ denote Jacobi’s $\vartheta$-functions and ${\sf Im}(\tau)>0$. In this way Eq. (29) is rewritten in terms of the new variable. The integrands are given by products of meromorphic modular forms, cf. SERRE ; COHST ; ONO1 , which can be written as a linear combination of ratios of Dedekind’s $\eta$-function

[TABLE]

Depending on the largest multiplier $k\in\mathbb{N}$, $k_{m}$, of $\tau$ in the argument of the $\eta$-function, the solution transforms under the congruence subgroup $\Gamma_{0}(k_{m})$. One can perform Fourier expansions in $q$ around the different cusps of the problem, cf. ZUDILIN ; BROADH18 .

In the case that the occurring modular forms are holomorphic, one obtains representations in Eisenstein series with character, while in the meromorphic case additional $\eta$–factors in the denominators are present. In the former case the $q$–integrands can be written in terms of elliptic polylogarithms in the representation Adams:2014vja ; Adams:2015gva

[TABLE]

and products thereof, cf. Adams:2015gva . The corresponding $q$–integrals can be directly performed. The solution (29) usually appears for single master integrals. Other master integrals are obtained integrating further other letters, so that finally representations by $\mathbb{H}(x)$ occur. Iterated modular forms, resp. Eisenstein series, have been also discussed recently in Adams:2017ejb ; Broedel:2018iwv . Efficient numerical calculations of modular forms based on $q$–series were obtained in Bogner:2017vim .

For systems which factorize only to 3rd and higher order much less is known.

3.10 Iterative Integrals of Functions with More Variables

The occurrence of several masses or additional external non–factorizing scales in higher order loop- and phase–space integrals leads in general to rational and root–valued letters with real parameter letters in the contributing alphabet, cf. Ablinger:2017xml ; Ablinger:2018brx ; Blumlein:2019srk ; Blumlein:2019qze ; Blumlein:2019zux . In the case of the loop integrals one obtains letters of the kind

[TABLE]

The iterative integrals and constants which appeared in Ablinger:2017xml ; Ablinger:2018brx could finally be all integrated to harmonic polylogarithms containing complicated arguments, at least up to one remaining integration, which allows their straightforward numerical evaluation.

In the case of phase space integrals with more scales, e.g. Blumlein:2019srk ; Blumlein:2019qze , also letters contribute, which may imply incomplete elliptic integrals and iterated structures thereof. Contrary to the functions obtained in Section 3.9 these are still iterative integrals, because the boundaries of the phase–space integrals are real parameters and not constants. The integrands could not by rationalized completely by variable transformations, see also Besier:2018jen . Contributing letters are e.g.

[TABLE]

with $k,z\in[0,1]$. The corresponding iterative integrals are called Kummer–elliptic integrals. They are derived using the techniques described in Refs. Ablinger:2014bra ; RAAB1 ; GuoRegensburgerRosenkranz .

4 A Series of Function Spaces

Intermediary and final results for zero– and single–scale multi–loop calculations have representations by special functions as polynomials over $\mathbb{Q}$. In the case of of zero–scale quantities these are special numbers. For single scale quantities one either uses finite nested sum representations in Mellin $N$–space or iterative integral representations in $x$–space. Here $x$ denotes a Lorentz invariant ratio of two physical quantities. Both spaces are related to each other by the Mellin transform (26), where $f(x)$ denotes an iterative integral. The zero–scale quantities can be obtained e.g. in the limit $N\rightarrow\infty$ of these Mellin transforms or by the values $f(x=1)$.

4.1 Classes of Nested Sums

The methods described in Section 3 very often lead to finite nested sum representations for which algorithms exist Schneider:2007a ; Schneider:2013a to cast these sums into indefinitely nested sums. They are given by

[TABLE]

with $\bar{\mathfrak{A}}$ the associated alphabet of functions. The sums obey quasi–shuffle relations HOFFMAN ; Blumlein:2003gb . The simplest structures are the finite harmonic sums Vermaseren:1998uu ; Blumlein:1998if , where $g_{b}(k)=({\rm sign}(b))^{k}/k^{|b|},~{}~{}b\in\mathbb{N}\backslash\{0\}$. A generalization is obtained in the cyclotomic case Ablinger:2011te . Here the characteristic summands are $g_{a,b,c}(k)=(\pm 1)^{k}/(ak+b)^{c}$, with $a,b,c\in\mathbb{N}\backslash\{0\}$. Further, the generalized harmonic sums have letters of the type $b^{k}/k^{c}$, with $c\in\mathbb{N}\backslash\{0\},b\neq 0,b\in\mathbb{R}$, Ablinger:2013cf . Another generalization are nested finite binomial and inverse-binomial sums, containing also other sums discussed before. An example is given by

[TABLE]

with $c_{1}\approx 0.10184720\dots$, cf. Ablinger:2014bra . Here the indices ${\sf w_{k}}$ label specific letters given in Ablinger:2014bra . Infinite binomial and inverse binomial sums have been considered in Davydychev:2003mv ; Weinzierl:2004bn . Given the general structure of (36) many more iterated sums can be envisaged and may still appear in even higher order calculations.

4.2 Classes of Iterated Integrals

Iterated integrals have the structure

[TABLE]

where $f_{c}$ are real functions and are the letters of the alphabet $\mathfrak{A}$. Iterated integrals obey shuffle relations REUTENAUER ; Blumlein:2003gb which allows to represent them over a multinomial basis of fewer terms.

The simplest iterative integral having been considered in quantum field theory are the Nielsen integrals for the two–letter alphabets $\{1/x,1/(1-x)\}$ or $\{1/x,1/(1+x)\}$ NIELSEN ; KMR70 ; Kolbig:1983qt ; Devoto:1983tc , covering also the polylogarithms LEWIN1 ; LEWIN2 ; Devoto:1983tc . This class has later been extended to the harmonic polylogarithms Remiddi:1999ew build over the alphabet $\{1/x,1/(1-x),1/(1+x)\}$. A further extension is to the real representations of the cyclotomic polylogarithms, with $\{1/x,1/\Phi_{k}(x)\}$ Ablinger:2011te , where $\Phi_{k}(x)$ denotes the $k$th cyclotomic polynomial. Another extension is given by Kummer–Poincaré iterative integrals over the alphabet $\{1/(x-a_{i}),~{}~{}a_{i}\in\mathbb{C}\}$, KUMMER ; POINCARE ; LADAN ; CHEN ; GONCHAROV . Properties of these functions have been studied in Refs. Moch:2001zr ; Ablinger:2013cf . In general one may have also more general denominator polynomials $P(x)$, which one can factor into

[TABLE]

in real representations. One then performs partial fractioning for $1/P(x)$ and forms iterative integrals out of the obtained letters. Further classes are found for square–root valued letters as studied e.g. in Ref. Ablinger:2014bra . In multi–scale problems, cf. e.g. Ablinger:2017xml ; Blumlein:2019srk ; Blumlein:2019qze ; Blumlein:2019zux and Section 3.10, further root–valued letters appear, like also the Kummer–elliptic integrals Blumlein:2019qze .

4.3 Classes of Associated Special Numbers

For the sums of Section 4.1 which are convergent in the limit $N\rightarrow\infty$ and the iterated integrals of Section 4.2 which can be evaluated at $x=1$ one obtains two sets of special numbers. They span the solution spaces for zero–scale quantities and appear as boundary values for single–scale problems. Examples for these special numbers are the multiple zeta values Blumlein:2009cf , associated to the harmonic sums and harmonic polylogarithms, special generalized numbers Ablinger:2013cf like ${\rm Li}_{2}(1/3)$, associated to generalized sums and to Kummer–Poincaré iterated integrals, special cyclotomic numbers Ablinger:2011te like Catalan’s number, special binomial numbers Ablinger:2014bra , as e.g. arccot($\sqrt{7}$), and special constants in the elliptic case Ablinger:2017bjx ; Laporta:2017okg . The latter numbers are given by integrals involving complete elliptic integrals at special rational arguments and related functions. In general these numbers obey more relations than the finite sums and iterated integrals. One may use the PSLQ–method to get a first information on relations between these numbers occurring in a given problem and proof the conjectured relations afterwards.

4.4 Numerical Representations

Physical observables based on single scale quantities can either be represented in Mellin $N$–space or $x$–space. Representations in Mellin $N$–space allow the exact analytic solution of evolution equations Blumlein:1997em and scheme-invariant evolution equations can be derived in this way Blumlein:2000wh ; Blumlein:2004xs . The $x$–space representation is then obtained by a single numerical integral around the singularities of the respective quantity for $N\in\mathbb{C}$, cf. Blumlein:1997em , requiring to know the complex representation of the integrand in $N$–space. In the case of harmonic sums semi–numerical representations were given in Blumlein:2000hw ; Blumlein:2005jg . Furthermore, it is known that basic harmonic sums, except of $S_{1}(N)$, which is represented by the Digamma function, and its polynomials, have a representation by factorial series FACT1 ; FACT2 , which has been used in Blumlein:2009fz ; Blumlein:2009ta for their asymptotic representation, see also Kotikov:2005gr . One uses then the recursion relations, which can be obtained from (36), to move $N\in\mathbb{C}$ from the asymptotic region to the desired point on the integration contour in the analyticity region of the problem. This can be done for the sums of the type being described in Refs. Ablinger:2011te ; Ablinger:2013cf ; Ablinger:2014bra as well, since also in this case asymptotic expansions can be provided, at least for certain combinations of sums occurring in the respective physical problem, cf. Ablinger:2014yaa ; Ablinger:2015tua . In the case that the corresponding relations are not given in tabulated form, they can be calculated using the package HarmonicSums Vermaseren:1998uu ; Blumlein:1998if ; Ablinger:2014rba ; Ablinger:2010kw ; Ablinger:2013hcp ; Ablinger:2011te ; Ablinger:2013cf ; Ablinger:2014bra ; Ablinger:2017Mellin . Relations for harmonic sums are also implemented in summer Vermaseren:1998uu , and for generalized harmonic sums in nestedsums Weinzierl:2002hv , Xsummer Moch:2005uc , and PolyLogTools Duhr:2019tlz .

In other applications one may want to work in $x$–space directly. Here numerical representations are available for the Nielsen integrals KMR70 , the harmonic polylogarithms Gehrmann:2001pz ; Vollinga:2004sn ; Maitre:2005uu ; Ablinger:2017tqs ; Ablinger:2018sat , the Kummer–Poincaré iterative integrals Vollinga:2004sn , and the cyclotomic harmonic polylogarithms Ablinger:2018zwz . These representations are also useful to lower the number of numerical integrations for more general problems, e.g. in the multi–variate case. The relations for the corresponding quantities are implemented for the harmonic polylogarithms in Remiddi:1999ew ; Maitre:2005uu and for all iterative integrals mentioned, including general iterative integrals, in the package HarmonicSums.

5 Conclusions

In parallel to the analytic higher–loop calculations in Quantum Field Theory the associated mathematical methods have been developed by theoretical physicists and mathematicians since the 1950ies. We witness a very fast development since the late 1990ies approaching difficult massive problems at two–loop and higher order and massless problems form three loops onward. The classical methods of polylogarithms and Nielsen–integrals which were standard means, turned out to be not sufficient anymore. Since then more and more special number– and function spaces have been revealed, studied and were brought to flexible practical use in very many applications. Moreover, a wide host of analytic integration and summation methods has been developed during a very short period. In this way very large physics problems could be solved analytically – a triumph of the exact sciences, also thanks to various groundbreaking methods in computer algebra. In this context the goal is to improve the accuracy of the fundamental parameters of the Standard Model of the elementary particles further. Within the present projects this concerns in particular the relative precision of the strong coupling constant $\alpha_{s}(M_{Z}^{2})$ to less than 1% and of the $\overline{\rm MS}$ mass of the charm quark to better than 1.5 %.

At even higher loop order and for more separated final state legs, introducing more masses and kinematic invariants, one expects further mathematical structures to contribute. Possible structures of this kind could be Abel–integrals NEUMANN and integrals related to K3–surfaces BSCH . More inclusive methods, like the method of differential equations, can certainly determine the degree of non–factorization of a physical problem. However, one would like to know in a closer sense the respective analytic solution. Here cutting methods can be of use since the underlying integrands can be systematically related to the final integral by (various) Hilbert-transform Hilbert:1912 ; KRONIG ; KRAMERS .444For a recent application to the one–loop case, see e.g. Abreu:2017mtm . In this way integrand structures are revealed, which are somewhat hidden in the case of differential equations. This method has been advocated early by M. Veltman Veltman:1963th , see also REMI1 .

This process to master highly complex Feynman integrals using analytic methods is of course just at the beginning and will develop further given the present and future challenges in the field. All of these results put experimental analyses in precision measurements at the high energy colliders into the position to analyze the data with much reduced theory errors and we will get far closer in our insight into the structure of the micro cosmos to reveal its ultimate laws. The interdisciplinary joined effort by mathematicians, theoretical and experimental particle physicists and experts in computer algebra makes this possible and allows to answer quite a series of fundamental scientific questions of our time.

I would like to give my warmest thanks to Peter Paule for his continuous collaboration and support to the DESY–RISC collaboration, starting with our first contacts in 2005, arranged by Bruno Buchberger. This scientific symbiosis has produced a large number of methods to tackle quite a series of difficult problems since, and is continuing to do so in the future. Physics, mathematics, and computer algebra profit from this and reach new horizons, which, not at all, could have been imagined. In this way we follow together the motto D. Hilbert has given to us:

Wir müssen wissen. Wir werden wissen.

Le but unique de la science, c’est l’honneuer de l’esprit humain.555Jacobi to Legendre, July 2nd, 1830.

Acknowledgment. I would like to thank J. Ablinger, D. Broadhurst, A. De Freitas, D. Kreimer, A. von Manteuffel, P. Marquard, S.-O. Moch, W.L. van Neerven, P. Paule, C. Schneider, K. Schönwald, J. Vermaseren and S. Weinzierl for countless fruitful discussions. This work was supported by the EU TMR network SAGEX Marie Skłodowska-Curie grant agreement No. 764850 and COST action CA16201: Unraveling new physics at the LHC through the precision frontier.