About calculation of massless and massive Feynman integrals

We report some results of calculations of massless and massive Feynman integrals particularly focusing on difference equations for coefficients of for their series expansions


Introduction
Calculation of Feynman integrals gives basic information about both experimentally investigated processes and the characteristics of physical models. The calculations of the matrix elements of the processes under study depend on the masses of the particles involved in the interaction and, strictly speaking, require the calculation of Feynman integrals including those with massive propagators. Depending on the kinematics of the studied processes, the values of some masses can be neglected.
The study of the characteristics of physical models (for example, critical parameters, anomalous dimensions of particles and operators) usually requires calculations of massless Feynman integrals, which have a much simpler structure. They allow one to obtain results for these characteristics in high orders of perturbation theory.
Note that when calculating Feynman integrals, it is advisable to use analytical methods whenever possible. The fact is that approximate methods for calculating Feynman integrals do not always have sufficiently high accuracy. Moreover, the calculation of Feynman integrals is ambiguous for numerical calculations both due to the singular nature of the integrals themselves and (especially for gauge theories) with strong mutual contractions between the contributions of different diagrams or even between parts of the same diagram.
Moreover, in many important situations it is necessary to know exactly the exact results. For example, in the framework of renormalization group calculations in theories with high internal symmetry, it is important to know [1] about the vanishing of β-functions in the order under consideration.
We note that when using dimensional regularization [2], i.e. for an arbitrary dimension of space, once found diagrams for any field model (or process) can be applied to other models (and processes), since the main ones are scalar diagrams. Consequently, the complexity of analytical calculations of Feynman integrals pays off for their universality in application to various quantum-field models.
We also note the fact that the calculation of complicated diagrams may have some independent interest. So the use of non-trivial identities such as the uniqueness relation [3,4] can give information (see [5,6]) about the properties of some integrals and series that are not yet in the reference literature. For example, calculations of the same Feynman integral carried out in Refs. [5] and [7] using different calculation methods made it possible to find a previously unknown relation between hypergeometric functions with arguments 1 and −1. This relation has been proven only just recently [8].
tions for the coefficients of their decomposition are used. 1 Solving these recursive relationships allows you to get accurate results for these most complex parts. Obtaining these results is discussed in detail in Sections 3 and 6 for massless and massive diagrams, respectively.
A popular property of maximal transcendentality is shown in Section 7. It was introduced in [22] for the Balitsky-Fadin-Kuraev-Lipatov (BFKL) kernel [23,24] in the N = 4 Supersymmetric Yang-Mills (SYM) model [25], is also applicable for the anomalous dimension matrices of the twist-2 and twist-3 Wilson operators and for the coefficient functions of the "deep-inelastic scattering" in this model. The property gives a possibility to recover the results for the anomalous dimensions [26,27,28] and the coefficient functions [29] without any direct calculations by using the QCD corresponding values [30].
The very similar property appears also in the results of calculation of the large class of Feynman integrals, mostly for so-called master integrals [31]. The results for most of them can be reconstructed also without any direct calculations using a knowledge of several terms in their inverse-mass expansion [32]. Note that the properties of the results are related with the ones of the amplitudes, form-factors and correlation functions (see [33]- [41] and references therein) studied recently in the framework of the N = 4 SYM.
In Section 6, we demonstrate the existence of the propertiy of maximal transcendentality (or maximal complexity) in the results of two-loop two-and three-point Feynman integrals (see also [42,43]).

