The SAGEX Review on Scattering Amplitudes, Chapter 4: Multi-loop Feynman Integrals

The analytic integration and simplification of multi-loop Feynman integrals to special functions and constants plays an important role to perform higher order perturbative calculations in the Standard Model of elementary particles. In this survey article the most recent and relevant computer algebra and special function algorithms are presented that are currently used or that may play an important role to perform such challenging precision calculations in the future. They are discussed in the context of analytic zero, single and double scale calculations in the Quantum Field Theories of the Standard Model and effective field theories, also with classical applications. These calculations play a central role in the analysis of precision measurements at present and future colliders to obtain ultimate information for fundamental physics.


Introduction
The present and upcoming high luminosity results at the Large Hadron Collider (LHC) at CERN, with input of the results measured at the ep-collider HERA at DESY, yield a big amount of precision data which require further fundamental precision calculations in perturbative Quantum Chromodynamics (QCD). This also applies to projects in the Future like the EIC [1], the LHeC [2,3], the ILC [4][5][6][7] or CLIC [8][9][10], the FCC ee [11], and the proton version of the FCC [11], also for Quantum Electrodynamics (QED). Many of these outstanding problems can be formulated by large expressions in terms of hundred thousands (and even millions) of sophisticated Feynman integrals at higher loop order of presently up to three and four loops, with no, a single and two scales, or even multiscale problems. † Up to one scale, and in much fewer cases at two scales, technologies have been designed to calculate these integrals analytically over sets of basic functions, the properties of which have been studied to a certain extent. Furthermore, numerical representations of these building blocks have been derived.
The central objects are s-fold multiple integrals of the form F (n, ε) = where the discrete parameter n stands for the Mellin variable, and ε = D − 4, ε ∈ R, |ε| 1, is the dimensional parameter. A crucial property is that the integrand f is hyperexponential ‡ in each of the integration variables x i (1 ≤ i ≤ s) and hypergeometric in the discrete parameter n. In particular, one is interested in calculating the first coefficients of their Laurent series expansion w.r.t. ε: F (n, ε) = F l (n)ε l + F l+1 (n)ε l+1 + · · · + F r (n)ε r + O(ε r+1 ).
(2) x n−1 f (x, ε)dx (3) or its power series representation During the last decades more and more significant methods have been derived to simplify such Feynman integrals. Based on the representation (for instance (2)-(4)), we will present important tools that are currently used to perform such challenging calculations and discuss the associated special function spaces. We further relate these aspects to different precision calculations. It is needless to say that we had to leave out the description of a series of techniques, which are also important. This concerns a series of aspects, which have been surveyed in Ref. [12], appearing in the same volume, and has been agreed between the different authors. It concerns e.g. the use of the symbol [13] and specific Hopf algebra structures [14][15][16], which are omnipresent in quantum field theoretic calculations. Related to this, many methods found in algebraic and arithmetic geometry are applicable [17][18][19][20][21][22][23][24]. For mostly numerical methods in use in multi-leg calculations we refer to [25]. †For zero scale problems results are available also at the five loop level.
‡h(x) is hyperexponential (or hypergeometric) in ) is a rational function in x.
One way to classify the emergence of new mathematical structures in quantum field theories is given by the study of their differential equations. This is first of all a practical issue at the respective loop-level, where these structures are recognized and tied up with the respective graph topologies. There is in general no all order statement possible ab initio. However, as has been found out recently, Calabi-Yau motives play an important role here, cf. e.g. [24]. A selection criterion on what we are focusing on in the following are key technologies for multi-loop calculations in the massive case, presently to threeloop order in the zero-, single-, and two-scale cases. These technologies do synonymously apply to the corresponding massless calculations. There, clearly simplifications can be obtained using even other technologies, unlike the case in the massive case.
In the following we discuss the following key research topics: • guessing methods (see Section 2), • solving linear recurrences and differential equations (see Section 3), • solving coupled systems (see Section 4), • transformation to special integral and sum representations (see Section 5), • symbolic summation (see Section 6), • symbolic integration (see Section 7), • the large moment method (see Section 8), • special function tools (see Section 9) • and concrete calculations in the Quantum Field Theories of the Standard Model and within effective field theories (see Section 10).
We emphasize that each of the different techniques cannot be considered as a standalone toolbox. Contrary, they all have to be applied in non-trivial interactions. In particular, based on a concrete problem, one has to choose the best tactic among the conglomeration of tools. For further and supplementary aspects we refer also to [26][27][28] and Chapter 3 [12] of the SAGEX review. We will conclude this survey on multi-loop tools and multi-loop calculations in Section 11.

Guessing methods
Often physical quantities can be evaluated up to a certain precision and one seeks for a mathematical representation that allows one to represent the data in a more compact fashion, to gain further insight and to support further calculations that depend on these quantities. Here we will emphasize two crucial tactics: (1) to predict from a given floating point number (that approximates a real number to very high precision) an alternative representation in terms of special constants and (2) to guess from a finite set of evaluations at integer points a linear recurrence or linear differential equation that satisfies all evaluations at integer points of the physical quantity.

Guessing integer relations
Using the LLL algorithm [29], or the PSLQ algorithm introduced in [30] and substantially improved in [31,32], one can try to solve the following problem: Given a finite set of finite floating point numbers a 1 , . . . , a n with high precision (say l fractional digits), find integers z 1 , . . . , z n ∈ Z as small as possible in its absolute value § such that z 1 a 1 + · · · + z n a n < 10 −m , where m is large.
In order to obtain further confidence one may apply this method with further precision (i.e., more digits l of the input) and checks if the obtained result remains the same but m gets larger. For instance, suppose that we are given the finite floating point number a 1 = 5.68700407989058207630312605688168433418849655155997 with l = 50 fractional digits that approximates a real number r 1 ∈ R. Then one can use, e.g., the Mathematica implementation of the PSLQ algorithm to search for an alternative representation in terms of ζ(2), ζ(3) and ζ (5). Namely, by activating FindIntegerNullVector[{a 1 , N[Zeta [2], 50], N[Zeta [3], 50], N[Zeta [5], 50]}] one obtains the result (z 1 , z 2 , z 3 , z 4 ) = (−2, 2, 5, 2) with m = 49. Thus we may conjecture that holds. We remark that the PSLQ method finds this relations already with l = 5 fractional digits with precision m = 4. Of course, if one is given even more digits of r 1 , one may check if even more digits m agree. Similarly, one may apply PSLQ again for this improved data in order to check if there is even a smaller relation (with smaller z i ). Summarizing, the method can be very efficiently used if one knows the set of numbers, by which the result of a calculation is finally spanned or if one wants to find a linear combination of such numbers. A highly non-trivial example is, e.g., the calculation of the 5-loop β-function in QCD in Ref. [33], where these guessing tools were instrumental. For a recent survey on these techniques (covering not only PSLQ but also the LLL approach) and further applications we refer to [34].