Basic formulas
Let us briefly consider the rules for calculation of massless diagrams. All calculations are carried out in momentum space.
The following formulas hold. 2 A. For simple chain: q µ 1 ...q µn q 2α 1 q ν 1 ...q νm q 2α 2 = q µ 1 ...q µn q ν 1 ...q νm q 2(α 1 +α 2 ) , or graphically i.e. the product of propagators is equivalent to a new propagator with an index equal to the sum of the indices of the original propagators. The number of momentums in the product in the numerator is equal to the sum of the products of the impulses in the original propagators.
As it was noted already, indeed, we use the transless product q µ 1 ...µn in the r.h.s. but really we need only the first term, i.e. q µ 1 ...q µn , because the rest (i.e. the terms containing g µ i µ j ) is exactly reconstracted from the exact form of traceless production.
The results have the following form Note that we use everything µ i belonging to the traceless product, i.e. we consider the case of a scalar diagrams with a traceless product. In real theories such as QCD, there are still other indices λ j , which correspond to the numerators of propagators. In such cases, we can no longer neglect the terms g λ i λ j and g µ i λ j and, therefore, the integration rules are complicated. They can be found in Refs. [15,16].
So, all diagrams, which can be expressed as combinations of loops and chains can be evaluated immediately. However, starting already with the two-loop level, there are diagrams, which cannot be expressed as combinations of loops and chains (simplest example is shown below in Fig.1). For these cases there are additional rules, which will be shown only graphically with a purpose to increase an illustation power.
C. When α i = d, there is so-called uniqueness ratio [3,4,5] for the triangle with indices α i (i = 1, 2, 3) where The results (10) can be exactlly obtained in the following way: perform an inversion q i → 1/q i (i = 1, 2, 3), k → 1/k in the subintegral expression and in the integral measure. The inversion keeps angles between momentums. After the inversion, one propagator is cancelled because α i = d and the l.h.s. becomes to be equal to a loop. Evaluating it using the rule (7) and returning after it to the initial momentums, we recover the rule (10). An extension of the rule (10) to the case with two traceless products can be found in [16].
D. For any triangle with indices α i (i = 1, 2, 3) there is the following relation, which is based on integration by parts (IBP) procedure [4,48] Eq. (12) can been obtained by introducing the factor (∂/∂k µ ) (k − q 1 ) µ to the subintegral expression of the triangle, shown below as [...], and using the integration by parts procedure as follows: The first term in the r.h.s. becomes to be zero because it can be represented as a surphase integral on the infinite surphase. Evalutiong the second term in the r.h.s. we reproduce Eq. (12).
As it is possible to see from Eqs. (12) and (13) the line with the index α 1 is distingulished. The contributions of the other lines are same. So, we will call below the line with the index α 1 as a "distingulished line". It is clear that a various choices of the distingulished line produce different tipes of the IBP relations.
Using equation (12) allows you to change the indices of the line diagrams by an integer. One can also change line indices using the point group of transformations [4,49]. The elements of the group are: a) the transition to impulse presentation, b) conformal inversion transformation p → p ′ = p/p 2 , c) a special series of transformations that makes it possible to make one of the vertices unique, and then apply relation (3) to it. An extension of the group of transformations for diagrams with the traceless product can be found in Ref. [16].

Basic massless two-loop integrals
The general topology of the two-loop two-point diagram, which cannot be expressed as a combination of loops and chanins is shown on Fig.1. Below in the present analysis we will concentrate mostly on two particular cases, which can be taken from the diagram shown in Fig. 1, for α 3 + n 3 = α, n 3 = n, α j (j = 3) = 1, n j (j = 3) = 0 (we denote I 1 (α, n)) and for α 5 + n 5 = α, n 5 = n, α j (j = 5) = 1, n j (j = 5) = 0 (we denote I 2 (α, n)) It is convenient to calculate the diagrams I 1 (α, n) and I 2 (α, n) using the functional relations 3 similar to those obtained in [6,50], which reduces the amount of computations.
Repeating analysis done in [6], we obtain the following functional relations: where the inhomogeneus terms are We see that the inhomogeneus terms in the functional equations (15) and (16), i.e. the diagrams I 11 (α, n) and I 21 (α, n), can be represented as combinations of loops and chains and, thus, they can be evaluated by rules (4) and (7).
We would like to note that for the massless two-point diagrams the subject of the study is so-called coefficient functions, which expression C i... (α, n) (i = 1, 2) for the considered diagrams I i... (α, n) can be introduced in the following form The result (19) contains the fact that we consider the two-loop diagrams. In general, L-loop diagram I L (α 1 , ..., α N , n) having propagators with indices α i (i = 1, ..., N) and one traceless product of momentums can be respresented as where α = N 1 α i . The results for C i... (α, n) can be obtained directly from rules (4) and (7). They have the following form: where and the result for A n,m (α 1 , α 2 ) is given in (6). Thus, the coefficient functions C i,1 (α, n) (i = 1, 2) can be represented as combinations of Γ-functions.
3.1 I 1 (0, n) and I 2 (0, n) The diagrams I 1 (0, n) and I 2 (0, n) can be considered a s natural boundary conditions for the functional equations (15) and (16). Moreover, in a sence, they have a special property: their results can be obtained with help of Eqs. (4) and (7) but with an additional resummatiomn.
Indeed, expanding for I 1 (0, n) and I 2 (0, n), respectively, the corresponding products of momentums in the following way: we see that the diagrams I 1 (0, n) and I 2 (0, n) can be represented as where So, the diagrams I 1 (0, n) and I 2 (0, n) can be expressed as combination of loops and chains and, thus, their coefficient functions can be calculated using rules (4) and (7).
So, we have for the coefficient functions C 1 (0, n) and C 2 (0, n): As we noted already at the beginning of the subsection, the results for C i (0, n) (i = 1, 2) are very important because they give a possibility to recover all results C i (m, n) using Eqs. (15) and (16) with α = m. But the calculation of C i (0, n) needs performing additional series (see Eqs. (29) and (30)). So, we present these calculations in some details, separetely for I 1 (0, n) and I 2 (0, n) in the following two subsections.

I 1 (0, n)
Now we consider C 1 (0, n) form Eq. (29), which can be represented as It is conveninet to rewrite the last term in the r.h.s. as that leads to the following form for C 1 (0, n): with where the normalization N 2 and the factors K 1 and K 2 are The result for C 1 (0, n) can be found exactly as where B(n + 1, a, b) = Γ(n + 1 + aε)Γ(1 + bε) Γ(1 + aε)Γ(n + 1 + bε) .
Evaluating the series Φ 1 (a, n) is quite difficut, and we will show this in a separate subsection. The calculation of more complicated series can be found in a famous paper [51].

Φ 1 (a, n)
It is convenient to consider the recursive relation between Φ 1 (a, n) and Φ 1 (a, n − 1). Indeed, we can represented Φ 1 (a, n) in the following form Taking partial fraction we can rewrite (47) in the following way The series in the r.h.s. can be calculated usng results forΦ 1 (α, β), studied in Appendix A, as It is convenient to introduce the new function Φ 1 (a, n) as which has the following relation: The last relation can be solved as Since Φ 1 (a, n = 0) = 0, then Φ 1 (a, 0) = 0 and we find finally and for a = 0 Note that the case Φ 1 (0, n) can be evaluated directly in the way similar to the one (37) for C 1 (0, n) and there is a full agreement of such calculations with (54).
After little algebra, we have where we used the formulas from the part C of Appendix A.
We would like to note about an appearence the nested sum S 2,1 (n), which cannot be obtained from expansions of products of Γ-functions.

C 1 (0, n)
Taking the results for φ
(114) below), we can replace in the r.h.s. the normalization and put the normalization N 1 in the definition of µ 2 g -scale of g-scheme [52], which relates with the usual MS one as µ 2 g = K 1 µ 2 M S (see discussions in Ref. [53]).