Guessing recurrences and differential equations
In Section 8 below we will introduce a method that enables one to compute many moments F (n, ε) in (3) or coefficients in (4) for n = 0, 1, 2, . . . . More precisely, if we §Since any finite floating point number can be written as a rational number, this problem can be always solved if the integers z i can be arbitrarily large. Thus a solution to the problem might indicate a proper relation among the approximated real numbers if l is large but the values z i are small.
Within multi-loop calculations these moments depend linearly also on special constants, such as the multiple zeta values [35], with rational coefficients. This finally leads to several finite sequences, F (0), F (1), . . . , F (µ), of rational numbers. Then given these numbers, one can try to guess a linear recurrence a 0 (n)F (n) + a 1 (n)F (n + 1) + · · · + a λ (n)F (n + λ) = 0 (5) of order λ with polynomial coefficients a i (n) ∈ Q[n] that contains this finite sequence as solution. Namely, fixing the order λ and assuming that the degrees of the polynomials a i (n) are less than or equal to δ, one searches for the r = (δ + 1)(λ + 1) unknown coefficients. More precisely, by setting n = 0, . . . , r − 2 in (5) and plugging in the rational numbers F (0), . . . , F (r + λ − 2) one gets r − 1 equations in r unknowns over the rational numbers which can be solved by linear algebra. In many cases this yields solutions that do not hold for n ≥ r − 1. Thus one usually takes an over-determined system (by more evaluations, say 0 ≤ n ≤ r + 100). In this way one can exclude basically all wrong solutions. Finally, given a found solution one usually checks at many extra points if the recurrence is still valid. This gives further evidence that the guessed recurrence is reliable. This tactic implemented, e.g., in the Maple package gfun [36] or the Mathematica package GeneratingFunction [37] is surprisingly simple and can be carried out in this naive fashion for small examples (i.e., for recurrences of small orders λ and small degree bounds δ). For large examples within QCD calculations this straightforward procedure utterly fails. Here highly efficient computer algebra technologies, such as homomorphic image calculations and rational/polynomial reconstructions, are essential [38]. Using in addition gcd-calculations to determine recurrences with minimal order, the Mathematica implementation Guess.m by Kauers could be utilized with about µ = 5000 moments to guess all the recurrences that determine the massless unpolarized 3-loop anomalous dimensions and Wilson coefficients in deep-inelastic scattering [39][40][41] in Ref. [42], see also [43][44][45][46][47]. For even larger problems, the highly efficient Sage implementation in ore algebra [48] (utilizing among other smart techniques the fast integer arithmetic of Flint) was instrumental to guess linear recurrences with minimal order. E.g., for the massive form factor [49,50] about µ = 10000 moments were needed to obtain recurrences up to order λ = 55 and degree δ = 1300. In the case of a massive operator matrix element 8000 moments [51] could be calculated and difference equations were derived for all contributing color and ζ-value structures. Recently, also the 3loop splitting functions [44], the anomalous dimensions from off shell operator matrix elements [44,52,45,46], and lately the polarized transition matrix element A gq (N ) [53] and the logarithmic contributions to the polarized O(α 3 s ) asymptotic massive Wilson coefficients [54] have been derived by guessing the underlying recurrence relations.
Further we note that one can guess in a similar fashion a linear differential equation x. Both, the Mathematica implementation in Guess.m and the Sage implementation in [48] cover this extra feature.
Given such recurrences, one succeeds in many cases to solve the recurrences in terms of special functions that are most relevant in QCD calculations. Further details on these solving aspects will be given in the next section.

Solving linear recurrences and differential equations
As already motivated in Section 2.2 above and further emphasized in Sections 6-8, one can derive a linear recurrence (linear difference equation) or a linear differential equation which contains the given multi-loop Feynman integral (1) or a given physical expression in terms of such Feynman integrals as a solution. Then a natural strategy is to apply the available toolboxes to compute all solutions of the derived equations that can be represented in terms of certain classes of function spaces that will be introduced in more detail in Section 9. In the case that one finds sufficiently many (linearly independent) solutions one may obtain an alternative representation of the physical problem in terms of these solutions.
In the following we describe different algorithms that can provide solutions of linear difference and differential equations that occur in QCD calculations.