I 2 (0, n)
The diagram I 2 (0, n) s zero for odd n values. So, we can calculate it firstly at the even n values and recover its genetal form at the end. So, now we consider C 2 (0, n) from Eq. (30), which can be represented as where the normalization factor N 1 is given in Eq. (63). The last part of r.h.s. can be written as (for even n) Moreover, with the required accuracy O(ε 0 ), we can rewrite the product B(k + 1, −1, −2)B(n − k + 1, −1, −2) in the following way: because Then, C 2 (0, n) can be splitted to the four different parts: The parts C 2 (0, n) and C 2 (0, n) can be evaluated directly. Indeed, As in the case of C 1 (0, n), using Eq. (32) it is convenient to split the part C 2 (0, n) in two parts C where As in the case of C 1 (0, n), the part C 2 (0, n) can be found exactly: Since then C 2 (0, n) has the following form for even n values The result forC (2) 2 (0, n) can be expressed as The last part C It is easy to show that the part C After replacement k → n + 1 − k the last series obtains the additional factor (−1) n+1 and, thus, it is zero for even n values. It is really the case: the results for C 2 (0, n) are exactly evaluated in the part B of Appendix A. Now we return to the coefficient function C 2 (0, n) given in Eq. (69). It is convenient to split the result (64) into two parts: where The part C 2 (0, n) can be evaluated exactly as Using expansions (41) of B(n + 1, a, b) with respect of ε, we have The result forC 2 (0, n) can be expressed as To evaluate the r.h.s it is convenient to use the result (45) for φ 1 (a, n), to calculate the series Φ 2 (a, n): and to differentiate several times these series with respect of a and put a = 0. Indeed, where Φ (m) We consider the evaluation of the series Φ 2 (a, n) in the next subsection.
As in the case of Φ 1 (a, n), to obtain the results for Φ 2 (a, n) it is convenient to consider the difference relation between Φ 2 (a, n) and Φ 2 (a, n − 1). Indeed, we can rewrite (93) in the following way Taking partial fraction (48) we can represent (93) in the following way The series in the r.h.s. can be calculated usng results forΦ 1 (α, β) studied in Appendix A as It is convenient to introduce the new function Φ 2 (a, n) as which has the following reccursion: It can be solved as Since Φ 2 (a, n = 0) = 0, then Φ 2 (a, 0) = 0 and we find finally (101) We can calculate Φ 2 (0, n) directly in the way similar to the evaluation of Φ 1 (0, n). The direct calculation is a full agreement of such calculations with (102).
We have for Φ After little algebra, we obtain the final resuls We would like to note about an appearence the sum S −2 (n), which cannot be obtained from expansions of products of Γ-functions.
Taking the results for φ 1 (n) and Φ 2 (n) together, we can obtain the following results forC 2 (0, n) Taking the results for C 2 (0, n) andC 2 (0, n) together, we have for C 2 (0, n) Changing the nested sums S i (n) to S i (n + 1), we can simplify these results: We note, that we can put the normalization in the definition of µ 2 g -scale of g-scheme [52], which relates with the usual MS one as µ 2 g = K 1 µ 2 M S . Using the fact that the above result (107) is obtained for even n values (and C 2 (0, 2m + 1) = 0), we can recover the full resul for C 2 (0, n) in the following form: 3.4 I 1 (1, n) and I 2 (1, n) Now we evaluate the particular cases I 1 (1, n) and I 2 (1, n), which are important for the future studies. Indeed, we will see that the integral are finite and have very compact form.
The results for the C 1,1 (1, n) and C 2,1 (1, n) are where the normalization factors N 1 and N 2 and also the factors K 1 and K 2 are given in Eqs. (63) and (36), respectively. Below we consider the results for I 1 (1, n) and I 2 (1, n) is the different subsections.
Using the fact that the above result (121) is obtained for even n values (and C 2 (1, 2m + 1) = 0) , we can recover the full resul for C 2 (1, n) in the following form 4 Examples of calculation of some type of massless diagrams In the previous section, we studied the basic diagrams I 1 (α, n) and I 1 (α, n) with α = 0 and α = 1. I would especially like to note the case α = 1 where the coefficient functions of the integrals I 1 (1, n) and I 1 (1, n) have a very compact form.
In this section, we consider the expansion coefficients of scalar diagrams (we call them the "moments" of diagrams), arising in the study of forward elastic scattering. These moments are extracted from the initial diagrams by the method of "projectors" [13,14]. Some basic diagrams and getting their points are discussed in Appendix B.
Here we look at the diagrams I 1 (α, n) and I 1 (α, n) with α = n + 1 and α = n + 2. 4.1 I 1 (n + 1, n) and I 2 (n + 1, n) Let us first consider the two simplest diagrams: J 1 (α = 1, q, p) and J 2 (α = 1, q, p), shown in Fig. 2. As it was discussed already, using the method of "projectors" (see Appendix B) we can study the so-called moments J i (α = 1, n) (i = 1, 2). The moments of the diagrams shown in Fig. 2 are represented by the diagrams shown in Eq. (180) for α = n + 1, i.e. J i (α = 1, n) = I i (n + 1, n) To evaluate the diagrams it is convenient to use the so-called momentum transformation. To show it's effectivity it is useful to consider more complicated diagrams I i (α + n, n) (i = 1, 2), which have the coefficient functions C i (α + n, n) which is similar to (19). Now we introduce the following Fourier transforms [9] d d p e ipx 1 p 2α = 2 2α π d/2 a 0 (α) d d p e ipx p µ 1 ...p µn p 2(α+n) = (−i) n 2 2α π d/2 a n (α + n) where we neglect the terms of the order g µ i µ j . The Fourier transform (126) is the usual one (see, for example, the recent review [9]) and the Fourier transform (126) can be obtained from (126) with the help of projector Applijing the Fourier transforms to l.h.s. we will come to the following diagrams in x-space: x .
If we replace x and all internal coordinats by the corresponding moment q and the internal moments we obtain the corresponding diagram in the momentum space. We will call them as I i (α, n). Such procedure is called by "dual transformation" (see discussions in Ref. [9]) and denoted as d = (see aslo Introduction). So, we have Lets to denote the coefficient functions of the last integral as C i (α, n), i.e.
Consideration of the more complicated examples can be found in Refs. [15,16]. We only note here that for diagrams containing several propagators depending on the momentum p (see, for example,Ĵ 1 (α, β, q, p) in Appendix B), their moments will contain the sum of two-point diagrams (see the momentĴ 1 (α, β, n) in Appendix B).
Calculation of moments of this type is a much more serious problem compared to I i (n + 1, n) and I i (n + 2, n) with (i = 1, 2). Nevertheless, as was shown in Ref. [16], it is almost always possible to separate the contributions of complex integrals and complex sums. Moreover, ε-singularities remain only in the simplest parts, where it is usually possible to sum in all orders in ε.

Calculation of massive Feynman integrals
Feynman integrals with massive propagators are significantly more complex objects compared to the massless case. The basic rules for calculating such diagrams are discussed in Section 2, which are supplemented by new ones containing directly massive propagators. These additional rules are discussed in the next subsection.

Basic formulas
Let us briefly consider the rules for calculating diagrams with the massive propagators.
Propagator with mass M will be represented as The following formulas hold.
A. For simple chain of two massive propagators with the same mass: i.e. the product of propagators with the same mass M is equivalent to a new propagator with the mass M and an index equal to the sum of the indices of the original propagators.

B. Massive tadpole is integrated as
C. A simple loop of two massive propagators with masses M 1 and M 1 can be represented as hypergeometric function, which can be calculated in a general form, for example, by Feynman-parameter method. It is very conveninet, using this approach to represent the loop as an integral of a propagator with the "effective mass" µ [32,54,55,56]: It is useful to rewrite the equation graphically as D. For any triangle with indices α i (i = 1, 2, 3) and masses M i there is the following relation, which is based on integration by parts (IBP) procedure [48,4,54] As it was in the massless case (see Eq. (12)), Eq. (146) can been obtained by introducing the factor (∂/∂k µ ) (k − q 1 ) µ to the subintegral expression of the triangle and using the integration by parts procedure as in Eq. (13).
As it was in the massless case, the line with the index α 1 is distingulished. The contributions of the other lines are same. So, we will call below the line with the index α 1 as a "distingulished line". It is clear that a various choices of the distingulished line produce different tipes of the IBP relations.

Basic massive two-loop integrals
The general topology of the two-loop two-point diagram, which cannot be expressed as a combination of loops and chanins is shown on Fig.4. Below in the present analysis we will concentrate mostly on two-loop two-point and threepoint diagrams, which can be taken from the diagram shown in Fig. 4. We will call them as: Application of the IBP procedure [48] to loop internal momenta leads to relations between various Feynman integrals and, therefore, to the necessity of calculating only some of them, which in a sense are independent. These independent diagrams (which were chosen completely arbitrarily, of course) are called master integrals [31].
Applying the IBP procedure [48] to the master-integrals themselves leads to differential equations [18,54] for them with the inhomogeneous terms containing less complex diagrams. 5 Applying the IBP procedure to diagrams in inhomogeneous terms leads to new differential equations for them with new inhomogeneous terms containing even more less complex diagrams (≡ less 2 complex ones). By repeating the procedure several times, in the last step we can obtain inhomogeneous terms containing mainly tadpoles, which can be easily calculated in-turn (see also the discussions in Section 7 below).
Solving the corresponding differential equations in this last step, diagrams for the inhomogeneous terms of the differential equations in the previous step can be reproduced. Repeating the procedure several times, me can get the results for the original Feynman diagram.

Evaluation of series
As the first example, we consider P 126 integral, whous graphical representation has the following form This integral was calculated in Ref. [58] and presented also in [32,59]. It has the following series representation: 6 where Here the normalizationN = (µ 2 /m 2 ) 2ε , where µ = 4πe −γ E µ is in the standard MS-scheme and γ E is the Euler constant.
Eq. (149) contains the function F (n): 6 In Minkowski space all ln n x → ln n (−x) and there is an imaginary part which is important for cuts of the diagram. Such property exists also in other examples (see subsection 6.3 below).
Our purpose here to show the calculatio of the series F i (n). The most important point is the following: if some series has the following form Technically, the series F i (n) can be evaluated by considering a connection between F i (n+ 1) and F i (n), which can be expressed as a difference equation of the following type: with some coefficient functionsC(n). Usually the new functionF (1) (n) can be summed by standard formulas from techbooks. If it is not the case, then it is necessary to repear the above procedure for the functionF (1) (n), then we will come to new function, e.g.,F (2) (n). So, repeating the above procedure for the functionsF (i) (n), we will come to the oneF (i+1) (n), which can be summed by standard formulas from techbooks.
For all the above functions F i (n) (i = 1, 2, 3, 4) we should repeat the procedure twicely. Is it not so convenient. More simple to consider the new functions F i (n, a): which can be used to reproduce F i (n) as The new functions F i (n, a) can be evaluated with the usage of the above procedure only one time. We will havẽ where the new functions F i (1) (n) can be summed by standard formulas from techbooks.
Here we will present the exact evaluation only the function F 1 (n, a). The other functions can be evaluated similiry and we will show only the basic steps.