Ordinary linear equations
We start with equations in one variable, i.e., with ordinary linear difference equations of the form and ordinary linear differential equations of the form where D x = d dx denotes the differentiation w.r.t. x.
3.1.1. Ordinary linear difference equations The first major contribution for recurrence solving is elaborated in [55] and is substantially improved in [56]. Given rational functions a 0 (n), . . . , a λ (n), r(n) ∈ K(n) (K denotes a computable field that contains the rational numbers) it finds all rational solutions F (n) ∈ K(x) of (6). More generally, using the algorithms from [57] and the more efficient version given in [58] one can compute all hypergeometric solutions of (6), this means one can compute all solutions that can be written as hypergeometric products where l is an integer and f (k) is a rational function in k; here l is chosen such that the evaluation f (k) for k ∈ N with k ≥ l has no pole and is nonzero. In particular, the solutions can be described in terms of a product of Γ-functions, Pochhammer-symbols, factorials, binomial coefficients and rational functions. Even more generally, using the algorithms described in [59,60] one can search for all d'Alembertian solutions, i.e., all solutions that can be expressed in terms of iterative sums defined over hypergeometric products. Special cases of this class of sums are harmonic sums [61,62], cyclotomic harmonic sums [63], generalized harmonic sums [64,65] and finite binomial sums [66]; infinite binomial sums have been also studied in [67,68,66]. Further details and extra properties of such sums are presented in Section 9.
Finally, one can search in addition for all Liouvillian solutions [69] which cover in addition the interlacing of expressions in terms of iterated sums over hypergeometric products. Basically all these tools have been generalized to the setting of difference fields [70] and rings [71] (utilizing results from above and [72][73][74][75][76]) that allows one to find such solutions for difference equations (6) where the coefficients a i (n) and the inhomogeneous part r(n) are not just rational functions but can be built again by indefinite nested sums over hypergeometric products. E.g., using the summation package Sigma [77][78][79], that contains this general toolbox, one can compute for the recurrence 1 + S 1 (n) + nS 1 (n) 2 3 + 2n + 2S 1 (n) + 3nS 1 (n) + n 2 S 1 (n) 2 F (n) −(1 + n)(3 + 2n)S 1 (n) 3 + 2n + 2S 1 (n) + 3nS 1 (n) + n 2 S 1 (n) 2 F (n + 1) +(1 + n) 2 (2 + n) 3 S 1 (n) 1 + S 1 (n) + nS 1 (n) F (n + 2) = 0 the complete solution set here S 1 (n) = n k=1 1 k denotes the nth harmonic number. Internally, the recurrence operator is factorized as much as possible into linear factors. Then each extra factor provides one extra linearly independent solution which is constructed by one extra indefinite sum. In other words, finding ν linear factors (ideally ν = λ) yields ν linearly independent solutions where the most complicated solution is built by an iterative nested sum over hypergeometric products of nesting depth ν − 1; the particular solution will lead to a nested sum of depth ν. Then a key task is to simplify these sum solutions further such that the nesting depth is minimal; further aspects on such simplifications will be given in Section 6.1. We note that all solutions of a linear recurrence can be given in terms of d'Alembertian solutions if the linear recurrence operator factors completely into first-order linear factors. We call such a recurrence also first order factorizing. If this is not the case, i.e., if only parts of the recurrence can be factored into linear righthand factors then it is called non-first order factorizing. The recurrences coming from QCD calculations usually have polynomial coefficients a i (n) and the right-hand side r(n) is either 0 or is built by indefinite nested sums over hypergeometric products. One of the largest homogeneous recurrences (r(n) = 0) that have been solved with Sigma were of order λ = 35 and the degree of the coefficients of a i (n) was up to 1000 and the occurring integers required up to 1400 decimals digits; for details see, e.g., [80,42]. The largest inhomogeneous recurrences were up to order λ = 12 where r(n) may be built up to hundreds of highly nested indefinite nested sums.
In most cases Feynman diagrams or physical expressions in terms of such integrals depend on the dimensional parameter ε. In particular, this parameter ε occurs in the coefficients a i (n) and the inhomogeneous part r(n) of the recurrence (6). In some special cases, the solution F (n) can be given in terms of indefinite nested sums over hypergeometric products where ε occurs inside of the sums and products. In such situations, the above methods implemented in Sigma can find the complete solution in n and ε. However, in most instances such a closed form solution does not exist and one seeks for closed form solutions of the first coefficients F i (n) (free of ε) of the εexpansion (2). In order to accomplish this task, one can apply the algorithm from [81] implemented in Sigma in order to constructively decide if the coefficients F i (n) can be represented in terms of nested sums over hypergeometric products.
The more complicated multi-loop Feynman integrals are considered, the more complicated function spaces arise. Thus further techniques are extremely desirable that extend the class of indefinite nested sums over hypergeometric products. In this regard, one should mention the special case of factorial series [82][83][84] solutions of the form f (n) = ∞ k=0 k! (n+k)! a k . Namely, given a linear recurrence in f (n), an operator method is described in [85] and further considered in [86], to provide a linear recurrence for the sequence a k . Precisely here one can utilize the recurrence solver of Sigma to decide, if a k can be written in terms of d'Alembertian solutions. E.g., for the recurrence Petkovšek proposed new ideas in [87] to find solutions of truncated binomial sums: instead of k! (n+k)! one can choose certain products of binomial coefficients and the upper bound should be integer-linear in n.
3.1.2. Ordinary linear differential equations In various instances one is interested in a power series solution f (x) = ∞ n=0 F (n)x n of a linear differential equation (7). If the coefficients a i (x) are rational functions in x and the inhomogeneous part r(x) itself can be given in form of a power series representation, one can utilize holonomic closure properties [36,37,88] as follows. Plugging the power series ansatz into the differential equation and comparing coefficients w.r.t. x n yield a linear recurrence of the form (7) (with updated a i (n) and r(n)) for the desired coefficients F (n). In a nutshell, one can activate the available recurrence solver introduced in Section 3.1.1 to compute closed form representations of the coefficients F (n).
Alternatively, there are also direct algorithms available, similarly to the difference equation case, that can solve linear differential equations in terms of rather general classes of special functions. Namely, using the algorithms from [55] one can find all rational solutions. More generally, using [89] and, e.g., the improved versions given in [56,90] one can find all hyperexponential functions f (x). In general, the functions can be given in the form e x l h(x)dx for some rational function h and lower bound l; special cases are, e.g., rational functions or roots over such functions. More generally, one can use these algorithms to compute all d'Alembertian [59,60], i.e., all solutions that can be given in terms of iterated integrals over hyperexponential functions. Special cases of these integrals, are harmonic polylogarithms [91], cyclotomic polylogarithms [63], generalized multiple polylogarithms [64,65] but also root-valued nested integrals [66]; further details are given in Section 9. As for the recurrence case the corresponding differential operator is factorized as much as possible into linear factors.
Then each factor yields one extra linearly independent solution by introducing one extra indefinite integration quantifier. These d'Alembertian solutions can be computed with the package HarmonicSums [92]. Similarly to the recurrence case we note that all solutions of a linear differential equation can be given in terms of d'Alembertian solutions if the differential operator factors completely into first-order linear factors. We call such a differential equation also first order factorizing. If this is not the case, i.e., if only parts of the linear differential equation can be factored into first-order linear right-factors then it is called non-first order factorizing. More generally, also Liouvillian solutions [93] can be calculated partially with HarmonicSums by utilizing Kovacic's algorithm [94]. For instance, given HarmonicSums finds the general solution where √ 1 + x is hyperexponential and 3 1 − √ 1 + x is algebraic over a field generated by x and √ 1 + x. More generally, in [93] an algorithm has been described that finds all Liouvillian solutions of a homogeneous linear differential equations, i.e., all solutions that can be given by iterated integrals over hyperexponential function and functions that are algebraic over the extension below. More generally, an algorithm has been proposed in [95] that can find all Liouvillian solutions of linear differential equations whose coefficients are given in terms of functions that are Liouvillian. In some sense, this highly general solver can be considered as the continuous version of the recurrence solver [70] implemented within the package Sigma.
As already emphasized in Section 3.1.1, also the dimensional parameter ε appears in the coefficients a i (x) and the inhomogeneous part r(x) of the linear differential equation (7) when one deals with Feynman integrals. In some special cases one can use the above algorithms directly where ε may arise inside of d'Alembertian and Liouvillian solutions. However, similarly to the recurrence case, this approach usually does not work and one aims at finding closed forms of the first coefficients of the ε-expansion In this regard, the package HarmonicSums can decide constructively if the first coefficients f i (x) (free of ε) can be given in terms iterated integrals over hyperexponential functions; for the underlying algorithm we refer to [96] which is based on ideas given in [81].
By looking at more and more complicated Feynman integrals also the class of Liouvillian solutions is not sufficient. For second order linear differential equations van Hoeij proposed algorithms in [97] that can find hypergeometric series solutions ( p F q 's) in terms of certain rational function arguments. These advanced tools turned out to be instrumental to deal with the ρ-parameter in [98,99] and related quantities.

Partial linear equations
Solving partial linear difference and differential equations is a hard problem. It has been shown in [100] based on [101] that already the task to solve such equations in terms of polynomial solutions is an unsolvable problem. Recently, new methods have been introduced in [102,103] that enable one to search for (not necessarily all) rational solutions of partial linear difference equations of the form (s 1 ,...,sr)∈S a (s 1 ,...,sr) (n 1 , . . . , n r )F (n 1 + s 1 , . . . , n r + s r ) = 0, where the coefficients a (s 1 ,...,sr) are rational functions in the variables n 1 , . . . , n r and S ⊂ Z r is a finite set. In [104] further ideas coming from Section 3.1.1 have been incorporated to hunt also for solutions in terms of a given set of nested sums. E.g., suppose that we are given the partial linear difference equation − (n + 1) 2 k + n 2 + 2 4k 2 − 3kn 2 + 5kn + 12k − 2n 3 − 2n 2 + 8n + 8 F (n, k + 1) + (n + 1) 2 k + n 2 + 3 2k 2 − 2kn 2 + 2kn + 6k − n 3 − n 2 + 4n + 4 F (n, k + 2) + kn 2 (n + 2) 2 k + n 2 + 2n + 3 F (n + 1, k + 1) = 0 and the set W = {S 1 (k), S 1 (n + k), S 2,1 (n + k)} in terms of the harmonic numbers and the harmonic sum S 2,1 (n) = n k=1 S 1 (k) k 2 ; compare Section 9. Then fixing the total degree bound 5 or the arising objects in the numerator, one can compute with the package SolvePLDE introduced in [104] the 37 solutions p (1+n) 2 (1+k+n 2 ) where p is taken from the set In particular, the method for scalar linear difference equations in [80] has been carried over in this new package to search also for closed form solutions of the first coefficients of an ε-expansion.
We emphasize that this new package enables one also to attack partial linear differential equations and to find solutions in its multivariate power series expansion f (x 1 , . . . , x r ) = (n 1 ,...,nr)∈N r F (n 1 , . . . , n r )x n 1 1 . . . x nr r . Namely by plugging the power series ansatz into the partial differential equation and comparing coefficients w.r.t. x n 1 1 . . . x nr r produce a partial linear difference equation of the form (8). Thus one can utilize the tools described above to search for closed form representations of F (n 1 , . . . , n r ). In Section 4.2 this tactic will be refined further to find solutions for certain classes of coupled systems of partial linear differential equations.

Solving coupled systems of linear differential equations
In order to solve open problems at the forefront of elementary particle physics, millions of complicated Feynman integrals have to be tackled. As a preprocessing step one often applies integration-by-parts (IBP) methods [105][106][107] ¶ that crunch these integrals to a few hundred (or thousand) so-called master integrals; for a recent survey, possible refinements and applications see, e.g., [107,110,111]. Then the main task is to simplify only these master integrals to expressions in terms of special functions and to assemble the original problem with these sub-results. Most of these master integrals f i (x, ε) can be determined as solutions of coupled systems of linear differential equations. For single-variate systems they are of the form . . .
with A being a λ × λ matrix with entries from K(x, ε) where the right-hand sides are given in terms of simpler master integrals. They are either determined by other coupled ¶Here the method of syzygies [108,109] from computational algebraic geometry helps to reduce the number of contributing scalar products. systems or have to be tackled by tools presented, e.g., in Sections 6 and 7. Here we elaborate the most relevant approaches. Before one considers to solve such systems, one may also analyze them further as exemplified in [112] in order gain further insight or to find further relations among them.

Uncoupling algorithm and scalar solvers
In the last years a general toolbox has been elaborated that finds all solutions that can be given in terms of iterated integrals (or sums) as follows. By uncoupling algorithms [113,114] available, e.g., in the package OreSys [115], one first decouples the system (9) to a scalar linear differential equation in one of the unknowns. Using the differential equation solver in HarmonicSums [92] (based on [93,94,59]) one finds, whenever possible, a closed form representations of the unknown functions f 1 , . . . , f λ in terms of d'Alembertian (and partially of Liouvillian) solutions. Based on this strategy we recalculated the 2-loop form factors [116] and obtained first results for the 3-loop case [117]. Another fruitful approach [118] is based on recurrence solving. Here one assumes that the arising Feynman integrals can be given in the power series representations Then the machinery proceeds as summarized in Fig. 1. After the decoupling of the system (step 1) one takes the scalar differential equation of f 1 (x) and calculates by means of holonomic closure properties [88] a scalar linear recurrence of F 1 (n) (step 2). Activating the recurrence solver of Sigma (based on [55, 57-59, 75, 81]) one can decide algorithmically in step 3 if F 1 (n) can be represented in terms of indefinite nested sums. If yes, one plugs this representation into the decoupled system and gets closed forms of the remaining F 2 (n), . . . , F λ (n) in step (4)

Direct solver
So far there are only few algorithms available that can compute directly (i.e., without uncoupling) the desired set of solutions of a given coupled system of the form (9). For instance, with the algorithms in [124] one can find all hyperexponential solutions of (higher order) coupled systems but so far no algorithms are available to compute all d'Alembertian or Liouvillian solutions. First steps have been elaborated in [125] within the general difference field setting. However, in various instances it has been observed in [126] that the arising coupled systems of the form (9) can be transformed to a system of the form for a matrixÃ(x) which is free of ε and wheref andg are defined by the multiplication of an invertible matrix T with f = (f 1 , . . . , f λ ) and g = (g 1 , . . . , g λ ), respectively. Under the assumption that such a transformed system exists, algorithms are available to compute such a transformation matrix T . Furthermore, methods exist, cf. [127][128][129][130] that can hunt for such a basis transformation for the multivariate case, i.e., for systems of partial linear differential equations. Given such a transformed system one obtains the benefit that one can read off the coefficients of the ε-expansion in terms of indefinite nested integrals. However, such a transformation does not hold for more complicated systems. As elaborated in [104], there are other special cases that enable one to solve coupled systems of partial linear differential equations if one considers the solution as a multivariate power series solution where the coefficients satisfy a nicely coupled system of linear difference equations. For instance, take the coupled partial system x n x m the coupled system of homogeneous firstorder difference equations Given such a first-order homogeneous system, it follows that its solution can be expressed in terms of hypergeometric products. Namely, using the algorithm given in [104,Sec. 4.1] (which is a simplified version of the algorithm given in [131]) and implemented in the package HypSeries one obtains the solution in terms of hypergeometric products or equivalently in terms of factorial and Pochhammer symbols. As a consequence the derived solution of the original coupled system of differential equations can be given in the form Using this sum representation one can deploy the summation tools in Section 6 to calculate the first coefficients of its ε-expansion. We note that these solutions (coming from homogeneous first-order difference systems) are closely related to special functions that are introduced in the next section.

Transformation to special integral and sum representations
In the simplest cases, the integrands of Feynman integrals (1) exhibit Euler Betafunction structures and by clever rewriting, cf. also e.g. [132], the integral can be rewritten in terms of hypergeometric functions and their generalization [133][134][135][136]. More precisely, one may rewrite simple Feynman integrals in terms of the following hierarchy of p+1 F p functions, the first of which read Here the parameters a i , b i are such, that the corresponding integrals exists, [136]. We note that computer algebra can be used non-trivially to explore further properties on these special functions. E.g., using symbolic summation (see also Section 6) one can compute all arising contiguous relations of a finite set of sums [139]. More generally, it is possible to represent Feynman integrals by Mellin-Barnes [140][141][142] representations [81].
When one succeeds in detecting that the the given Feynman integrals can be rewritten in integral representations that can be connected to p+1 F p or Appell-like functions, one can utilize the essential property that all the p+1 F p functions have a single infinite sum representation, while the Appell-functions are represented by two infinite sums. For instance, we get Similarly, there are also other classes of higher transcendental functions, which obey multi-sum representations [146,147,149,104]. Up to the level of massless and single mass two-loop integrals, cf. [153] and in some cases in the three loop case, cf. [118], these representations are usually sufficient. For more complicated integrand-structures, however, one has to apply other techniques. Applying successively Newton's binomial theorem and Mellin-Barnes [140][141][142] decompositions on the integrand, implemented in different packages [154][155][156][157], enables one to carry out all integrals by introducing Mellin-Barnes integrals. Finally, carrying out the remaining Mellin-Barnes integrals with the residue theorem yields definite multiple sums Here the upper bounds L 1 (n), . . . , L v (n, k 1 , . . . , k v−1 ) are integer linear (i.e., linear combinations of the variables over the integers) in the dependent parameters or ∞, and f is hypergeometric in n and the summation variables k i . For further details on this rewriting process we refer, e.g., to [81,28]. In physical applications the dimensional parameter ε arises in the parameters a i , b i , c i , ... of the p+1 F p and Appell-type representations. Moreover, the hypergeometric summand f in (12) may also depend on ε. In all these cases one seeks for an εexpansion where the coefficients are represented in sum representations that are as simple as possible. In order to accomplish this task, highly general summation methods introduced in Section 6 can be applied.

Symbolic summation
Following the strategy sketched in Section 5 one ends up at (thousands or even millions) of definite multiple sums where the summand is built by hypergeometric products and indefinite nested sums, like harmonic sums [61,62], cyclotomic harmonic sums [63], generalized harmonic sums [64,65]; these sums may pop up in particular if one expands the summand w.r.t. the ε-parameter (i.e., if one applies the differential operator w.r.t. ε to the hypergeometric products; for a detailed description see, e.g., [104]). Producing such sum representations without making the original problem more complicated is highly non-trivial. However, if one succeeds in getting an appropriate sum representation, one can apply various symbolic summation algorithms to simplify these sums.

Simplification of indefinite nested sums
The simplification of indefinite nested sums defined over hypergeometric products started with Gosper's and Karr's summation algorithms [158,72] and has been enhanced significantly within the last 20 years to a strong summation machinery based on difference field and ring theories [73,[159][160][161][162][163]. Using our summation package Sigma [77,78] it is now possible to design completely automatically appropriate difference rings in which one can represent such indefinite nested sums fulfilling various optimality criteria: e.g., the number of nested summation quantifiers or the degrees in the denominators are minimized; see [164][165][166][167]. Furthermore, employing our contributions to a refined Galois theory of difference rings [168] (see also [69,[169][170][171]), the used sums do not admit any algebraic relations. As a consequence, one obtains canonical (unique) product-sum representations [79].
Furthermore, these algorithms can be accompanied with quasi-shuffle relations [172][173][174][175][176] for the discovery of such relations in a very efficient way; for further details we refer to Section 9.

The WZ-summation approach
The treatment of single nested definite hypergeometric sums started with Zeilberger's creative telescoping paradigm [177][178][179][180][181] and has been enhanced to multi-summation with the WZ-summation approach due to [182] and its refinements given, e.g., in [183][184][185]. Given a multiple sum F (n, ε) over a hypergeometric summand, like on the left-hand side of * * one can search for a recurrence/difference equation of order λ of the form λ i=0 a i (n, ε) F (n + i, ε) = r(n, ε) with polynomials a i (n, ε) in n, ε and r(n, ε) being an expression in terms of multiple sums of simpler type than F (n, ε). By further tricks one can compute even a homogeneous recurrence. With the package MultiSum [183] one obtains, for instance, for the triple * * For a ∈ Z \ {0} we define the generalized harmonic numbers S a (n) = n k=1 (sign(a)) k k |a| . sum in (13) a homogeneous linear recurrence with polynomial coefficients in n of order λ = 4 in about 2 days. Given this recurrence, one can utilize algorithms from [55,57,58,69,75,70] (see Section 3.1.1) encoded in our package Sigma that find all d'Alembertian solutions, i.e., all solutions that can be expressed in terms of indefinite nested sums defined over hypergeometric products. More precisely, Sigma computes 4 linearly independent solutions (i.e., their linear span generates all solutions) where is the most complicated sum solution. Finally, with four initial values of the triple sum one finds an alternative representation of it in terms of indefinite nested sums.
In general, these are highly nested, and the summands might consist of ugly polynomials in the denominator (like in (14)) that do not factorize nicely. However, employing our sophisticated difference ring algorithms introduced in Section 6.1, one can simplify the found representation further and obtains the right-hand side in (13). In total, the solving and simplification steps need around 10 seconds.
Summarizing, combining the WZ-approach (recurrence finding) and solving tools, one obtains a summation machinery that can transform a definite nested sum to expressions in terms of indefinite nested sums. When the input sum depends furthermore on the dimensional parameter ε, this machinery has been generalized in [81] to determine the coefficients of the ε-expansion of (2) whenever they are expressible in terms of indefinite nested sums defined over hypergeometric products. This toolbox is very general, but has a substantial drawback: it reaches already with such simple sums like in (13) its limit. With the difference ring approach described next, this situation can be improved substantially.

The difference ring approach
With the difference ring and field theories worked out in [72,159,165,167,161,168] one can simplify not only indefinite nested sums, but one can also apply Zeilberger's creative telescoping paradigm [177]. This means that one can try to compute a linear recurrence of order λ for a definite sum, say S(n) = n k=0 f (n, k), where f (n, k) is given in terms of indefinite nested sums defined over hypergeometric products w.r.t. the summation variable k. Given such a recurrence, one can solve it in terms of indefinite nested sums defined over hypergeometric products by the algorithms given in Section 3.1.1. If one succeeds in combining the solutions accordingly (matching λ initial values), one obtains an alternative representation of S(n). If this expression itself is summed over n, one can repeat this process w.r.t. another variable (over which one may sum later again). In a nutshell, one can apply the summation spiral illustrated in Fig. 2 iteratively with the goal to transform a given multi-sum from inside to outside to a representation purely in terms of indefinite nested sums.   This interplay has been automated in the package EvaluateMultiSums [78] based on Sigma's difference ring algorithms and produces the right-hand side in (13) in about 70 seconds. If a sum depends also on the dimensional parameter ε, one can first expand the summand of the multi-sum w.r.t. ε and can apply afterwards the summation quantifiers to each of the coefficients being free of ε. A clear drawback of this approach is that the summands blow up when higher ε-orders are calculated. Nevertheless, the pure difference ring approach produces simplifications that currently no other toolbox can achieve. E.g., while treating the 2-mass 3-loop integral given in Fig. 3 triple and quadruple sums between 0.4 to 1.6 GB of memory arose in [187] that could be simplified to expressions in terms of binomial sums using 8.4 MB of memory only. Further challenging calculations based on the difference ring/field approach can be found, e.g., in [118,[187][188][189].

The holonomic-difference ring approach
Another prominent branch of symbolic summation is the holonomic system approach which has been introduced in [190] and pushed further, e.g., in [191,192] to determine recurrence relations. Here the summand of a definite sum is described by a system of homogeneous recurrences with polynomial coefficients. Then given such a system, one can try to compute a linear recurrence system by introducing the next definite summation quantifier. Applying these algorithms iteratively from inside to outside yields a linear recurrence in n for the input sum. However, the underlying recurrence systems may grow heavily and the holonomic approach usually fails due to time and memory limitations. In [193] a hybrid strategy has been introduced and developed further in [194,195] that brings the holonomic and difference ring/field approach under a common umbrella. This new approach allows one to deal with recurrence systems with inhomogeneous parts in terms of indefinite nested sums that covers the pure holonomic and difference ring approaches as special cases. So far this approach has been nontrivially applied to obtain the first computer assisted proof [196] of Stembridge's TSPP theorem [197] and to provide the first proof of a non-trivial identity in [198] that is connected to irrationality proofs of zeta-values. In QCD-calculations this new approach has been explored further to evaluate, e.g., bubble topologies [199].

Symbolic integration
In the following we will present some of the most relevant tools of symbolic integration that have been used (at least in parts) in particular for multi-loop calculations in the case of a few number of external legs in elementary particle physics. Other tools suited for lower loop multi-leg calculations are described in part e.g. in [12].

The hyperlogarithm approach
If a Feynman diagram of the form (1) has no pole terms in (2) (i.e., l = 0) or can be made finite by certain transformations splitting off its pole terms [200], it can be calculated under certain conditions by using the method of hyperlogarithms [201]. Since here the denominator of the integral (1) is a multinomial in the Feynman parameters x i ∈ [0, 1], one may seek a sequence of integrations, such that the denominator is always a linear function in the integration variable. In this case the Feynman integral can be found as a linear combination of Kummer-Poincaré iterated integrals (also known as Goncharov iterated integrals), [202][203][204][205][206]. The method has been first devised for massless scalar integrals in [201], for a corresponding code see [207], and it has been generalized to massive diagrams [208], dealing even with cases with no thorough multi-linearity, which is an extension to [201,207].

The multivariate Almkvist-Zeilberger approach
Similar to the WZ summation approach its continuous version, the multivariate Almkvist-Zeilberger algorithm [185], can compute a linear recurrence/difference equation for a Feynman integral of the form (1). Likewise, if the Feynman integral depends on a continuous parameter x and the integrand is hyperexponential in x, one can search for a linear differential equation of the form ε).
A refined and improved method for the input class of Feynman integrals has been developed [209,118,96] which can hunt efficiently for homogeneous recurrences or differential equations. E.g., a recurrence in n of order 5 can be calculated in about 8 hours for the master integral dx dy du dv that arose in the context to tackle the highly non-trivial 3-loop Feynman diagram given in Fig. 4; see [118]. Then using the the linear difference equation solver of Sigma one can compute the first coefficients of the ε-expansion (2) in terms of harmonic sums and generalized harmonic sums. More generally, one can utilize the algorithms from [81,96] to solve linear difference and differential equations in terms iterated sums and integrals; for further details see Section 3.

The differential field and holonomic approach
Risch's algorithm [210] for indefinite integration (for details see [211]) allows as input an integrand from the class of elementary functions (they are recursively built by compositions of algebraic, exponential or logarithmic functions and the standard operations +, −, ×, /) and one can decide, if the indefinite integral defined over the input function can be written again in terms of elementary functions. Inspired by this result many further extensions have been derived. In particular, with [212,211] it is possible to deal with special classes of Liouvillian integrands (recursively built by indefinite integrals and hyperexponentials). In this regard, e.g., the package Integrator [213] enables one to treat not only indefinite integration problems, but also to compute difference/differential equations if the integrand depends on a discrete/continuous parameter. These tools have been exploited, e.g., to study root-valued integrals in [66] that arise within massive 3-loop Feynman integral calculations.
In particular, the holonomic system approach [190][191][192] can be applied not only to multi-sums, but also to multi-integrals of the form (1) to determine a linear recurrence in a discrete parameter n or a linear differential equation in a continuous parameter x. Analogously to the sum case, one can compute stepwise systems of linear differential/difference equations working from inside to outside of the multi-integral and ending up at a scalar equation of the free parameter n or x. First examples have been elaborated in [192,214] using the package HolonomicFunctions that illustrate further possibilities in QCD-calculations.

The method of arbitrarily large moments
One is often interested in the calculation of a certain number of moments in the Mellin variable, say n = 0, 1, 2, . . . , µ, to predict extra properties of physical quantities in terms of Feynman integrals. Standard procedures, like Mincer [215] or MATAD [216], allow the calculation of a comparable small number of Mellin moments, e.g., µ = 20. Recently, a new method has been worked out in Ref. [217] and implemented within the package SolveCoupledSystem [218,50,219] to compute thousands of such moments.
In general, this new method assumes that we are given a coupled system (9) with (10) where already µ moments for the inhomogeneous part in (9) are computed (by applying this method recursively). Then given such an input, one follows the calculation steps in Fig. 1 but instead of solving the recurrence in step (3), one uses the recurrence together with a small number of initial values of F 1 (n) (bounded by the order of the recurrence) to compute in linear time the moments F 1 (0), . . . , F 1 (µ), and finally the corresponding moments for F 2 (n), . . . , F λ (n). If the F i (n) depend also on ε, one can calculate the moments of the coefficients of the ε-expansions by exploiting refined ideas from [81].
More generally, using IBP methods [105,106], we suppose that a physical expression f (x) is given in terms of master integrals that are described in terms of recursively defined coupled systems of differential equations. Then using the large moment method iteratively one can calculate for a very large µ the moments of the master integrals. Assembling all the building blocks in the physical expressionf (x), one finally derives at the coefficients F (0), . . . , F (µ) of its power series (4) in terms of rational numbers (if ζ-values and other constants arise linearly, they are separated accordingly).
While in traditional solving methods very complicated function spaces might arise in intermediate steps, the large moment method deals simply with rational numbers and one can represent physical quantities with such sequences without entering any structural challenges.
Given a large number of moments, one may follow various strategies illustrated in Fig. 5. First, one can try to obtain interpolation expressions, e.g., by using orthogonal polynomials [220] that provide numerical data of sufficient high precision being relevant for the experiments at the LHC and other future colliders.
Second, one can apply the guessing methods from Section 2.2 in order to produce linear recurrences with minimal order for the physical quantities. In short, analyzing this quantity amounts precisely to the exploration of the computed recurrence. Next, one can try to solve the recurrences in terms of special functions by using the tools form Section 3. This strategy is particularly successful if the final result (but not necessarily the intermediate results) can be given in terms of indefinite nested sums over hypergeometric products. As demonstrated in [80], we could calculate from about µ = 5000 moments all the recurrences that determine the massless unpolarized 3-loop anomalous dimensions and Wilson coefficients in deep-inelastic scattering [39][40][41] by solving the recurrences. Similarly, we could calculate, e.g., the 3-loop splitting functions [44], the massive 2-and 3-loop form factor [49,50], the anomalous dimensions from off shell operator matrix elements [45,52,46], the polarized transition matrix element A gq (N ) [53] and others, the logarithmic contributions to the polarized O(α 3 s ) asymptotic massive Wilson coefficients [54], and the two-loop massless off-shell QCD operator matrix elements [47]. We remark that the found recurrences may also contribute substantially in the case that one fails to find closed form solutions. For instance, one may extract the asymptotic behavior of the physical quantities by using methods described in [221,222].

Special functions and their algorithms
The representation of the results of calculations in QCD and QED are characterized by special constants and functions. The former ones appear in zero scale calculations and as boundary conditions in single and more scale problems. Since in particular in QCD and QED the Mellin transform (see (3) where in the following x is replaced by z) relates nested sums at the one hand to nested integrals at the other hand, and vice versa, two principle classes of special single scale functions emerge: indefinitely nested sums over hypergeometric products and iterated integrals over certain alphabets of letters. Both in the limit n → ∞ of the sums and at z = 1 for the iterated integrals special numbers are obtained. Examples on different classes of functions are given in Tab. 1. All these function spaces obey (quasi) shuffle relations, cf. [172,173], implying algebraic relations, which allow to reduce to the respective algebraic bases [173,176].
A systematic description in terms of harmonic sums started in 1998 with [61,62]. They are defined by S b, a (n) = n k=1 (sign(b)) k k |b| S a (k), S ∅ = 1, b, a i ∈ Z\{0}, n ∈ N\{0}. (15) Related to that, the iterative integrals are the harmonic polylogarithms, Binomial Sums root-valued iterated integrals associated numbers with the alphabet [91]. A special defintion is required for the case H a (z), ∀a i = 0, i = 1...n, which has no integral representation, but is defined as ln n (z)/n!, for completeness. In the case of infinite sums we also allow for the symbol σ ∞ := ∞ k=1 (1/k), which is not a number, but simplifies various algebraic relations and is therefore useful.
The special numbers are multiple zeta values in both cases. Their representations at high order can be found in [35].
At the next level, generalized harmonic sums [64,65] contribute, e.g. in the case of the pure singlet 3-loop massive Wilson coefficients in the asymptotic region [43]. These quantities are given by The corresponding iterated integrals are also called Kummer-Poincaré iterated integrals [202][203][204][205][206] and are given by with the alphabet Further, cyclotomic harmonic sums and polylogarithms [63] contribute. The letters of the alphabet forming the iterated integrals are those of the harmonic polylogarithms extended with letters of the type with k labeling the cyclotomic polynomials and a ∈ [0, ϕ(k)], and ϕ(k) is Euler's totient function. The associated cyclotomic harmonic sums iterate monomials of the type Finite binomial sums [66] contribute for a series of topologies in the massive OMEs A Qg [118]. The corresponding sums are generalized sums with an additional factor of 2k k in the numerator or denominator. The associated iterated integrals, obtained by a Mellin inversion, are formed out of letters containing square root valued structures, as e.g. shown in Tab. 1. Another example is dz.
In particular, iterative application of integration by parts yield the asymptotic expansion Such expansions are extremely useful for limit calculations and for analyzing the expression behavior for large values of n. Moreover, the sum and integral representations equipped with their shuffle and quasi-shuffle algebras [172][173][174][175] give rise to algebraic relations of infinite sums. In particular, attaching special constants to sums that cannot be simplified further, one can discover evaluations such as (2i+1) 2 denotes the Catalan constant; for further details see [230,231]. For general classes, like nested binomial sums, more flexible methods were developed recently to map between n-and z-space, cf. [187]: given a recurrence of F (n) (resp. a differential equation of f (z)), compute a differential equation for f (z) (resp. a recurrence for F (n)). In particular, using the introduced solvers from Section 3.1, one can check, if the Mellin transform (resp. inverse Mellin transform) can be given in terms of indefinite nested sums (resp. integrals). Infinite (inverse) binomial sums have been also studied in [67,68]. For the simpler cases efficient rewrite rules have been developed to switch between the sum and integral representations via the (inverse) Mellin transform.
In more general cases, in particular in two-scale problems, the so called G-functions appear, which are iterated integrals over larger alphabets, partly with root-valued letters. They are given by Actually f c (x) even denotes in general a differentiable function, up to regularizations in special cases. At higher and higher orders in perturbation theory, new building blocks arise that cannot be represented in terms of indefinite nested sums or iterated integrals. In particular, one ends up at linear difference/differential equations, that cannot be solved completely in terms of d'Alembertian/Liouvillian solutions. For this reason, the class of iterative non-iterative integrals have been introduced in 2016 [232].
Probably the first case in which complete elliptic integrals emerged in a quantum field theoretic calculation has been the fourth order spectral functions for the electron propagator by Sabry 1962 [233]. For complete physical processes more recently elliptic integrals were needed. This is the case in massive three-loop calculations for the QCD corrections of the ρ parameter [98,99] 2017 † † and in massless three-three loop calculations 2018 and later [235,236]. There is a series of well-known examples of individual integrals of a certain structure in the literature as the sun-rise integral, cf. e.g. [237][238][239] and the kite-integral [233,240,241]; for a collection of recent surveys see Ref. [242]. ‡ ‡ In general, these classes of integrals form iterative non-iterative integrals, cf. [98]. Beyond this level one has Abel-integrals [246] and Calabi-Yau structures, cf. [247,24,248]. Even more involved structures will occur at higher topologies. In Mellin space they have the common characteristics of difference equations with rational coefficients which are not factorizing at first order. Any of the corresponding solutions also needs efficient numerical representations, as e.g. [249][250][251], for phenomenological and experimental applications. This also applies to Mellin space representations for n ∈ C, [252-254, 174, 175]. Feynman integrals will imply a multitude of new function spaces in the future.

Calculations in Quantum Field Theory
Our major topic concerns analytic Feynman diagram calculations. As has been shown, this is deeply rooted in solving large systems of differential or difference equations. The single scale cases are mathematically widely understood and one may project to the zero scale case, i.e. to special numbers. However, just by this one will probably not be able to find all relations between these special numbers by using, e.g., the techniques described in Section 2.1, beginning at a certain level of complexity, which requires further advanced methods in these cases. On the other side, as experience shows, certain twoscale problems can still be solved analytically, as we will discuss in Section 10.4. But already starting at that level, one has to deal with partial differential and difference † †The same differential equations rule the non first-order factorizing cases in the calculation of the massive three-loop operator matrix element A Qg in the single mass case [234], found together with the calculation [99].
‡ ‡Very naturally, as now the technical aspects on complete elliptic integrals are very well known in particle physics, many applications find these contributions, cf. e.g. [243][244][245], the reason being the occurrence of the corresponding Heun and 2 F 1 -type differential equations, cf. e.g. [98].
equations, on which is much less known, cf. [255,256]. In the single scale case, going to higher and higher orders, one will face non-first order factorizing differential and difference equations of higher and higher order [257,258,246] § §, for which only the properties of very few concrete classes have been studied so far and a very wide field of future mathematical investigation is opening up.
For more scales, one probably will have to rely on using numerical precision methods in the first place, because of the wide variety of structures [25]. Computational Quantum Field Theory (QFT) is urged to invest much more efforts to obtain fast and highly reliable methods in this direction to be able to cope with the challenges in future precision measurements. Developments of this kind may take quite a long time and need intense collaboration with experts in the field of numerical mathematics.

Zero Scale Calculations
Zero scale quantities in QFTs, as QED and QCD, are characterized by color factors, rational coefficients and special numbers like multiple zeta values [35]. Examples are fixed moments for massive three-loop OMEs [259] and massless four-loop anomalous dimensions [260]. Particularly for massive problems more special numbers contribute, as those related to generalized harmonic sums [65], cyclotomic harmonic sums [63], binomial sums [66], and those related to elliptic integrals [261], see also Tab. 1. More and more different sets will emerge including even higher topologies. One also may calculate moments of single scale quantities, which depend on an integer parameter n, by obtaining sequences of rational numbers. These numbers incorporate thus an essential part of the more involved single scale dependence for general values of n. It is sometimes of advantage to first work with these moments, despite the fact that the general n relation is determined by a difference equation, which does not factorize at first order. This is often the case for master integrals in the massive case. However, the corresponding recurrences for anomalous dimensions are factorizing at first order. One inserts first the master integrals for fixed moments and then determines the difference equation for the anomalous dimension, see [44,52].

Massless Single Scale Calculations
These quantities are the anomalous dimensions, currently known to three-loop order [39,40,262,43,44,52,45,46], the massless Wilson coefficients for deep-inelastic scattering [263] up to O(α 3 s ) [41], the Drell-Yan process and Higgs production to two-loop order [264][265][266][267]. All these quantities can be expressed by harmonic sums [62,61] in Mellin space and by harmonic polylogarithms in z-space [91]. For Higgs production and the Drell-Yan process at three-loop order [235,236] also elliptic integrals contribute. It is generally expected that a further nesting in the Feynman diagram topologies leads to new mathematical structures also in the massless case in higher orders of the coupling constant.
For massless single scale calculations one may very efficiently apply the method of arbitrarily high moments [217], together with guessing to obtain the recurrences, which may either be solved or reduced, by factoring of the first order factors, using Sigma. This also applies to the case of massive single scale calculations, to which we turn now.

Massive Single Scale Calculations
The method of massive OMEs [268] allowed to calculate single scale quantities, such as the heavy flavor Wilson coefficients to three-loop order in the asymptotic region, obtaining all logarithmic contributions [269,54] and also the constant term. The method has also been applied to problems in QED, cf. [270][271][272], see Section 10.3.1. In some cases even full results have been obtained at two-loop order [268,[273][274][275]270], cf. Section 10.3.2.  [268,276,277,153,278,279]. At three-loop order all but the massive OME A Qg have been calculated analytically in complete form both in the unpolarized and polarized case. In Mellin n space they can be expressed by harmonic sums for all N F -terms [280,122], and for A Qg receives also contributions by complete elliptic integrals [234].
Another case belonging to this class of integrals are the contributions to the massive form factor at three-loop order. The first order factorizing contributions can be expressed by harmonic and cyclotomic harmonic polylogarithms in the variable x q 2 with q 2 the virtuality and m the heavy quark mass, [49,283,117,284,285,50] The method of massive OMEs has also been applied for the calculation of the initial state radiation to the process e + e − → Z * /γ * . These corrections are of importance for planned future high luminosity measurements at the ILC, CLIC, and FCC ee. The results at O(α 2 ) [270] showed disagreement with an earlier direct calculation [286]. The only way to find the correct results has been a complete diagrammatic calculation, without expansion in the small parameter ρ = m 2 e /s, with m e the electron mass and s the cms energy. This has been performed in Ref. [287]. Furthermore, we expanded in ρ through different steps, controlled by high precision numerics, and confirmed the results of [270]. The fermionic integrals could be represented using iterated Kummerelliptic integrals over larger alphabets. Numerical results were presented in [288,289].
The method of massive OMEs has then been extended to calculate the first three logarithmic series up to O(α 6 L 5 ), where L = ln(s/m 2 e ) in [271]. Here in Mellin space also generalized harmonic sums contribute. Higher order corrections for the forward-backward asymmetry were calculated in [272], where also cyclotomic harmonic polylogarithms contribute to the radiators. ¶ ¶ 10.3.2. Massive Single Scale Calculations: including also power corrections In some cases, massive two-loop problems can be integrated analytically in the whole kinematic region. This applies to the flavor non-singlet contributions [268,273] and the pure singlet contributions [274,275]. While in the non-singlet case classical polylogarithms with root-valued arguments suffice for the representation, in the pure singlet case iterated integrals over alphabets containing Kummer-elliptic letters are necessary. Part of them integrates to incomplete elliptic integrals, which do not destroy the iterated integral structure, unlike the case in the iterative non-iterative integrals [291]. In establishing the contributing alphabet also rationalization of roots is performed as far as possible; for other investigations see also [292]. Examples for these letters are, cf. [274] and z is the momentum fraction variable. Depth-three iterated integrals over the contributing integrals emerge.
The fact that one finds analytic integral representations in these cases is related to the tree-like structure of the contributing diagrams. At higher orders or for other processes, correspondingly, one has to perform corresponding expansions in m 2 /Q 2 , to successively obtain analytic results, improving the logarithmic and constant orders obtained in the region Q 2 m 2 . The possibility to analytically calculate the puresinglet corrections, conjectured by van Neerven and J.B. around 2000, turned out to be correct, however, the necessary technologies for this became only available with [274] later and Nielsen integrals with whatsoever complicated argument are not sufficient to represent the final result.

Massive Double Scale Calculations
In the case of deep-inelastic scattering from the level of three loop onward, diagrams contribute, which contain charm and bottom quarks. This leads to a double scale problem, as similarly also in the case of the massive from factor and for other processes. In the following we will consider the case of deep-inelastic scattering. In all cases but the massive OMEs A Qg and its polarized counterpart ∆A (3) Qg , complete analytic solutions are possible either in Mellin n or momentum fraction z-space. The different OMEs ¶ ¶For a recent survey on the QED corrections see [290].
have been calculated in Refs. [293,294,187,295,189,188]. They can be expressed by G-functions, cf. (24). One example is [187] where η = m 2 1 /m 2 2 and χ = (1/η)( In the case of the pure singlet two-mass contributions [294,189] we work in z-space by using Mellin-Barnes integrals [140][141][142]. One also obtains G-functions and in part integrals over them, with a different support than usual, expressed by Heaviside functions. These problems cannot be solved in Mellin n space.

Classical Gravity
The classical kinematics of massive astrophysical objects, such as black holes and neutron stars, can be calculated by using methods of effective field theory developed for Quantum Field Theory. Concepts like the path integral [296] and Feynman diagrams are also applicable at the classical level. This is an enormous bonus to the field of general relativity and classical gravity, since very advanced computation technologies, starting from Feynman diagram generation [297], effective performance of Lorentz algebra [298,299], integration-by-parts reduction [107], and the calculation of master integrals already exist. One expands Einstein-Hilbert gravity in terms of auxiliary fields [300]. Furthermore, one performs the classical limit using the method of expansion by regions, cf. [301,302], where only the potential and radiation modes are contributing. These methods can be applied for the inspiraling process of the massive objects [303][304][305][306][307][308], as well as for their scattering process, cf. [309,310]. The main challenge is here to deal with the ever growing effective vertex structures and the integration by parts reduction, which can be performed by packages like Crusher, cf. [107]. The bound state kinematics is described by the post-Newtonian (PN) approach, having reached now 6PN order [311,308,306] and the scattering process by the post-Minkowskian approach, now available at O(G 4 N ), where G N denotes Newton's constant [309,310]. After expanding the potential contributions of the post-Minkowskain results, a part of the PN results is re-obtained, which has been shown to 6PN in [311,306]. This applies to the potential contributions. In principle, the method of guessing could be used obtaining post-Minkowskian results, again from potential contributions, [312]. While agreement has been reached for the 4PN level between various approaches, the level of 5PN is still under discussion because of the non-potential contributions. For their description the different methods proposed in the literature do not lead to consistent results as of yet, requiring both more clear theoretical foundations and also more work to obtain the final result for bound state problems.

Conclusion
For less than the last quarter century, a technological revolution has happened in the field of perturbative analytic calculations in renormalizable Quantum Field Theories, which is accompanying this field since. Considering single scale Feynman diagrams, in the time until 1998, the analytic integration of these amplitudes has been an art based on hypergeometric function structures and maximally Nielsen integrals dealing with sets of up to O(50) Feynman diagrams mostly to two-loop order. Before about this time computer algebra has been inspired and motivated often by its own discipline or by other mathematical research areas, such as combinatorics, number theory, or special functions. In this survey article we introduced recent methods from both communities and showed how they can be combined in non-trivial ways to new methods that may be instrumental for current and future precision calculations in particle physics. For a graphical summary of the different interactions presented in this article we refer to Fig. 6.
Elements of the revolution in Quantum Field Theories were also (quasi)shuffle algebras and the discovery of a hierarchy of function spaces both in Mellin n and momentum fraction z space. For all quantities considered one may find recurrences by applying the methods of arbitrary large moments and guessing, which delivers closed form equations in the first place. In the moment technology presented in Section 8 one simply deals with rational numbers ignoring completely the possible function spaces that may arise there. At the end of the day, one can apply the computer algebra tools introduced above and obtains from this data numerical representations or even symbolic representations of the final physical problem. In general, we feel that such new strategies will be crucial for future calculations and we are curious to see how these techniques can be developed further or can be outperformed with new ideas and strategies.
Systematic mathematical methods, like difference ring theory, allowed to reveal various new structures. Nowadays first order factorizing difference and differential equations (or systems thereof) are fully understood. Non-first order factorizing systems, indefinite summation y y special function algorithms q q Figure 6. The computer algebra and special function tools in interaction leading to 2 F 1 -solutions, complete elliptic integrals and modular forms are understood as well and steps in the direction of equations related to Calabi-Yau manifolds are done. Yet these are rather special systems only and Feynman diagrams can in principle cause more general, yet unknown structures also belonging to non-first order factorizing equations. They are fascinating as such and their complex analysis is a highly interesting topic. One may intend to derive general characteristics for these quantities [313]. Many of the present massless and massive three-loop problems of single and double scales could be solved by the technologies described in the present survey and new structures challenge further innovative mathematical solutions and efficient computeralgebraic implementations. All present achievements have in various instances been achieved by sophisticated computer algebra algorithms. Another challenge comes from the experimental possibilities at future colliders, operating at high luminosity, with which the theoretical results have to cope. All methods described do not only apply to relativistic renormalizable Quantum Field Theories, but also to effective field theories, e.g. dealing with (non-linear) Einstein general relativity in post-Newtonian and post-Minkowskian expansions at the classical level and various applications more, e.g. also to solid state physics. Problems with more scales do still escape complete analytic solutions at present and require more research in the future. The close collaboration of theoretical physicist, mathematicians and researchers in the field of computer algebra led both to the use of known algorithms from quite different fields in Quantum Field Theory, but have also triggered new mathematical and algorithmic research. The success reached has only been possible due to this symbiosis. This process will continue in full strength in the future.