Properties of series
This scheme was successfully used to calculate the two-loop two-point [18,54,57] and threepoint diagrams [58,32,55] with one nonzero mass. This procedure is very powerful, but rather complicated. However, there are some simplifications based on representations of series of Feynman integrals. Indeed, the inverse-mass expansion of two-loop two-point and three-point diagrams 7 with one nonzero mass (massless and massive propagators are shown by dashed and solid lines, respectively), can be considered as The diagrams are complicated two-loop Feynman integrals that do not have cuts of three massive particles. thus, their results should be expressed as combinations of Polylogarithms. Note that we consider only three-point diagrams with independent upward momenta q 1 and q 2 , which satisfy the conditions q 2 1 = q 2 2 = 0 and (q 1 + q 2 ) 2 ≡ q 2 = 0, where q is downward momentum. where x = q 2 /m 2 (as it was shown in (149), η = 1 or −1 and α = 1 and 2 for two-point and three-point cases, respectively. The normalizationN is defined below Eq. (149). Moreover, for diagrams with two-massive-particle-cuts (2m-cuts). For the diagrams with one-massiveparticle-cuts (m-cuts) C n = 1.
For m-cut case, the coefficients F N,k (n) should have the form where S ±a,... ≡ S ±a,... (j−1) are harmonic sums defined in (42) and ζ(±a) are the Euler-Zagier constants For 2m-cut case, the coefficients F N,k (n) can be more complicated The terms ∼ V a,... and ∼ W a,... can come only in the 2m-cut case. The origin of the appearance of these terms is the product of series (174) with the different coefficients C n = 1 and C n =Ĉ n .

Other examples
As an example, consider two-loop two-point diagrams I 1 and I 12 studied in [32].
Their results are From (181) one can see that the corresponding functions F N,k (n) have the form if we introduce the following complexity of the sums (Φ = (S, V, W )) The number 3 − N determines the level of transcendentality (or complexity, or weight) of the coefficients F N,k (n). The property greatly reduces the number of the possible elements in F N,k (n). The level of transcendentality decreases if we consider the singular parts of diagrams and/or coefficients in front of ζ-functions and of logarithm powers. Thus, finding the parts we can predict the rest using the ansatz based on the results already obtained, but containing elements with a higher level of transcendentality.
Other two-loop two-point integrals in [32] have similar form. They were exactly calculated by differential equation method [18,54]. Now we consider two-loop three-point diagrams, P 5 and P 12 : Their results are (see [32]): Now the coefficients F N,k (n) have the form The diagram P 5 (and also P 1 , P 3 and P 6 in [32]) was calculated exactly by differential equation method [18,54]. To find the results for P 12 (and also all others in [32]) we have used the knowledge of the several n terms in the inverse-mass expansion (174) (usually less than n = 100) and the following arguments: • If a two-loop two-point diagram with a "similar topology" (for example, I 12 for P 12 , etc.) was already calculated, we should consider a similar set of basic elements for corresponding F N,k (n) of two-loop three-point diagrams but with a higher level of complexity.
• Let the diagram under consideration contain singularities and/or powers of logarithms.
Since the coefficients are very simple before the leading singularity, or the largest degree of the logarithm, or the largest ζ-function, they can often be predicted directly from the first few terms of the expansion.
Moreover, often we can calculate the singular part using a different technique (see [32] for extraction of ∼ W 1 (n) part). Then we should expand the singular parts, find the main elements and try to use them (with the corresponding increase in the level of complexity) in order to predict the regular part of the diagram. If we need to find ε-suppressed terms, we should increase the level of complexity of the corresponding basic elements.
Later, using the ansatz for F N,k (n) and several terms (usually less than 100) in the above expression, which can be exactly calculated, we obtain a system of algebraic equations for the parameters of the ansatz. Solving the system, we can obtain the analytical results for Feynman integrals without exact calculations. To check the results, we only need to calculate a few more terms in the above inverse-mass expansion (174) and compare them with the predictions of our anzatz with the fixed coefficients indicated above.
Thus, the considered arguments give a possibility to find results for many complicated two-loop three-point diagrams without direct calculations. Several process options have been successfully used to calculate Feynman diagrams for many processes (see [58,32,55,60]).
Note that properties similar to (183) and (186) were recently observed [35] in the so-called double operator-product-expansion limit of some four-point diagrams.

Modern technique of massive diagrams
Coefficients have the structure (183) and (186) with the rule (184). Note that these conditions greatly reduce the number of possible harmonic sums. In turn, the restriction is associated with a specific form of differential equations for the Feynman integrals under consideration. Differentials equations can be formally represented as [42,43] (x + a) d dx − kε FI = less complicated diagrams(≡ FI 1 ), with some number a and some function k(x). This form is generated by IBP procedure for diagrams including an inner n-leg one-loop subgraph, which in turn contains the product k µ 1 ...k µm of its internal momenta k with m = n − 3. Indeed, for ordinary degrees α i = 1 + a i ε with arbitrary a i of subgraph propagators, the IBP relation (12) gives the coefficient d − 2α 1 − p i=2 α i + m ∼ ε for m = n − 3. Important examples of applying the rule are the diagrams I 1 , I 12 and P 5 , P 12 , P 126 (for the case n = 2 and n = 3) and also the diagrams in in Ref. [61] (for the case n = 3 and n = 4). However, we note that the results for the non-planar diagrams (see Fig. 3 of [32]) obey the Eq. (186) but their subgraphs do not comply with the above rule. The disagreements may be related to the on-shall vertex of the subgraph, but this requires additional research.
Taking the set of less complicated Feynman integrals FI 1 as diagrams having internal n-leg subgraphs, we get their result stucture similar to the one given above (186), but with a lower level of complexity.
So, the integrals FI 1 should obey to the following equation Thus, we will have the set of equations for all Feynman integrals FI n as with the last integral FI n+1 contains only tadpoles. Moreover, following [62], we can reconstruct the above set of inhomogeneous equations as the the homogeneous matrix equation 8 (see Ref. [64] containg methods to obtain the equation) for the vector where the matrix K contains the functions k j /(x + a j ) as its elements. The form (190) is called as the "canonic basic". It is now very popular (see, for example, recent papers in Ref. [65]). Please note that for real calculations of FI n it is convenient to replace where the term FI n obeys the corresponding homogeneous equation The replacement simplifies the above equation (189) to the following form having the solution Usually there are some cancellations in the ratio FI n+1 /FI n and sometimes it is equal to 1. In the last case, the equation (193) coincides wuth the definition of Goncharov Polylogariths (see [66] and the references therein).
where the quantity a+b (a = i=1 |a i |) is fixed for each order of expansion in ε. For integrals I 1 (1, n) and I 2 (1, n) (and also I 1 (n + 1, n) and I 2 (n + 1, n)), this property is even more strict, since in this case b = 0. Unfortunately, calculations in QCD (or another theory) mix orders of magnitude over ε. Moreover, in expansions of diagrams containing propagators with non-trivial numerators (for example, as propagators of quarks and gluons) this property is also violated.
It is an amazing fact that the property (194) in its most strict form b = 0 is restored for diagonal elements of anomalous dimensions [26,27,28] and DIS coefficient functions [29] within N = 4 SYM.
The anomalous dimensions of the twist-2 Wilson operators govern the Bjorken scaling violation for parton distributions (≡ matrix elemens of the twist-2 Wilson operators) in the framework of Quantum Chromodynamics.
Balitsky-Fadin-Kuraev-Lipatov [23] and Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [70] equations resum, respectively, the most important contributions ∼ α s ln(1/x B ) and ∼ α s ln(Q 2 /Λ 2 ) in different kinematical regions of the Bjorken variable x B and the "mass" Q 2 of the virtual photon in the lepton-hadron DIS and, thus, they are the cornerstone in analyses of the experimental data from lepton-nucleon and nucleon-nucleon scattering. In the supersymmetric generalization of QCD the equations are simplified drastically [22]. In the N = 4 SYM the eigenvalues of the matrix of anomalous dimesion contain only one universal function with shifted arguments [71,26].
The three-loop result 9 for the universal anomalous dimension γ uni (j) for N = 4 SYM is [28] γ uni (n) =âγ The expression (199) is the analytical continuation (to real and complex n) [72] of the harmonic sums S −a,b,c,··· (n) (see discussions andout the analytical continuation and its applications in Appendix C) The results for γ (3) uni (n) [73,74], γ uni (n) [75], γ uni (n) [76] and γ (6) uni (n) [77] are obtained from the long-range asymptotic Bethe equations [78] for twist-two operators and the additional contribution of the wrapping corrections. The similar calculations for the anomalous dimension in the twist-three case can be found in [79].
Then, the basic functions γ uni (n) and γ (2) uni (n) are assumed to be of the types ∼ 1/n i with the levels i = 1, i = 3 and i = 5, respectively. A violation of this property may be obtained from contributions of the terms appearing at a given order from previous orders of the perturbation theory. Such contributions could be generated and/or removed using an appropriate finite renormalization and/or redefinition of the coupling constant. But these terms do not appear in the DR-scheme [80].
It is known, that at the first three orders of perturbation theory (with the SUSY relation for the QCD color factors C F = C A = N c ) the most complicated contributions (with i = 1, 3 and 5, respectively) are the same as in QCD [30]. This property allows one to find the universal anomalous dimensions γ   uni (n) without knowing all elements of the anomalous dimension matrix [26], which was verified for γ (1) uni (j) by the exact calculations in [27].
Note that in N = 4 SYM some partial cases of anomalous dimension are also known for the large couplings from string calculations and AdS/QFT correspondence [81]. We would like to note that if the property of the maximal transcendentality exists at low coupling, then sometimes it appeares at large couplings (see, for example, the results for the cusp anomalous dimension at low [82] and large [83] couplings, both of which are based on the Beisert-Eden-Staudacher equation [84]). However, this is not true for Pomeron intercept, which results lose the property of maximal transcendentality at large couplings (see [28,85,86]). The reason of the difference in the results for the cusp anomalous dimension and Pomeron intercept is currently not clear. More research is required.

Conclusion
In this review, we presented the results of the calculation of some massless and massive two-loop Feynman integrals. In the massless case, we considered scalar two-point diagrams with a traceless product in the numerator of one propagator and diagrams depending on two momenta q and p when p 2 = 0. In the massive case, we studied the 1/m 2 -expansion of Feynman integrals. The similarity of the structure of the expansion coefficients of massless and massive diagrams is shown.
The main calculation method for Φ i (n) and F j (n) is to build recurrence relations (i.e. the relationships between Φ i (n) and Φ i (n−1) and also F j (n) and F j (n−1)), which ingomogeneous parts contains only simpler amounts. The calculation of the series in the ingomogeneous parts allows us to get complete information about the recurrence relations and, as a result, solve them and get the desired values for the studied Φ i (n) (i = 1, 2) and F j (n) (j = 1, 2, 3, 4).
Note that this approach is close to the method of differential equations for calculating massive diagrams. This is not surprising, since the differential equations for any function correspond to recurrence relations for the coefficients of its expansion.
For the massless and massive diagrams under consideration, the level of transcendentality (or complexity, or weight) remains unchanged in any order of ε (see Subsections 6.3 and 6.4). Moreover, it decreases in the presence of logarithms or ζ-functions. It is called as transcendentality principle. Its application leads to the possibility to get the results for most of integrals without direct calculations.
This property is violated in physical models, such as QCD, where propagators (both quarks and gluons) contain momenta in their numerators, which are responsible for mixing levels of complexity. However, this property is preserved (see Section 7) for diagonalized quantities, such as diagonal anomalous dimensions and coefficient functions in N = 4 SYM, which is an amazing but poorly studied property.
Author thanks the Organizing Committee of Helmholtz International Summer School "Quantum Field Theory at the Limits: from Strong Fields to Heavy Quarks" for invitation. He is grateful to Andrei Pikelner for his help with the Axodraw2. .
So, we see that C and expess them trough standard sums S i (n) and the one S 2,1 (n). For the first series S 1,1 (n) we have Comparing l.h.s. and r.h.s., we see By analogy with the case S 1,1 (n) we can obtain Taking the results (A17), we have Comparing l.h.s. and r.h.s., we see The evaluation of the more complicated sums can be found in Ref. [87].
Next, we will neglect the symmetrizatorŜ. Note that this transformation from the diagram to its moment remains valid for arbitrary indices of the diagram lines, as well as the presence of additional momentums in the propagators of the diagram (if the latter are located on a differentiable line, then some changes will be required).
The second consider diagram is By analogy with the previous diagram we have for its moments: A similar conclusion can be drawn for the following diagram Its moment has the form J 1 (α, β, n) q ν 1 ...q νn q 2(n+α+2ε) = n k=0 Γ(k + β)Γ(n − k + α) (n − k)!k!Γ(α)Γ(β) Note that there is another technique to calculate the considered diagrams: the gluing method [12]. Using the orthogonality of traceless product, on can obtain the moment of the diagram by an additional integration on the q momentum with a propagator, which has an index δ and the additional tranceless product. This additional integration leads to very complicated three-loop diagrams. For the considered Feynman integrals J 1 (α, q, p), J 1 (α, q, p),Ĵ 1 (α, β, q, p) these gluing three-loop diagrams have the following form The evaluation of these complicated diagrams is above of the slope of the paper. Some example of application of the gluing method can be found in Ref. [15]. As a conclusion of Appendix B, we would like to note, when applying the method of "projectors" [13,14], the expression obtained for the nth moment of the initial diagram has alwais much simpler form than using the gluing method [12].

Appendix C, about analytic continuation
Here we demonstrate the direct calculation the particular case C 2 (1, n = 0) from the general expression C 2 (1, n) in (122) by using an analytic continuation (from even n values) of the sum Indeed, using the simple sum, it is very convenient to show the basic steps of the analytic continuation itself. In the case of more general nested sums S ±a,±b,... (n) formulas are more complex, which may make it difficult to understand the procedure.
The basic idea of the analytic continuation is very simple: try to translate the argument n from the upper limit of the sum to the summed expression. After the procedure performed we have a possibility to expand and to differenciate with respect of n and so on.
Firstly we represent the sum S −2 (n) in (C1) as and we see the unpleasant factor (−1) n in the front of the last term in r.h.s. Considering the variable (−1) n S −2 (n): we move the factor (−1) n to the front of the first term. Now we introduce the new n-dependent function S −2 (n) as S −2 (n) = (−1) n S −2 (n) + (1 − (−1) n )S −2 (∞) , which coincides with the initial one S −2 (n) at even n values and have no the unpleasant factor (−1) n So, the function S −2 (n) can be considered as the analytic continuation (from even n values) of the sum S −2 (n). Now it is possible to consider the small-n limit of C 2 (1, n) using the corresponding limit of S −2 (n) as Thus, in the small-n limit we have for S −2 (n): and for the coefficient function C 2 (1, n = 0) of (122) that exactly coincides with C 1 (1, n = 0). This analytic continuation has many important uses. For example, to analyze the evolution of parton distributions and DIS structure functions, there is a popular approach [89] which is based on the Gegenbauer polynomials, which in-turn are associated with the moments of parton distributions and structure functions. Using in the analysis a simple evolution for moments, which is determined here by the simple differential DGLAP equations, at the last step the parton distributions and/or structure functions are restored by summing (to a certain value N MAX ) the Gegenbauer polynomials.
In this analysis, evolution should be performed for both even and odd moments, so the analytic continuation is necessary. Using this analytical continuation, a lot of analysis of experimental data was performed (see a review in [90]) by the method described here.
Another important application [91] is the study of parton distributions and structure functions in the region of small values of the Bjerken variable x, which directly relates with above study of the nested sums at n → 0. The approch includes, in particular, an extraction of gluon density and longitudinal structure function from data for structural function F 2 , the evolution of parton distributions at low x in nucleon and in nuclei, an ultrahigh-energy asymptotics of the neutrino-hadron interaction cross section. Some review of these studies can be found in Ref. [92].