format: text

text:
# High-precision calculation of multi-loop Feynman integrals by difference equations

## 1 Introduction

Nowadays the technique probably more popularly used for calculating the contribution of Feynman diagrams is that based on the integration-by-parts in $`D`$ dimensions . After some algebraic operations, like contracting Lorentz indices and calculating fermion traces, the contribution of the diagram is expressed as a combination of several Feynman integrals with different powers of numerators and denominators.

This expression composed of many integrals is reduced to a combination of a limited number of ‘master integrals’ using recurrence relations obtained by combining the identities obtained by integration-by-parts. Then, the master integrals are calculated numerically or analytically, with some other method.

As this technique is applied to more and more complicated diagrams, some difficulties appear.

First, a working general algorithm for identifying master integrals and for obtaining such recurrence relations is not known at present111 An algorithm which in principle may solve this problem has been recently proposed in , but up to now no practical application was shown.. Up to now, for each diagram laborious handwork was needed for obtaining such recurrence relations222See section 2.3..

Second, the number of master integrals grows rapidly with the number of loops and legs of the diagrams. Considering for example the reduction to master integrals of the contribution to the $`g`$-$`2`$ of the electron in QED, it is known that the one-, two-, and three-loop contributions are reduced by integration-by-parts to respectively, 1, 3 and 17 master integrals (see  for the analytical calculation of the three-loop contribution). At the four-loop level, to which integration-by-parts has been still not applied, one expects several hundreds of master integrals.

Third, the calculation of a single multi-loop master integral is a difficult problem for which a general method applicable to any diagram, with any values of masses and momenta, able to provide high-precision values and suitable for automatic calculation, is not known at present. Up to now a variety of different methods has been used, according to the topology and the values of masses and momenta of the particular diagram considered.

Therefore, to face problems like four-loop $`g`$-$`2`$, it arises the need of a completely automated approach to the calculation of Feynman integrals, and of finding out suitable new methods, algorithms and techniques of calculations which allow one to solve or avoid the above difficulties.

In this paper we describe the methods and techniques to be used in such automated approach. In particular we present:

1. A new method for determining the master integrals and for reducing generic Feynman integrals to master integrals, applicable to arbitrary Feynman diagrams.
2. A new method for calculating master integrals based on the numerical solution of the recurrence relations provided by the integration-by-parts method, seen as linear difference equations in one index, applicable to arbitrary Feynman diagrams.
3. A new method of generation and solution of systems of differential equations in masses and momenta for master integrals, applicable to arbitrary Feynman diagrams.

The part more innovative and important of this work is the method of calculation of master integrals based on *difference equations*. As this mathematical topic appears (surprisingly) to be practically absent from the literature, in order to improve the intelligibility of the paper we will give an extensive discussion on this argument, including methods of solutions, techniques of calculations and examples of applications. Moreover, study of boundary conditions of difference equations will lead us to a detailed discussion of asymptotic behaviour of Feynman integrals *for large powers of one denominator*, another topic which has received very scarce attention in the literature.

Another important point is the automation of the approach. All methods, algorithms and techniques developed in this work have been implemented in a comprehensive program called $`\mathrm{𝚂𝚈𝚂}`$, written on purpose, which, among other things, contains a simplified algebraic manipulator. The aim of this program is to calculate the value of a generic Feynman integral in completely automatic way, by reducing it to master integrals and by calculating the master integrals. The program turned out to be able to deal with diagrams with hundreds of master integrals and to obtain very high precision values (for example 100-200 digits) expanded at will in $`ϵ=(4D)/2`$, with all coefficients in numerical form, divergent terms included. Limitations of the current implementation will be described.

A crucial point is the test of the approach. By using $`\mathrm{𝚂𝚈𝚂}`$ we have analyzed vacuum, self-energy, vertex and box diagrams up to three or two loops and we have calculated the values of master integrals for some values of masses and momenta, comparing them with already known results, when possible. We will show the results obtained. In particular we have calculated all the previously unknown single-scale massive three-loop self-energy master integrals. These numerical calculations, involving the calculation of hundreds of master integrals at the same time, allowed us to prove the reliability of the approach in real cases and to accumulate a considerable experience of calculation.

The plan of the paper is as follows: In section 2 we illustrate an algorithm for the construction and solution of systems of integration-by-parts identities which allows one to reduce a generic Feynman integral to a combination of master integrals; this algorithm is the basis for the algorithms developed in the following sections. In section 3 we show how to obtain systems of difference equations in one variable satisfied by the master integrals. In section 4 we show how to find solutions of difference equations using expansions in factorial series. In section 5 we discuss the methods to find the values of the arbitrary constants appearing in the solutions of the difference equations. In section 6 we describe the techniques used to sum the factorial series expansions. In section 7 we describe in detail the application of our methods to the one-loop self-energy diagram. In section 8 we illustrate an alternative method of solution of difference equations based on the Laplace’s transformation. In section 9 we apply our methods to various diagrams from one to three loops. In section 10 we show how to deal with integrals with some particular values of masses and momenta, by solving systems of differential equations in masses and momenta. In section 11 we show some technical information on the computer program SYS used in the calculations. In section 12 we give our conclusions.

## 2 Systems of identities between Feynman integrals

After the introduction of the notation used for integrals and identities, in sections 2.3-2.4 we describe the new method used for solving the systems of integration-by-parts identities and for obtaining reduction to master integrals. The description is much detailed because the algorithm here shown is an essential part of this work, being the basis of all similar algorithms of solution of systems of identities which will be described in sections 3.2, 8.2 and 10.2. In section 2.5 we show how the algorithm can be used for determining the master integrals. In sections 2.6-2.7 we show how to use the algorithm for reducing integrals, and we describe its more interesting feature: the absence of intermediate integrals with increasing exponents of the denominators in the intermediate steps of the reduction to master integrals.

### 2.1 Generalities

Let us introduce some notations used throughout the paper. We consider a generic Feynman diagram with $`N_k`$ loops, $`N_e`$ external and $`N_d`$ internal lines. The loop momenta are $`k_i`$, $`i=1,\mathrm{},N_k`$; the independent external momenta are $`p_i`$, $`i=1,\mathrm{},N_p`$, where $`N_p=N_e1`$ (if $`N_e>0`$) for total momentum conservation. The denominators of the propagators are $`D_i=q_i^2+m_i^2`$, where $`q_i`$ is the momentum flowing in the internal line $`i`$ and $`m_i`$ is the mass of the internal line; each $`q_i`$ is a linear combination of the momenta $`\{p_j\}`$ and $`\{k_j\}`$. The usual prescription $`q_i^2+m_i^2=q_i^2+m_i^2i0`$ is implied if necessary. A generic Feynman integral regularized in $`D`$-dimensional euclidean momentum space has the form

$$[d^Dk_1][d^Dk_2]\mathrm{}[d^Dk_{N_k}]V_{\gamma \delta },$$
(1)
where $`[d^Dk_i]=d^Dk_i/\pi ^{D/2}`$ and $`V_{\gamma \delta }`$ is the generic integrand

$$V_{\gamma \delta }=\frac{_{i=1}^{N_p}_{j=1}^{N_k}(p_ik_j)^{\delta _{ij1}}_{i=1}^{N_k}_{j=i}^{N_k}(k_ik_j)^{\delta _{ij2}}}{_{i=1}^{N_d}D_i^{\gamma _i}},\gamma _i0,\delta _{ijl}0,$$
(2)
$`\gamma =\{\gamma _1,\mathrm{},\gamma _{N_d}\}`$ and $`\delta =\{\delta _{ijl}\}`$. The numerator of Eq.(2) is a product of powers of all the possible scalar products involving the loop momenta $`k`$; the total number $`N_{sp}`$ of such scalar products is

$$N_{sp}=N_pN_k+N_k(N_k+1)/2.$$
(3)
### 2.2 Algebraic and integration-by-parts identities

Let us consider the generic integrand (2). For each denominator $`D_j`$ we write the identity

$$\frac{(pk)_j}{D_j}=\frac{1}{C_j}\left(1\frac{D_jC_j(pk)_j}{D_j}\right),j=1,\mathrm{},N_d,$$
(4)
where $`(pk)_j`$ indicates one scalar product involving loop momenta which appears in the expression of $`D_j`$, and $`C_j`$ is its coefficient in the expression. The scalar products $`(pk)_j`$ must be chosen all different. These algebraic identities are applied in sequence to $`V_{\gamma \delta }`$ and to the terms subsequently generated, more times if it is necessary. As a result, the original integrand $`V_{\gamma \delta }`$ is transformed into a sum of new terms, each one not containing $`(pk)_j`$ and $`D_j`$ simultaneously, with the general form

$$V_{ni\alpha \beta }^{}=\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _j}}{_{j=1}^nD_{i_j}^{\alpha _j}},nN_d,\alpha _j,\beta _j0,$$
(5)
where the subscript $`ni\alpha \beta `$ shows the dependence on the number $`n`$ of denominators, their particular combination $`i=\{i_1,\mathrm{},i_n\}`$, and the exponents $`\alpha =\{\alpha _1,\mathrm{},\alpha _n\}`$ and $`\beta =\{\beta _1,\mathrm{},\beta _{N_{sp}n}\}`$. The symbols $`(pk\text{ irred.})_j`$, $`j=1,\mathrm{},N_{sp}n`$ indicate the $`N_{sp}n`$ ‘irreducible’ scalar products which cannot be simplified further on by Eq.(4) with one of the denominators $`D_{i_1}`$, $`\mathrm{}`$, $`D_{i_n}`$ of $`V_{ni\alpha \beta }^{}`$; if $`n=N_{sp}`$ the numerator of Eq.(5) is unity. We stress that one different set of irreducible scalar products corresponds to each different combination of denominators.

Integrating by parts in $`D`$ dimensions  one can write the identities

$$\begin{array}{cc}& [d^Dk_1]\mathrm{}[d^Dk_{N_k}]\frac{}{(k_j)_\mu }\left((p_l)_\mu V_{ni\alpha \beta }^{}\right)=0,j=1,\mathrm{},N_k,l=1,\mathrm{},N_p,\hfill \\ & [d^Dk_1]\mathrm{}[d^Dk_{N_k}]\frac{}{(k_j)_\mu }\left((k_l)_\mu V_{ni\alpha \beta }^{}\right)=0,j,l=1,\mathrm{},N_k,\hfill \end{array}$$
(6)
where $`V_{ni\alpha \beta }^{}`$ is defined in Eq.(5). For each different $`V_{ni\alpha \beta }^{}`$ Eq.(6) gives $`N_k(N_p+N_k)`$ identities. The ratio $`V_{ni\alpha \beta }^{}`$ contains only irreducible scalar products (relative to the particular combination of denominators); the calculation of the derivative and the contraction of the index $`\mu `$ form terms also containing reducible scalar products, which must be transformed into irreducible scalar products using the algebraic identities (4). As a final result, the identity will contain a linear combination of integrals of two kinds: integrals containing all the $`n`$ denominators $`\{D_{i_1}`$, $`D_{i_2}`$, $`\mathrm{}`$, $`D_{i_n}\}`$ appearing in $`V_{ni\alpha \beta }^{}`$ and integrals with one denominator missing as effect of algebraic identities.

### 2.3 A new method for solving the system of identities

Each identity obtained from Eq.(6) is a linear combination of integrals like

$$[d^Dk_1]\mathrm{}[d^Dk_{N_k}]V_{ni\alpha \beta }^{},$$
(7)
with polynomials of degree zero or one in the number of dimensions $`D`$ as coefficients, and where the integrand has the general form of Eq.(5). A generic set of these identities forms a homogeneous linear system of equations, with the integrals as unknowns. As it is well-known, the system is under-determined, and some integrals exist, the “master integrals”, whose values cannot be determined from the system.

In the literature333 For a review see . In the program BUBBLES reduction of systems of identities is performed by semi-automatic means. A different approach to solving integration-by-parts identities was developed in .  it is usual to look for the general solution of the *infinite* system (6) in the form of combination of identities which, lowering and raising exponents, transform integrals of the form (1) or (7) into linear combinations of carefully chosen master integrals. Up to now, the methods used to find these identities have been based on laborious human analysis; as the number of these identities (and the difficulties encountered) grow rapidly with the number of denominators and loops, it becomes very difficult to analyze in this way diagrams beyond a certain limit.

On the contrary, our approach is different, as it consists in the solution of systems made up of a *finite* number of identities. The identities (6) are generated explicitly using suitable $`V_{ni\alpha \beta }^{}`$, with parameters $`n`$, $`i`$, $`\alpha `$, $`\beta `$ taken from a large (but finite) set of values carefully chosen. The set of the generated identities forms a linear system with integrals as unknowns, which is solved using the well-known Gauss elimination method; the solution gives the expressions of the integrals as linear combinations of the master integrals with coefficients rational in $`D`$. The advantage of this approach is that it is simple, applicable without modifications to diagrams with any topology, and suitable for completely automatic calculations.444 A similar method, in very preliminary form, was developed by the author in  to reduce to master integrals the triple-cross vertex diagrams contribution to the $`g2`$ of the electron. In that work a system of about 100000 identities was built and solved, and this required some months of computer time; with our new method the same calculation may be performed in a fraction of hour. Some similar technique has been recently used in , to solve small systems of identities in order to obtain differential equations for two-loop massless box diagrams. It does not consist only in a mere “mechanical” approach because, as we will see in section 2.6, the solution of the whole system allows us to discover very useful identities of a kind *a priori* not expected; these identities allow one to avoid the appearing of intermediate integrals with increasing exponents of the denominators in the intermediate steps of the reduction to master integrals.

Let us now consider more in detail the solution of the system. The identities are generated and inserted in the system one at a time. Let $`_jc_jW_j=0`$ be an identity obtained from Eq.(6), where $`W_j`$ are integrals and $`c_j`$ the coefficients. The identities already existing in the system, expressing some of the integrals $`W_j`$ in terms of other integrals are substituted in the new identity, which becomes $`_jc_j^{}W_j^{}=0`$. One particular integral $`W_l^{}`$ is chosen between the integrals $`\{W_j^{}\}`$, and the identity itself is rewritten as $`W_l^{}=_{jl}c_j^{\prime \prime }W_j^{}`$ in order to express the particular integral in terms of the remaining integrals. Then the new identity is added to the system and the chosen integral $`W_l^{}`$ is substituted in the rest of the system.

The choice of the integral $`W_l^{}`$ is carried out following an ordering555 The importance of a total ordering for the solution of systems of differential equations derived from integration-by-parts identities was remarked in . of all the integrands $`V_{ni\alpha \beta }^{}`$ as functions of the parameters: the number $`n`$ of denominators, the combination $`i`$ of denominators, and the exponents. The ordering defines the priority of extraction between integrals. Priorities are arranged so that integrals with a higher number of denominators are extracted first, and expressed as linear combinations of the integrals with a lower number of denominators. Remaining details of ordering used are shown in the algorithm of the next section. The form of the master integrals depends on the choice of the ordering: see section 2.5.

Now we consider the choice of the order of the generation and processing of the identities. The final solution of the system is independent of the order of processing, but the computing time is not. Each addition of a new identity to the system implies a substitution of an integral in all identities of the system which contains it; therefore the order must be carefully chosen in order to minimize the number of substitutions required. A bad choice may cause the computing time to blow up.

A good choice of the ordering of the ratios $`V_{ni\alpha \beta }^{}`$ appearing in Eq.(6) is the inverse of the above considered ordering of integrands used for the extraction of integrals. The identities corresponding to ratios $`V_{ni\alpha \beta }^{}`$ with the lowest priority of extraction (as integrands), therefore with the lowest number of denominators, are generated and processed first; then the identities with a higher number of denominators are processed. The last processed identities will have the highest number of denominators. This choice is efficient as long as the number $`n`$ of denominators of the ratio $`V_{ni\alpha \beta }^{}`$ used in Eq.(6) is equal to the number $`n^{}`$ of denominators of the integral $`W_l^{}`$ extracted from the identity generated. In some cases (relatively rare, but of great importance, see section 2.6) one finds $`n^{}<n`$; the identities obtained with ratios $`V_{n^{}i^{}\alpha ^{}\beta ^{}}^{}`$, containing integrals with $`n^{}`$ denominators, have already been inserted into the system, so that the “late” integral $`W_l^{}`$ must be substituted in a large number of identities, and that may require a considerable amount of time.

It is useful to split the system of identities (6) into several *subsystems*. Each subsystem is made up of all the identities obtained by inserting in Eq.(6) terms $`V^{}`$ containing one particular combination $`i`$ of denominators $`\{D_{i_1}`$, $`D_{i_2}`$, …, $`D_{i_n}\}`$. The integrals with a number of denominators $`n`$ less than the number of loops $`N_k`$ of the diagram are null in the framework of the dimensional regularization, and must be not considered, so that the total number $`S`$ of subsystems is

$$S=\underset{n=N_k}{\overset{N_d}{}}\left(\genfrac{}{}{0pt}{}{N_d}{n}\right)<2^{N_d}.$$
(8)
Let us now define the non-negative quantities $`M_p`$, $`M_d`$,

$$\underset{j=1}{\overset{N_{sp}n}{}}\beta _j=M_p,\underset{j=1}{\overset{n}{}}\left(\alpha _j1\right)=M_d,$$
(9)
which are the total sum of the powers of the scalar products in the numerator and the sum of the powers of the denominators minus the number of denominators of the generic integrand $`V_{ni\alpha \beta }^{}`$ of Eq.(5), respectively. Using this definition we can split the infinite set of the integrands $`V_{ni\alpha \beta }^{}`$ into finite sets of integrands with equal $`n`$, $`M_p`$ and $`M_d`$ that we indicate with the symbol $`[n;\frac{[}{0}pt]M_pM_d]`$. The set of the integrands with $`0M_pa`$ and $`0M_db`$ will be indicated with the symbol $`[n;\frac{[}{0}pt]0\mathrm{}a0\mathrm{}b]`$. In the following sections we will say in short that an integral (7) or an identity (6) belongs to a specified set, understanding that it is the corresponding $`V_{ni\alpha \beta }^{}`$ which belongs to that set.

### 2.4 The algorithm of solution of the system

The algorithm of solution is the following:

###### Algorithm 1

1. Let $`n=N_k`$.
2. Let $`i_1=1`$, $`i_2=2`$, …, $`i_n=n`$.
3. Consider the combinations of $`n`$ different denominators, chosen in the set $`\{D_1`$, $`\mathrm{}`$, $`D_{N_d}\}`$; let $`D_{i_1},D_{i_2},\mathrm{},D_{i_n}`$ be one of these combinations.
4. Choose two integer non-negative constants $`a_i`$ and $`b_i`$, $`i=\{i_1,\mathrm{},i_n\}`$.
5. Let $`M_d=0`$.
6. Let $`M_p=0`$.
7. Consider the ratios of the kind of Eq.(5), containing the $`n`$ denominators

$`D_{i_1},D_{i_2},\mathrm{},D_{i_n}`$ raised to some power; let $`W`$ be one of these ratios:

$$W(n,i,\alpha ,\beta )=\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _j}}{D_{i_1}^{\alpha _1}D_{i_2}^{\alpha _2}\mathrm{}D_{i_n}^{\alpha _n}},$$

where the non-negative exponents $`\alpha _j`$, $`\beta _j`$ are such that

$$\underset{j=1}{\overset{N_{sp}n}{}}\beta _j=M_p,\underset{j=1}{\overset{n}{}}\left(\alpha _j1\right)=M_d,$$

i.e., they are such that $`W`$ belongs to the set $`[n;\frac{[}{0}pt]M_pM_d]`$.
8. Set $`V_{ni\alpha \beta }^{}=W`$ in Eq.(6) and generate all the integration-by-parts identities.
9. Let $`[d^Dk_1]\mathrm{}[d^Dk_{N_k}]_jc_jW_j=0`$ be one of the identities of Eq.(6), where $`W_j`$ is a generic integrand with $`n`$ or $`n1`$ denominators, then:

1. Substitute all the integrals already known in the left-hand side of the new identity;

    let $`[d^Dk_1]\mathrm{}[d^Dk_{N_k}]_jc_j^{}W_j^{}=0`$ be the result. If the new identity is a linear combination of other identities of the system, go to step 10.
    2. If the new identity is linearly independent, choose an integrand $`W_l^{}`$ from the identity,

    $$W_l^{}(n^{},i^{},\alpha ^{},\beta ^{})=\frac{_{j=1}^{N_{sp}n^{}}(pk\text{ irred.})_j^{\beta _j^{}}}{D_{i_1^{}}^{\alpha _1^{}}D_{i_2^{}}^{\alpha _2^{}}\mathrm{}D_{i_n^{}}^{\alpha _n^{}}},$$

    belonging to $`[n^{};\frac{[}{0}pt]M_p^{}M_d^{}]`$, $`N_kn^{}n`$, following the order of priority:

    1. the greatest number of the denominators $`n^{}`$;
        2. the greatest $`M_d^{}`$;
        3. the greatest $`M_p^{}`$;
        4. the greatest $`i_1^{}`$, the greatest $`i_2^{}`$, …, the greatest $`i_n^{}`$;
        5. the greatest $`\alpha _1^{}`$, the greatest $`\alpha _2^{}`$, …, the greatest $`\alpha _n^{}`$;
        6. the greatest $`\beta _1^{}`$, the greatest $`\beta _2^{}`$, …, the greatest $`\beta _{N_{sp}n^{}}`$.
    3. Substitute and add the following identity to the system:

    $$[d^Dk_1]\mathrm{}[d^Dk_{N_k}]W_l^{}=[d^Dk_1]\mathrm{}[d^Dk_{N_k}]\underset{jl}{}(c_j^{}/c_l^{})W_j^{}.$$
10. Generate a new integration-by-parts identity among the $`N_k(N_p+N_k)`$ possible identities of Eq.(6) and go to step 9, otherwise continue.
11. Choose a new integrand $`W`$ with different exponents $`\alpha `$ and $`\beta `$, belonging to the set $`[n;\frac{[}{0}pt]M_pM_d]`$, among the $`\left(\genfrac{}{}{0pt}{}{N_{sp}n+M_p1}{M_p}\right)\left(\genfrac{}{}{0pt}{}{n+M_d1}{M_d}\right)`$ elements of the set and go to step 8, otherwise continue.
12. $`M_p=M_p+1`$; if $`M_pa_i`$ go to step 7.
13. $`M_d=M_d+1`$; if $`M_db_i`$ go to step 6.
14. Choose a new combination of indices $`i_1<i_2<\mathrm{}<i_n`$ among the $`\left(\genfrac{}{}{0pt}{}{N_d}{n}\right)`$ possible combinations of numbers $`\{1`$, $`2`$, $`\mathrm{}`$, $`N_d\}`$ and go to step 3, otherwise continue.
15. $`n=n+1`$; if $`nN_d`$ go to step 2, otherwise end.

In this algorithm two arbitrary integer constants, $`a_i`$ and $`b_i`$, must be chosen for each different combination of denominators $`\{D_{i_1},\mathrm{},D_{i_n}\}`$. These constants are cutoffs and define which identities must be included in the system and which must be excluded; $`a_i`$ limits the sum of the powers of scalar products in the numerator and $`b_i`$ limits the sum of the exponents of the denominators.

With a suitable choice of the parameters $`a_i`$ and $`b_i`$ (see section 2.6), this algorithm allows one to reduce any given integral $`I`$ to a sum of master integrals $`B_l`$

$$I=\underset{l=1}{\overset{L}{}}r_lB_l,$$
(10)
where $`r_l`$ are rational functions of $`D`$, masses and scalar products of external momenta.

### 2.5 Determination of master integrals

Our algorithm of solution of the system of identities provides a general method for determining the master integrals: it suffices, for each different combination of $`n`$ denominators $`\{D_{i_1}`$, $`D_{i_2}`$, …, $`D_{i_n}\}`$ to build the subsystem with suitable $`a_i`$ and $`b_i`$, to solve it using the algorithm, and to find the integrals belonging to $`[n;\frac{[}{0}pt]0\mathrm{}a_i0\mathrm{}b_i]`$ which are not reduced to terms with $`n1`$ denominators. These should be master integrals. Of course, since we are dealing with a system made up with a finite number of identities, there is the possibility that a master integral is erroneously identified or not found because some essential identities have not been included in the system as they exceed the limits of generation. Explicit solutions of subsystems with increasing values of $`a_i`$ and $`b_i`$ show that empirical ‘thresholds’ $`\overline{a}_i`$ and $`\overline{b}_i`$ exist, depending on the combination of denominators, such that all the subsystems built with $`a_i\overline{a}_i`$ and $`b_i\overline{b}_i`$ give a stable and correct identification of the master integrals; for the diagrams so far examined (see section 9) one finds $`0\overline{a}_i3`$ (it is equal to the maximum number of scalar products appearing in the master integrals) and $`\overline{b}_i=0`$. Corresponding number of identities varies from some tens to some hundreds.

The particular order of selection of the integrals used in step 9b of the algorithm favours for each different combination of denominators the choice of master integrals with all the denominators $`D_j`$ raised to the first power and with the numerator containing increasing numbers of scalar products: first the scalar integral, with 1 as numerator, then integrals with one scalar product, and so on.

As result of the procedure just described we find the set of $`L`$ master integrals $`B_l`$ with the form

$$B_l=[d^Dk_1]\mathrm{}[d^Dk_{N_k}]\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _{jl}}}{D_{i_1}D_{i_2}\mathrm{}D_{i_n}},$$
(11)
where the combination of indices $`i`$ and exponents $`\beta `$ depends on the index $`l`$. The ordering $`B_1`$, $`B_2`$, …, $`B_L`$ follows the ordering of integrals provided by the algorithm.

### 2.6 Choice of constants $`a_i`$ and $`b_i`$: the “golden rule”

The minimal values of $`a_i`$ and $`b_i`$ necessary to reduce one particular integral strongly depend on the structure of the integral and are not easy to determine. Therefore we limit ourselves to give some general rules, based on the experience acquired by using the algorithm, rather than rigorous results. Actually, these rules turned out to be very effective.

Let us consider the system of identities obtained by inserting in Eq.(6) all the terms $`V_{ni\alpha \beta }^{}`$ belonging to the set $`G_{ab}`$

$$G_{ab}=\underset{n=N_k}{\overset{N_d}{}}[n;\frac{[}{0}pt]0\mathrm{}a0\mathrm{}b]=[N_k\mathrm{}N_d;\frac{[}{0}pt]0\mathrm{}a0\mathrm{}b]$$
(12)
containing all the possible $`S`$ subsystems of identities (see Eq.(8)), and choosing $`a_i=a`$ and $`b_i=b`$ for all the combinations of denominators. The number of elements $`V_{ni\alpha \beta }^{}`$ forming the set $`G_{ab}`$ is given by

$$N_{ele}\left(G_{ab}\right)=\underset{n=N_k}{\overset{N_d}{}}\left(\genfrac{}{}{0pt}{}{N_d}{n}\right)\left(\genfrac{}{}{0pt}{}{N_{sp}n+a}{a}\right)\left(\genfrac{}{}{0pt}{}{n+b}{b}\right),$$
(13)
while the number of identities of the system is

$$N_{ide}\left(G_{ab}\right)=N_k\left(N_p+N_k\right)N_{ele}(G_{ab}).$$
(14)
Explicit calculations suggest the following empirical “golden rule”: the solution of the system provides identities which reduce any integral belonging to the set $`G_{ab}`$ to master integrals

$$\text{Integrals}G_{ab}\mathrm{@}>>\text{Identities }G_{ab}>\text{Master Integrals}.$$
(15)
In order to better understand this rule, we have analyzed in some test cases the reduction of single integrals belonging to $`[N_d;\frac{[}{0}pt]0\mathrm{}a0\mathrm{}b]`$, by modifying the steps 1 and 15 of the algorithm so that the various subsystems of identities $`[n;\frac{[}{0}pt]0\mathrm{}a0\mathrm{}b]`$ are generated starting from $`N_d`$ denominators and decreasing $`n`$ up to $`N_k`$ denominators (note that the solution of the whole system or parts of it with such ‘inverted’ order is extremely time consuming). For each value of $`n`$, intermediate expressions contain master integrals with a number of denominators from $`n`$ to $`N_d`$, and integrals belonging to the set $`[n1;\frac{[}{0}pt]0\mathrm{}a0\mathrm{}b+1]`$, that is, integrals with one exponent of denominators increased by one compared with the original integral *whichever is $`n`$*. This is very different from what one obtains by performing the reduction one denominator at a time with recurrence relations which lower one index and raise another one, where the exponents of denominators of the integrals are observed to increase whenever the number of denominators decreases, causing the number of intermediate integrals to blow up. Our algorithm avoids this blowing-up.

The quite unexpected behaviour of our algorithm seems to be due to the fact that solving the whole system we consider all the identities generated, not only those which reduce single integrals, but also those many “additional” identities apparently useless which reduce *combinations* of integrals with the same number of denominators. Their effect is to reduce systematically all the generated combinations of integrals with increased exponents of the denominators. The reason for this behaviour is not known; however, its systematic presence suggests that (hopefully) it may have some simple explanation. The presence of these “additional” identities is important, particularly if $`N_d`$ is large, because without them the number of identities to consider (and the number of intermediate integrals) would be orders of magnitude larger.

The “golden rule” (15) turns out to be valid if $`aa_0`$ and $`bb_0`$, where $`a_0`$ and $`b_0`$ are some empirical ‘thresholds’ depending on the structure of denominators of the diagram; for the diagrams examined in section 9 we found $`1a_03`$ (values almost always equal to the maximum number of scalar products appearing in the master integrals) and $`b_0=0,1`$. Values of the thresholds are unfortunately not known *a priori*; if the chosen values of $`a`$ and $`b`$ are below the threshold, the result of the solution of the system is that some integrals belonging to $`G_{ab}`$ are not completely reduced to master integrals, and some non-master integrals with a few denominators still survive in the final expressions. In this case one may increase $`a`$ or $`b`$, or alternatively one may suitably enlarge only the subsystems of identities corresponding to residual integrals.

We found some very rare exceptions to the “golden rule” (15), when some denominators have zero mass; in these cases, whatever the values of $`a`$ and $`b`$ are, there are always few integrals belonging to $`[N_d;\frac{[}{0}pt]0b]`$ (note, changing with the value of $`b`$) which are not completely reduced to master integrals. In the example of the next section we will show just one of these exceptions.

Now, let us consider an important application: the reduction to master integrals of one combination of many integrals, for example, expressing the contribution of a diagram to some physical quantity. Let us suppose that all the integrals of this combination belong to a set $`G_{\overline{a}\overline{b}}`$ with some (minimal) $`\overline{a}`$ and $`\overline{b}`$; if the values of $`\overline{a}`$ and $`\overline{b}`$ are over thresholds, according Eq.(15) the solution of the system will provide identities which reduce to master integrals all the integrals of the set.

### 2.7 An example

Let us consider the self-energy diagram with 5 denominators $`D_1=(pk_1)^2+1`$, $`D_2=(pk_1k_2)^2+1`$, $`D_3=(pk_2)^2+1`$, $`D_4=k_1^2`$, $`D_5=k_2^2`$ and $`p^2=1`$. Following the notation of section 2.1, this diagram has $`N_p=1`$, $`N_k=2`$, $`N_d=5`$ and $`N_{sp}=5`$. For instance, we want to transform the integral

$$J=\frac{[dk]}{D_1^2D_2D_3D_4D_5},$$
(16)
where $`[dk]=[d^Dk_1][d^Dk_2]`$, into a combination of master integrals. The pairs between scalar products and denominators used in the algebraic identities are $`(pk_1,D_1)`$, $`(k_1k_2,D_2)`$, $`(pk_2,D_3)`$, $`(k_1k_1,D_4)`$ and $`(k_2k_2,D_5)`$. We must identify the master integrals; therefore, following section 2.5, we build the system of identities with $`a_i=1`$ and $`b_i=0`$ for all the combinations of denominators, and we look for the integrals which are not reduced. The master integrals turn out be: $`B_{123}=[dk]/D_1D_2D_3`$, $`B_{345}=[dk]/D_3D_4D_5`$, $`B_{12}=[dk]/D_1D_2`$, $`B_{13}=[dk]/D_1D_3`$ and $`B_{23}=[dk]/D_2D_3`$.

We consider the whole system made up of the identities of the set $`G_{11}`$, choosing $`a=b=1`$; according to Eq.(14) it contains 1776 identities. Solving it with the algorithm 1 (by using the program SYS described in section 11) we find 1122 independent identities. Examining them, we find the identities which reduce to master integrals all the 291 integrals $`[2\mathrm{}5;\frac{[}{0}pt]0\mathrm{}10\mathrm{}1]`$, integral $`J`$ included, with the exception of the two integrals $`[dk]/D_1D_2D_3D_4^2D_5`$ and $`[dk]/D_1D_2D_3D_4D_5^2`$, whose reduction still contains the integrals $`[dk]/D_2D_4^2D_5^2`$ and $`[dk]/D_3D_4^2D_5^2`$; in order to complete the reduction of these two integrals, we may add the two sets of 30 identities $`[3;\frac{[}{0}pt]0\mathrm{}10\mathrm{}2]`$ with these combinations of three denominators to the system. The remaining 831 identities reduce to master integrals complicated combinations of integrals $`[2\mathrm{}5;\frac{[}{0}pt]0\mathrm{}22]`$.

We want to illustrate how the mechanism of “additional” identities works. Let us consider the partial reduction of the integral $`J`$ from 5 to 4 denominators. We consider the subsystem of the identities $`[5;\frac{[}{0}pt]00\mathrm{}1]`$, which is formed by 36 identities; solving it, we find:

1. 6 identities reducing the integrals $`[5;\frac{[}{0}pt]00\mathrm{}1]`$ to 4 denominators;
2. 9 identities reducing combinations of integrals $`[5;\frac{[}{0}pt]02]`$ to 4 denominators;
3. 20 identities between integrals $`[4;\frac{[}{0}pt]01]`$ and $`[4;\frac{[}{0}pt]02]`$ with different combinations of denominators;

The identities of the groups 2 and 3 are the “additional” ones. Among the identities of the group 1 we find

$$J=[dk]\left(\frac{1}{4}W_1+\frac{1}{2}W_2\right),$$
(17)
$$\begin{array}{cc}\hfill W_1& =\frac{3}{D_2^2D_3D_4D_5}\frac{3}{D_1D_2^2D_3D_4}+\frac{2}{D_2D_3D_4^2D_5}\frac{1}{D_1D_2D_3^2D_4}+\hfill \\ & +\frac{1}{D_1D_2D_3^2D_5}+\frac{1}{D_2D_3D_4D_5^2}\frac{1}{D_1D_3^2D_4D_5}\frac{1}{D_1D_3D_4D_5^2},\hfill \\ \hfill W_2& =\frac{1}{D_1D_3^2D_4^2D_5}\frac{1}{D_1D_2D_4^2D_5^2}\frac{1}{D_1D_2^2D_4^2D_5}+\frac{1}{D_1D_3D_4^2D_5^2}\frac{1}{D_2D_3D_4^2D_5^2}.\hfill \end{array}$$

We see that the original integral, which contains a square denominator (set $`[5;\frac{[}{0}pt]01]`$), is transformed into a combination of integrals with 4 denominators, with one or two square denominators (sets $`[4;\frac{[}{0}pt]01]`$ and $`[4;\frac{[}{0}pt]02]`$).

Considering the reduction from 4 denominators to master integrals, if we peruse the solution of the subsystems $`[2\mathrm{}4;\frac{[}{0}pt]0\mathrm{}10\mathrm{}1]`$ we find 8 relations which transform each one of the terms of $`W_1`$ into master integrals; on the contrary we do not find relations which transform *single* terms of $`W_2`$ into master integrals. But if we also consider the previous identities of the group 3, we find one relation which transforms the *whole* combination $`W_2`$. The presence of this relation makes it unnecessary to search for identities which singly transform the integrals of $`W_2`$.

## 3 Difference equations for master integrals

By using the algorithm shown in the above sections, a generic Feynman integral can be reduced to master integrals. Now we consider their calculation.

We consider the generic integral (7) as a function of the exponents $`U(\alpha ,\beta )`$; the integration-by-parts identities (6) form a system of recurrence relations satisfied by the function $`U`$. The key observation is that such recurrence relations are *linear difference equations*666 Other authors are also encountering other kinds of difference equations in the evaluation of diagrams  and solving them, by using techniques different from the standard methods of  used in this paper.  satisfied by the function $`U(\alpha ,\beta )`$. This observation has two important implications:

1. Theory of linear difference equations is a well established (but not very well-known) mathematical topic for which various useful mathematical tools exist.
2. Difference equations can be solved numerically, and therefore this establishes a new method of calculation of the integrals $`U(\alpha ,\beta )`$.

Another important observation is that the integration-by-parts provides recurrence relations in *several* variables for the functions $`U`$, that is, a system of partial difference equations. In general, the numerical solution of such a system is, from the point of view of the numerical calculation, a task comparable to the solution of a system of partial derivative equations or to the multidimensional integration, so that high-precision results may be very difficult to obtain if the number of variables is high. Therefore, if we desire to obtain high-precision results, we must consider only difference equations in only *one* variable.

In section 3.1 we show where to introduce such single parameter; in section 3.2 we show how to obtain recurrence relations (that is, difference equations) in this parameter from the integration-by-parts identities (6).

### 3.1 Difference equations in one exponent

Let us consider, for instance, a master integral of the form

$$B=\frac{[d^Dk_1]\mathrm{}[d^Dk_{N_k}]}{D_1D_2\mathrm{}D_{N_d}}.$$
(18)
We raise one denominator to power $`x`$. If we choose $`D_1`$, we define a function $`U_{D_1}`$

$$U_{D_1}(x)=\frac{[d^Dk_1]\mathrm{}[d^Dk_{N_k}]}{D_1^xD_2\mathrm{}D_{N_d}},$$
(19)
where the subscript shows that the exponent has been introduced in the denominator $`D_1`$. The value of the integral (18) is recovered as $`U_{D_1}(1)=B`$.

Using integration-by-parts identities, we look for an identity which relates $`U_{D_1}`$ with integrals with a smaller number of denominators. As $`B`$ is a master integral, an identity which directly transforms $`U_{D_1}(x)`$ into simpler integrals cannot be found; instead, as we will see in the next section, one finds identities of the form

$$\underset{j=0}{\overset{R}{}}p_j(x)U_{D_1}(x+j)=F(x),$$
(20)
where $`p_j(x)`$ are polynomials in $`x`$, and $`F(x)`$ is a known function. $`U_{D_1}`$ appears with arguments shifted by an integer. This identity is a linear difference equation of order $`R`$ in the variable $`x`$ satisfied by the function $`U_{D_1}(x)`$.

The right-hand side function $`F(x)`$ is in general a linear combination of functions analogous to $`U_{D_1}`$, derived from master integrals containing $`D_1`$, but with some of the denominators $`D_2,D_3,\mathrm{},D_{N_d}`$ missing.

We note that it is possible to raise to $`x`$ another denominator, $`D_2`$, $`D_3`$, …, etc. (or equivalently to permute the numbering of the lines of the diagram); for each denominator $`D_j`$ a function $`U_{D_j}(x)`$ will be defined. In general the functions $`U_{D_j}(x)`$ are different each one and satisfy different difference equations, but for $`x=1`$ they all have the same value $`U_{D_j}(1)=B`$. This fact is particularly useful for checking the consistency of the calculations.

### 3.2 Construction of systems of difference equations

Let us suppose that, by means of the procedure described in section 2.5, we establish that there are $`L`$ master integrals, with a number of denominators ranging from $`N_k`$ to $`N_d`$.

We note that some master integrals may have the first denominators missing. Then we group the master integrals in $`N_dN_k+1`$ sets $`S_m`$ such that each set includes the integrals which contain the denominators $`D_i`$ with $`im`$. Each set will contain $`L_m`$ master integrals $`B_{ml}`$,

$$B_{ml}=[d^Dk_1]\mathrm{}[d^Dk_{N_k}]\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _{jml}}}{D_mD_{i_2}D_{i_3}\mathrm{}D_{i_n}},m<i_2<\mathrm{}<i_n,$$
(21)
where $`m=1,\mathrm{},N_dN_k+1`$, $`l=1,\mathrm{},L_m`$ and $`_mL_m=L`$; the number of denominators $`n`$ ranges from $`N_k`$ to $`N_dm+1`$. Note that $`n`$, the indices $`i`$ and the exponents $`\beta `$ of the various integrals depend on $`m`$ and $`l`$. The ordering $`B_{m1}`$, $`B_{m2}`$, …, $`B_{mL_m}`$ follows the ordering of master integrals with the same first denominator $`D_m`$ provided by the algorithm 1. Changing the exponent of $`D_m`$ from 1 to $`x`$ in each master integral, we define the “master functions”

$$U_{ml}(x)=[d^Dk_1]\mathrm{}[d^Dk_{N_k}]\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _{jml}}}{D_m^xD_{i_2}D_{i_3}\mathrm{}D_{i_n}};$$
(22)
$`U_{ml}(1)=B_{ml}`$ recovers the original master integral. Now we must find the difference equations satisfied by the functions $`U_{ml}(x)`$. We follow a method quite similar to that used for the reduction of a generic integral to master integrals. We build a system of identities obtained from integration-by-parts using Eq.(6), modifying the integrand $`V_{ni\alpha \beta }^{}`$ of Eq.(5) by adding $`x1`$ to the exponent of the first denominator:

$$V_{ni\alpha \beta }^{}(x)=\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _j}}{D_{i_1}^{x1+\alpha _1}D_{i_2}^{\alpha _2}\mathrm{}D_{i_n}^{\alpha _n}},nN_d,\alpha _j,\beta _j0.$$
(23)
In this way all the integrals become functions of $`x`$

$$U_{ni\alpha \beta }(x)=[d^Dk_1]\mathrm{}[d^Dk_{N_k}]V_{ni\alpha \beta }^{}(x),$$
(24)
and the identities (6) become difference equations between these functions. Integrals which have a different first denominator $`D_m`$ never appear in the same identity, so that it is convenient to build and solve separately the systems of the identities with different values of $`m`$. We do this with the algorithm

###### Algorithm 2

Consider the algorithm 1 with the following modifications:

1. Set $`i_1=m`$.
2. Add $`x1`$ to the exponent of $`D_m`$.
3. Ignore the addition of $`x1`$ to $`M_d`$ in the definition of set (9).
4. Add the following conditions to the order of priority of the extraction of integrals:

1. the generic functions $`U_{ni\alpha \beta }`$ must be extracted before of the functions $`U_{ml}`$;
    2. the functions $`U_{ml}(x1)`$ (generated by the algebraic identities) must be extracted before of the functions $`U_{ml}(x+i)`$, where $`i0`$.

As result of the solution of the system, we find lots of identities analogous to Eq.(10), not interesting, which express generic functions $`U_{ni\alpha \beta }`$ as combinations of master functions $`U_{ml}`$ with argument possibly shifted, and rational functions of $`x`$, $`D`$, masses and scalar products of external momenta as coefficients; we also find much more interesting identities containing only master functions $`U_{ml}`$.

Let us suppose, for simplicity, that each master integral $`B_{ml}`$ contains a different combination of denominators. With a suitable choice of the constants $`a_i`$ and $`b_i`$, the solution of the system provides difference equations satisfied by the functions $`U_{ml}`$, with the structure

$$\underset{i=0}{\overset{R_l}{}}p_{il}(x)U_{ml}(x+i)=\underset{k=1}{\overset{l1}{}}\underset{j=0}{\overset{Q_{lk}}{}}q_{jkl}(x)U_{mk}(x+j),l=1,\mathrm{},L_m,$$
(25)
where $`p_{il}`$ and $`q_{jkl}`$ are polynomials in $`x`$, $`D`$, masses and scalar products of external momenta. The right-hand side of the $`l`$th equation contains the functions from $`U_{m1}`$ to $`U_{m,l1}`$. Therefore the set of $`L_m`$ equations (25) forms a triangular system of difference equations. The triangular structure is particularly useful for simplifying the numerical solution: the equations are solved in ascending order, one at a time, for $`l=1,2,\mathrm{},L_m`$; when the equation for $`U_{ml}`$ has to be solved, the functions $`U_{m1}`$,…,$`U_{m,l1}`$ in the right-hand side are already known. Note that some functions may be missing from the right-hand side; for example, if the master function $`U_{ml}`$ has $`N_k`$ denominators, the right-hand side is zero.

Now we consider the more general case where $`G`$ different master integrals have the same denominator $`(D_{i_1}\mathrm{}D_{i_n})`$ and different numerators. These master integrals correspond to master functions with contiguous indices $`U_{m,l^{}+1}`$, $`U_{m,l^{}+2}`$, …, $`U_{m,l^{}+G}`$ as effect of the ordering of master integrals. The solution of the system provides for these functions a subsystem of simultaneous $`G`$ difference equations of the kind

$$\underset{g=1}{\overset{G}{}}\underset{i=0}{\overset{R_{hg}}{}}p_{igl^{}h}^{}(x)U_{m,l^{}+g}(x+i)=\underset{k=1}{\overset{l^{}}{}}\underset{j=0}{\overset{Q_{hl^{}k}}{}}q_{jkl^{}h}^{}(x)U_{mk}(x+j),h=1,\mathrm{},G,$$
(26)
not in triangular form, where the left-hand side of each equation contains all the master functions $`U_{m,l^{}+1}`$, …, $`U_{m,l^{}+G}`$. We prefer not deal with the solution of the subsystem of simultaneous difference equations, so that we transform it into triangular form; this is not obligatory, is only convenient to simplify the subsequent numerical solution. We make the replacement $`xx+c`$ in the equations (26), with $`c=1,2,\mathrm{}`$, generating new identities which are inserted in the subsystem, and taking care of extracting the terms containing the function $`U_{m,l^{}+j}`$ before of the terms containing $`U_{m,l^{}+k}`$ if $`j>k`$. This procedure is repeated until one obtains a set of $`G`$ equations in triangular form,

$$\underset{i=0}{\overset{R_h^{}}{}}p_{il^{}h}^{\prime \prime }(x)U_{m,l^{}+h}(x+i)=\underset{k=1}{\overset{l^{}+h1}{}}\underset{j=0}{\overset{Q_{hl^{}k}^{}}{}}q_{jkl^{}h}^{\prime \prime }(x)U_{mk}(x+j),h=1,\mathrm{},G.$$
(27)
Another advantage of the transformation into triangular form is that one obtains a difference equation containing the function $`U_{m,l^{}+1}`$, but not containing the other functions $`U_{m,l^{}+2}`$, $`U_{m,l^{}+3}`$,…, so that $`U_{m,l^{}+1}`$, which corresponds to the master integral with 1 as numerator, can be found independently. Unfortunately the transformation into triangular form has a price: one obtains for $`U_{m,l^{}+1}`$ an equation of order $`R_1^{}>\mathrm{max}_hR_{h1}`$ which has coefficients $`q_{jkl^{}1}^{\prime \prime }(x)`$ much more complicated and cumbersome than the coefficients $`q_{jkl^{}h}^{}(x)`$ of the equations of the subsystem (26). If $`R_1^{}`$ is large (say $`R_1^{}>10`$) these coefficients become huge and difficult to obtain so that it may be more convenient in these cases to solve the system of simultaneous equations (26) directly; this will be described in some next paper.

Once all the subsystems of simultaneous equations corresponding to the various groups of master integrals with equal denominators are transformed into triangular form, the whole system takes the form (25).

Concerning the choice of $`a_i`$ and $`b_i`$, all the considerations done in section 2.6 remain valid; for each combination of $`n`$ denominators $`\{D_{i_1}\mathrm{}D_{i_n}\}`$ there is a minimal subsystem $`[n;\frac{[}{0}pt]0\mathrm{}\stackrel{~}{a}_i0\mathrm{}\stackrel{~}{b}_i]`$ whose solution allows one to obtain the equations (26). For the diagrams so far examined $`\stackrel{~}{a}_i`$ turns out to be always equal to the maximum number of scalar products appearing in the master integrals (typically $`1\mathrm{}3`$); typical values of $`\stackrel{~}{b}_i`$ are $`0,1,2`$.

## 4 Solutions of difference equations with factorial series

By using the algorithms of section 3, the triangular system of difference equations satisfied by the master functions is worked out. Now we consider the solution of each difference equation.

Let us suppose, for instance, that the master function $`U(x)`$ defined as

$$U(x)=[d^Dk_1]\mathrm{}[d^Dk_{N_k}]\frac{_{j=1}^{N_{sp}N_d}(pk\text{ irred.})_j^{\beta _j}}{D_1^xD_2\mathrm{}D_{N_d}};$$
(28)
satisfies a difference equation of order $`R`$

$$\underset{i=0}{\overset{R}{}}p_i(x)U(x+i)=F(x),$$
(29)
where $`p_i(x)`$ are polynomials and $`F(x)`$ is some known function. The solution of this nonhomogeneous equation can be written as

$$U(x)=U^{HO}(x)+U^{NH}(x),$$
(30)
where $`U^{NH}`$ is a particular solution of the nonhomogeneous equation (29) and $`U^{HO}`$ is the general solution of the homogeneous equation

$$\underset{i=0}{\overset{R}{}}p_i(x)U^{HO}(x+i)=0.$$
(31)
The general solution of Eq.(31) can be written as

$$U^{HO}(x)=\underset{j=1}{\overset{R}{}}\stackrel{~}{\omega }_j(x)U_j^{HO}(x),$$
(32)
where $`\stackrel{~}{\omega }_j(x)`$ are periodic functions of period 1, and $`\{U_1^{HO}(x),\mathrm{},U_R^{HO}(x)\}`$ is a fundamental system of independent solutions of the homogeneous equation. In the following, recalling from  for the ease of the unfamiliar reader the essential matter on the solution of linear difference equations with factorial series, we describe how to obtain factorial series expansions of $`U^{HO}`$ and $`U^{NH}`$.

### 4.1 Factorial series

The series

$$\begin{array}{c}\underset{s=0}{\overset{\mathrm{}}{}}\frac{a_s\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(xK+s+1)}=\frac{\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(xK+1)}\left(a_0+\frac{a_1}{xK+1}+\frac{a_2}{(xK+1)(xK+2)}+\mathrm{}\right)\hfill \end{array}$$
(33)
is known as factorial series of the first kind , or series of inverse factorials , or faculty series . We refer to it in brief as factorial series. The series converges for every point in a half-plane which is limited on the left by $`\mathrm{}x=\lambda `$ (excluding $`x=K1,K2,\mathrm{}`$). The number $`\lambda `$ is the *abscissa of convergence*. If $`\lambda =\mathrm{}`$ the series is everywhere divergent.

As coefficients of the series encountered in this work behave for large $`s`$ as $`|a_s|s!s^\alpha `$, it is useful (especially for numerical applications) to define the reduced coefficients $`a_s^{}=a_s/s!`$. For large $`s`$ the generic term of the sum tends to $`a_s^{}\mathrm{\Gamma }(x+1)s^{Kx}`$, so that the factorial series has the same convergence properties as the Dirichlet series $`_{s=0}^{\mathrm{}}a_s^{}s^{Kx}`$.

### 4.2 Operators $`𝝅`$ and $`𝝆`$

Given an arbitrary number $`m`$ one defines the operator $`𝝆`$ as777 In Ref. the operators $`𝝅`$ and $`𝝆`$ are defined more generically as $`𝝆^mU(x)=(\mathrm{\Gamma }(xr+1)/\mathrm{\Gamma }(xrm+1))U(xm)`$ and $`𝝅U(x)=(xr)(U(x)U(x1))`$, where $`r`$ is a fixed number. Here for simplicity we have set $`r=0`$.

$$𝝆^mU(x)=\frac{\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(xm+1)}U(xm).$$
(34)
This operator has the property

$$𝝆^m𝝆^nU(x)=𝝆^{m+n}U(x);$$
(35)
if the operand is the unity, we omit it and write

$$𝝆^m1=𝝆^m=\frac{\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(xm+1)}.$$
(36)
One defines the operator $`𝝅`$ as

$$𝝅U(x)=x(U(x)U(x1)).$$
(37)
The following properties can be proven:

$`\begin{array}{cc}\hfill [𝝅,𝝆]U(x)& =𝝆U(x),\hfill \\ \hfill p(𝝅)𝝆^mU(x)& =𝝆^mp(𝝅+m)U(x),\hfill \end{array}`$ (38)
$`\begin{array}{cc}\hfill xU(x)& =(𝝅+𝝆)U(x),\hfill \\ \hfill x^nU(x)& ={\displaystyle \underset{k=0}{\overset{n}{}}}\left({\displaystyle \underset{j=k}{\overset{n}{}}}(1)^{jk}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right)S_{jk}𝝅^{nj}\right)𝝆^kU(x),\hfill \\ \hfill p(x)U(x)& =\left[p(𝝅)+p_1(𝝅)𝝆+p_2(𝝅)𝝆^2/2!+\mathrm{}\right]U(x),\hfill \end{array}`$ (39)
where $`p`$ is a polynomial, $`p_n`$ is

$$p_n(\lambda )=\mathrm{\Delta }_np(\lambda )=\underset{j=0}{\overset{n}{}}(1)^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)p(\lambda j),$$

and $`S_{jk}`$ are the Stirling’s numbers of second kind . Using these operators an expansion in factorial series becomes an expansion in powers of $`𝝆^1`$ so that we will able to obtain solutions in factorial series of difference equations in the same manner as power series solutions of differential equations are obtained.

### 4.3 Solution of the homogeneous difference equation

Let us consider the homogeneous difference equation (31) of order $`R`$

$$p_0(x)U^{HO}(x)+p_1(x)U^{HO}(x+1)+\mathrm{}+p_R(x)U^{HO}(x+R)=0;$$
(40)
making the replacement $`xxR`$ and defining $`q_{Ri}(x)=p_i(xR)`$ the equation becomes

$$q_0(x)U^{HO}(x)+q_1(x)U^{HO}(x1)+\mathrm{}+q_R(x)U^{HO}(xR)=0.$$
(41)
Using the previous definitions of operators now we search for a formal solution in factorial series. We make the change of variable $`U^{HO}(x)=\mu ^xV^{HO}(x)`$ in Eq.(41), where $`\mu `$ is an unspecified parameter, obtaining

$$\mu ^Rq_0(x)V^{HO}(x)+\mu ^{R1}q_1(x)V^{HO}(x1)+\mathrm{}+q_R(x)V^{HO}(xR)=0.$$
(42)
Now we multiply the equation by $`x(x1)(x2)\mathrm{}(xR+1)`$ and we observe that $`xV(x1)=𝝆V(x)`$, $`x(x1)V(x2)=𝝆^2V(x)`$, etc.; the equation takes the form

$$\left[\varphi _0(x,\mu )+\varphi _1(x,\mu )𝝆+\mathrm{}+\varphi _R(x,\mu )𝝆^R\right]V^{HO}(x)=0,$$
(43)
where $`\varphi _j`$ are polynomials in $`x`$ and $`\mu `$. The multiplication by $`x`$ is equivalent to the multiplication by $`𝝅+𝝆`$, therefore substituting the relations (39) in Eq.(43) one obtains the *first canonical form* of the difference equation:

$$\left[f_0(𝝅,\mu )+f_1(𝝅,\mu )𝝆+f_2(𝝅,\mu )𝝆^2+\mathrm{}+f_{m+1}(𝝅,\mu )𝝆^{m+1}\right]V^{HO}(x)=0,$$
(44)
where $`f_i`$ are polynomials in $`𝝅`$ and $`\mu `$. In the case of the difference equations encountered in this work we find that $`f_{m+1}(𝝅,\mu )`$ is independent of $`𝝅`$, but not independent of $`\mu `$; therefore $`f_{m+1}(𝝅,\mu )=f_{m+1}(\mu )`$. The algebraic equation in $`\mu `$

$$f_{m+1}(\mu )=0$$
(45)
is the *characteristic equation*888 If the coefficients of the equation are $`p_i(x)=_{j=0}^gp_{ij}x^{gj}`$, the characteristic equation has the explicit expression $`_{i=0}^Rp_{i0}\mu ^i=0`$.. It turns out that the characteristic equation has always $`R`$ solutions different from zero. Let $`\mu _1`$, $`\mu _2`$, …, $`\mu _\lambda `$ be the $`\lambda `$ *distinct* solutions of Eq.(45). For each of these values, $`\mu =\mu _i`$, $`i=1,\mathrm{},\lambda `$, the first canonical form (44) takes the form

$$\left[f_0(𝝅)+f_1(𝝅)𝝆+f_2(𝝅)𝝆^2+\mathrm{}+f_m(𝝅)𝝆^m\right]V^{HO}(x)=0.$$
(46)
Now we try to satisfy the equation (46) in $`V^{HO}`$ with the factorial series

$$V^{HO}(x)=\underset{s=0}{\overset{\mathrm{}}{}}\frac{a_s\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(xK+s+1)}=\underset{s=0}{\overset{\mathrm{}}{}}a_s𝝆^{Ks}=a_0𝝆^K+a_1𝝆^{K1}+\mathrm{},$$
(47)
whose asymptotic behaviour for large $`x`$ is $`V^{HO}(x)a_0x^K`$. Substituting Eq.(47) in Eq.(46) one obtains the recurrence relations

$`\begin{array}{cc}& a_0f_m(K+m)=0,\hfill \\ & a_1f_m(K+m1)+a_0f_{m1}(K+m1)=0,\hfill \end{array}`$
$`\begin{array}{cc}& \mathrm{}\hfill \\ & a_sf_m(K+ms)+a_{s1}f_{m1}(K+ms)+\mathrm{}+a_{sm}f_0(K+ms)=0(sm).\hfill \end{array}`$ (48)
Supposing that $`a_00`$, the equation

$$f_m(K+m)=0$$
(49)
is the *indicial equation*. Let $`K_1,K_2,\mathrm{},K_\nu `$ the roots of this equation. In the case of the difference equations encountered in this work all the roots turn out to be distinct, and $`\nu `$ turns out to be the multiplicity of $`\mu _i`$. For each of these values of $`K`$ we solve the system of recurrence relations. If there are no roots differing by a positive integer, $`f_m(K+ms)0`$ for $`s1`$, and therefore the recurrence relation (LABEL:systa2) provides the coefficient $`a_s`$ of the factorial series for every $`s`$. If there are roots differing by a positive integer (so-called *congruent* roots) we find $`f_m(K+ms_0)=0`$ for some $`s_0`$, so that the term $`a_{s_0}f_m(K+ms_0)`$ vanishes from the relation (LABEL:systa2) used to obtain $`a_{s_0}`$. The remaining terms of this relation always vanish999 Otherwise $`V^{HO}(x)`$ would be a derivative with respect to $`K`$ of a factorial series, with asymptotic behaviour $`x^K\mathrm{log}^n(x)`$.  in the case of the difference equations encountered in this work, therefore the value of $`a_{s_0}`$ remains undetermined and can be chosen at will (usually one puts $`a_{s_0}=0`$). For details on the convergence of the factorial series so obtained, see section 6.1.

In order to obtain the general solution of the difference equation, for each distinct solution of the characteristic equation $`\mu _i`$, $`i=1,\mathrm{},\lambda `$, we find an indicial equation, whose solutions are $`K_{ij}`$, $`j=1,\mathrm{},\nu _i`$. For each solution of the indicial equation we solve the system of recurrence relations between the coefficients $`a_s^{(i,j)}`$ and we find $`\nu _i`$ solutions of Eq.(46),

$$V_{ij}^{HO}(x)=\underset{s=0}{\overset{\mathrm{}}{}}\frac{a_s^{(i,j)}\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(xK_{ij}+s+1)},j=1,\mathrm{},\nu _i.$$
(50)
Multiplying by $`\mu _i^x`$ we find solutions of the difference equation (40). In total we find $`_{i=1}^\lambda \nu _i=R`$ different solutions. The general solution of the homogeneous difference equation (40) will be a linear combination of the these solutions with periodic functions $`\stackrel{~}{\omega }_{ij}(x)`$ of period one as coefficients,

$$U^{HO}(x)=\underset{i=1}{\overset{\lambda }{}}\underset{j=1}{\overset{\nu _i}{}}\stackrel{~}{\omega }_{ij}(x)\mu _i^xV_{ij}^{HO}(x).$$
(51)
### 4.4 Solution of the nonhomogeneous difference equation

Let us consider the nonhomogeneous difference equation (29) of order $`R`$

$$p_0(x)U^{NH}(x)+p_1(x)U^{NH}(x+1)+\mathrm{}+p_R(x)U^{NH}(x+R)=F(x);$$
(52)
making the same replacements as Eq.(40) the equation becomes

$$q_0(x)U^{NH}(x)+q_1(x)U^{NH}(x1)+\mathrm{}+q_R(x)U^{NH}(xR)=F^{}(x),$$
(53)
where $`F^{}(x)=F(xR)`$. In the difference equations encountered in this work, the right-hand side can be written as a sum of factorial series expansions. We assume for simplicity that $`F^{}(x)=\mu ^xT(x)`$ and $`T(x)`$ has the factorial series expansion

$$T(x)=c_0𝝆^K+c_1𝝆^{K1}+c_2𝝆^{K2}+\mathrm{},$$
(54)
where $`\mu `$, $`K`$ and $`c_i`$ are known.

Making the substitution $`U^{NH}(x)=\mu ^xV^{NH}(x)`$ in Eq.(53) and using the relations (38-39) we obtain the canonical form

$$\left[f_0(𝝅)+f_1(𝝅)𝝆+f_2(𝝅)𝝆^2+\mathrm{}+f_m(𝝅)𝝆^m\right]V^{NH}(x)=T(x).$$
(55)
The expansion in factorial series of $`V^{NH}(x)`$ has the form

$$V^{NH}(x)=a_0𝝆^{Km}+a_1𝝆^{Km1}+a_2𝝆^{Km2}+\mathrm{};$$
(56)
inserting this expansion in Eq.(55) and equating the coefficients of the powers of $`𝝆`$ one finds the recurrence relations

$`\begin{array}{cc}& a_0f_m(K)=c_0,\hfill \\ & a_1f_m(K1)+a_0f_{m1}(K1)=c_1,\hfill \end{array}`$
$`\begin{array}{cc}& \mathrm{}\hfill \\ & a_sf_m(Ks)+a_{s1}f_{m1}(Ks)+\mathrm{}+a_{sm}f_0(Ks)=c_s(sm).\hfill \end{array}`$ (57)
Solving this system we can find the coefficients $`a_s`$. If $`f_m(Ks_0)=0`$ for some $`s_0`$, the term $`a_{s_0}f_m(Ks_0)`$ vanishes from the relation (LABEL:systb2) and, in the difference equations encountered in this work, the remaining terms form the trivial identity $`c_s=c_s`$, so that the value of $`a_{s_0}`$ remains undetermined and can be chosen at will. If the factorial series (56) converges, $`U^{NH}(x)=\mu ^xV^{NH}(x)`$ is a particular solution of the nonhomogeneous difference equation.

The case $`F^{}(x)=\mu ^xp(x)T(x)`$, where $`p`$ is a polynomial, can be reduced to the previous case by transforming $`x`$ into $`𝝅+𝝆`$ in the polynomial and letting the operators act on the expansion of $`T(x)`$ following Eq.(38).

## 5 Determination of arbitrary constants

The general solution of the difference equation (29) can be written as

$$U(x)=\underset{j=1}{\overset{R}{}}\stackrel{~}{\omega }_j(x)U_j^{HO}(x)+U^{NH}(x),$$
(58)
where $`\stackrel{~}{\omega }_j(x)`$ are periodic functions of period 1. If we consider only integer values of $`x`$, the values of $`\stackrel{~}{\omega }_j(x)`$ are independent of $`x`$ so that we can replace them with arbitrary constants $`\eta _j`$:

$$U(x)=\underset{j=1}{\overset{R}{}}\eta _jU_j^{HO}(x)+U^{NH}(x).$$
(59)
The unambiguous solution of the difference equation requires the determination of the $`R`$ constants $`\eta _j`$. The value of these constants can be determined

1. from the large-$`x`$ behaviour of the solution (59) and of the integral (28), by equating the first coefficients of the expansions in factorial series;
2. by equating the values at $`x=0`$ of Eq.(59) and Eq.(28).

As we will see in the following sections, the method 1 may provide the values of all the constants, but the determination of the coefficients of the factorial series turns out to be a simple task only for euclidean massive integrals (see sections 5.1, 5.2 and 5.3), where it involves integrals with one-loop less, and gets more complicated for non-euclidean integrals or integrals with zero masses (see section 5.4); on the contrary the method 2 has the limitation that it provides only one relation between the constants, but on the other hand it is not affected by the kind of the external momenta or masses (see section 5.5).

### 5.1 Large-$`x`$ behaviour of integrals: euclidean massive case

In this section we work out a relation between the coefficients of the expansion in factorial series of the master integral (28) and simpler integrals with one-loop less. We assume that all the external momenta are euclidean, that is, the matrix of the scalar products $`p_ip_j`$ is semidefinite non-negative, and no mass is zero.

Choosing the momentum routing of the integral (28) such that $`D_1=k_1^2+m_1^2`$, we write it as

$$U(x)=\frac{[d^Dk_1]}{(k_1^2+m_1^2)^x}g(k_1),$$
(60)
where $`g`$ includes the contribution of the remaining $`N_k1`$ loops

$$g(k_1)=[d^Dk_2]\mathrm{}[d^Dk_{N_k}]\frac{_{j=1}^{N_{sp}N_d}(pk\text{ irred.})_j^{\beta _j}}{D_2D_3\mathrm{}D_{N_d}}.$$
(61)
Note that the function $`g`$ also depends on the external momenta $`p_j`$ and the masses $`m_2`$, $`m_3`$, $`\mathrm{},`$ $`m_{N_d}`$ of the other denominators. From a graphical point of view, $`g(k_1)`$ corresponds to the original diagram with the line $`D_1`$ cut. Introducing hyperspherical polar coordinates for the integration over $`k_1`$

$$d^Dk_1=|k_1|^{D1}d|k_1|d\mathrm{\Omega }_D(\widehat{k}_1),$$
(62)
($`\mathrm{\Omega }_D=2\pi ^{D/2}/\mathrm{\Gamma }(D/2)`$ is the $`D`$-dimensional solid angle) and separating the angular and radial part, Eq.(60) becomes

$$U(x)=\frac{1}{\mathrm{\Gamma }(D/2)}_0^{\mathrm{}}\frac{dk_1^2(k_1^2)^{D/21}}{(k_1^2+m_1^2)^x}f(k_1^2),$$
(63)
where $`f`$ is the angular mean over $`k_1`$ of $`g`$

$$f(k_1^2)=\frac{1}{\mathrm{\Omega }_D}𝑑\mathrm{\Omega }_D(\widehat{k}_1)g(k_1).$$
(64)
For large $`x`$ the factor $`(k_1^2+m_1^2)^x`$ of Eq.(63) peaks strongly around $`k_1^2=0`$; because of the assumption on external momenta and masses, the function $`f(k_1^2)`$ has no singularities for $`k_1^20`$, and therefore we expect that the large-$`x`$ behaviour of the integral is controlled only by the behaviour of $`f(k_1^2)`$ near $`k_1^2=0`$. Making the change of variable $`k_1^2=m_1^2{\displaystyle \frac{u}{1u}}`$ in Eq.(63) we obtain

$$U(x)=\frac{(m_1^2)^{D/2x}}{\mathrm{\Gamma }(D/2)}_0^1𝑑uu^{D/21}(1u)^{x1D/2}\stackrel{~}{f}(u),$$
(65)
where $`\stackrel{~}{f}(u)=f(m_1^2u/(1u))`$. We expand $`\stackrel{~}{f}(u)`$ in $`u`$

$$\stackrel{~}{f}(u)=u^\alpha (1u)^\beta \underset{s=0}{\overset{\mathrm{}}{}}b_su^s;$$
(66)
$`\alpha `$ is an integer greater than or equal to zero ($`f`$ is regular in the origin), while the factor $`(1u)^\beta `$ has been introduced for convenience. Expressing the integrals over $`u`$ in terms of Beta function (which takes care of analytical continuation if the integral diverges at the endpoints), and choosing $`\beta =D/2+1`$, Eq.(65) takes the desired form of a factorial series

$$U(x)=\mu _0^x\underset{s=0}{\overset{\mathrm{}}{}}a_s𝝆^{K_0s}=\mu _0^x\underset{s=0}{\overset{\mathrm{}}{}}a_s\frac{\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(x+1K_0+s)},$$
(67)
where

$$\mu _0=1/m_1^2,K_0=D/2\alpha ,\text{and }a_s=b_sm_1^D\mathrm{\Gamma }(s+D/2+\alpha )/\mathrm{\Gamma }(D/2).$$
(68)
The coefficients $`b_s`$ of Eq.(66) can be easily expressed in terms of the coefficients $`f_s`$ of the expansion in $`k_1^2`$ of $`f(k_1^2)`$

$$f(k_1^2)=(k_1^2)^\alpha \underset{s=0}{\overset{\mathrm{}}{}}f_sk_1^{2s},$$
(69)
($`b_0=m_1^{2\alpha }f_0`$, $`\mathrm{}`$), so that the large-$`x`$ behaviour of $`U(x)`$ proves to be determined by the behaviour of $`f(k_1^2)`$ for small $`k_1^2`$. For large $`x`$ the leading behaviour is given by the first term of Eq.(67)

$$U(x)(m_1^2)^{D/2x+\alpha }x^{D/2\alpha }f_0\frac{\mathrm{\Gamma }(D/2+\alpha )}{\mathrm{\Gamma }(D/2)}.$$
(70)
In the frequent case of an integral with numerator equal to one, $`\alpha =0`$ and $`f_0=f(0)`$, it becomes

$$U(x)(m_1^2)^{D/2x}x^{D/2}f(0).$$
(71)
Now we equate Eq.(67) with the expansion in factorial series of the general solution (59) (see Eq.(51) and Eq.(56)), obtaining

$$\left(\frac{1}{m_1^2}\right)^x\underset{s=0}{\overset{\mathrm{}}{}}a_s𝝆^{D/2\alpha s}=\underset{j=1}{\overset{R}{}}\eta _j\mu _j^x\underset{s=0}{\overset{\mathrm{}}{}}\widehat{a}_{js}𝝆^{K_js}+\underset{l}{}(\mu _l^{NH})^x\underset{s=0}{\overset{\mathrm{}}{}}a_{ls}^{NH}𝝆^{K_l^{NH}s}.$$
(72)
In this equation only the $`R`$ constants $`\eta _j`$ are unknowns, each one corresponding to a solution $`U_j^{HO}`$ with a different pair of values $`\mu _j`$ and $`K_j`$, $`\mu _j`$ being one solution of the characteristic equation of the homogeneous equation and $`K_j`$ being the solution of the corresponding indicial equation; we have set $`\widehat{a}_{j0}=1`$. We have also assumed that the expansion of the nonhomogeneous term contains a sum of expansions with different pairs of values $`\mu _l^{NH}`$ and $`K_l^{NH}`$.

Fortunately, the number of constants $`\eta _j`$ to find can be drastically reduced. If the difference equation is homogeneous, the constants $`\eta _j`$ which may be different from zero are only those such that the corresponding $`\mu _j`$ and $`K_j`$ satisfy the condition

$$\mu _j=1/m_1^2,K_j+D/2+\alpha =\text{integer}0\eta _j0;$$
(73)
all the other constants must be zero.

In the case of nonhomogeneous equation, one must recall that it is part of a triangular system of difference equations; the nonhomogeneous term receives contributions from master integrals with a smaller number of denominators, which in their turn satisfy other nonhomogeneous difference equations, up to the simplest master integrals which satisfy homogeneous equations. Clearly $`\mu _l^{NH}=1/m_1^2`$ and $`K_l^{NH}+D/2`$ is always an integer, therefore $`\eta _j`$ may be different from zero if $`\mu _j`$ and $`K_j`$ satisfy Eq.(73), or the new condition

$$\mu _j=1/m_1^2,0<K_j+D/2+\alpha =\text{integer}\underset{l}{\mathrm{max}}K_l^{NH}+D/2+\alpha \eta _j0$$
(74)
which implies cancellations between the first coefficients of the homogeneous and nonhomogeneous expansions. We found that in the equations considered in this work, the condition (74) is never satisfied, and the condition (73) is satisfied by a small number of pairs of $`\mu _j`$ and $`K_j`$, often only one; if no pair satisfies the conditions, $`U(x)`$ is completely determined by the nonhomogeneous term.

The values of the non-zero constants $`\eta _j`$ are determined by comparing the coefficients of the same powers of $`𝝆`$ of the two sides of Eq.(72). The required few coefficients $`a_s`$ are easily calculated from the first coefficients $`f_s`$ of the expansion in $`k_1^2`$ of $`f(k_1^2)`$.

The required coefficients $`f_s`$ are calculated by expanding numerators and denominators of Eq.(61) for small $`k_1`$, and by performing the angular integration over $`k_1`$ of Eq.(64). Angular integrals are straightforward as they contain exclusively powers of scalar products containing $`k_1`$. As a result, the coefficients $`f_s`$ will be expressed by integrals over $`k_2`$,…$`k_{N_k}`$, belonging to diagrams with *one loop less* (for example, if the numerator of Eq.(61) is unity then $`f_0=f(0)=g(0)`$). These integrals can be expressed using algorithm 1 as combinations of the master integrals of the new diagrams; in their turn, these new master integrals can be calculated by inserting an exponent in a denominator and building and solving new difference equations. The arbitrary constants of their solutions can be found using the large-exponent behaviours, which can be expressed in terms of integrals with another loop less, and so on. In this way we can explicitly calculate the values of all the arbitrary constants $`\eta _j`$. This fact is very important.

### 5.2 Non-euclidean case: deformation of the radial path

Now we consider the non-euclidean case: the matrix of the scalar products $`p_ip_j`$ is definite negative. We write Eq.(63) as

$$U(x,P)=\frac{1}{\mathrm{\Gamma }(D/2)}_{l_0}\frac{dk_1^2(k_1^2)^{D/21}}{(k_1^2+m_1^2)^x}f(k_1^2,P),$$
(75)
where we explicitly show the dependence on the external momenta $`P=`$ $`\{p_1,`$ $`\mathrm{},p_{N_p}\}`$, and where $`l_0`$ indicates the path of the radial integration. For the moment the path is assumed to be the positive axis in the complex plane.

The integral $`U(x,P)`$, considered as a function of the external momenta $`P`$, is defined in the non-euclidean region by an analytical continuation through a generic path in the complex $`P`$-space. The path begins in some euclidean initial point $`P_{in}`$ and ends in the desired non-euclidean final point $`P_{end}`$. It is possible that in some intermediate point of the continuation path in the $`P`$-space, for a particular value of the external momenta $`P_s`$, $`f(k_1^2,P_s)`$ is singular in a point $`k_1^2=q^2>0`$, just on the radial integration path, breaking off the analytical continuation. If this singularity in $`k_1^2`$ cannot be avoided by modifying the continuation path in the $`P`$-space, then we are forced to deform101010 The deformation of the radial integration path was first discussed in  for self-energy diagrams in $`D=4`$.  the integration path in the $`k_1^2`$-space by turning around the singularity; the singularity in $`k_1^2`$ itself moves in the complex $`k_1^2`$-plane to a final point $`k_{end}^2`$ when we complete the analytical continuation from $`P_s`$ to $`P_{end}`$. The final integration path $`l_0`$ starts in $`k_1^2=0`$, turns around the singularity $`k_{end}^2`$ (usually on the negative axis) and comes back to $`k_1^2+\mathrm{}`$. In the general case of multiple singularities the path turns around all them. The deformation of the radial path must be performed when the values of the scalar products $`p_ip_j`$ get over some “deformation thresholds” in the non-euclidean region, the exact point of these thresholds being depending on the values of masses and the structure of the diagram.

One can show using Feynman parameters that the deformation thresholds are the thresholds (anomalous thresholds included) in $`P`$ of $`g(0,P)`$, the diagram obtained by eliminating the line $`D_1`$ and setting $`k_1=0`$ everywhere. Possible denominators depending only on $`k_1`$ give “additional” deformation thresholds of the kind $`p^2=m^2`$.

Once the deformation thresholds are determined, we can split the non-euclidean region of the $`P`$-space into two regions:

1. one region below the deformation thresholds which adjoins the euclidean region, where the integration path is the same as the euclidean region;
2. the remaining region above the deformation thresholds where the integration path turns around some singularities.

#### 5.2.1 Example: deformation of the path for the one-loop self-energy integral

As an example, now we consider the simplest case, the one-loop integral

$$I=\frac{[d^Dk]}{(k^2+m_1^2)((pk)^2+m_2^2)}=\frac{1}{\mathrm{\Gamma }(D/2)}_{l_0}\frac{dk^2(k^2)^{D/21}}{k^2+m_1^2}\frac{1}{\mathrm{\Omega }_D}\frac{d\mathrm{\Omega }_D(\widehat{k})}{(pk)^2+m_2^2}.$$
(76)
The angular integral is111111 Here we generalize to arbitrary $`D`$ dimensions the result obtained in  in the four-dimensional case.

$$\frac{1}{\mathrm{\Omega }_D}\frac{d\mathrm{\Omega }_D(\widehat{k})}{(pk)^2+m_2^2}=\frac{Z}{pk}\mathrm{\Phi }(Z^2),$$
(77)
where $`\mathrm{\Phi }(Z^2)`$ is the hypergeometric function $`F(1,2D/2;D/2;Z^2)`$,

$$Z^{\pm 1}=\frac{p^2+k^2+m_2^2R(p^2,k^2,m_2^2)}{2pk},p=\sqrt{p^2},k=\sqrt{k^2},$$
(78)
and $`R(x,y,z)`$ is the usual two-body phase space square root

$$R(x,y,z)=\sqrt{x^2+y^2+z^22xy2xz2yz}.$$
(79)
The angular integral (77) has two branching points in the complex $`k^2`$-plane, $`k^2=(p\pm im_2)^2`$. Considering Eq.(76) in the euclidean case, $`p^2>0`$, the path of integration is always the positive real axis. If $`p^2<0`$, Eq.(76) must be analytically continued in the complex $`p`$plane from a generic initial euclidean point $`p=p_{in}>0`$ to the non-euclidean final point $`p=i\sqrt{p^2}`$ which is on the imaginary axis. The initial point has $`\mathrm{}p_{in}=0`$. If $`p^2m_2^2`$ no radial singularity appears so that

$$I=\frac{1}{\mathrm{\Gamma }(D/2)}_0^{\mathrm{}}\frac{dk^2(k^2)^{D/21}}{k^2+m_1^2}\frac{Z}{pk}\mathrm{\Phi }(Z^2).$$
(80)
If $`p^2<m_2^2`$, the final point has $`\mathrm{}p>m_2`$ so that any path connecting the initial and final point must cross the line $`\mathrm{}p=m_2`$ in some point. Let $`p_c=im_2+c`$ be this point; if we consider Eq.(77) with momentum $`p_c`$, one of the branching points falls in $`k^2=c^2>0`$ which is exactly on the path of integration. Therefore the path of integration $`l_0`$ must be deformed avoiding the branching point $`(pim_2)^2`$ which is on the negative real axis. One chooses as path of integration first the segment $`((pim_2)^2,0)`$ above the negative real axis, taking the square root $`R`$ with the plus sign and using $`Z^1`$ in the place of $`Z`$, then the straight line $`((pim_2)^2,\mathrm{})`$ below the axis. Therefore, if $`p^2<m_2^2`$ Eq.(76) becomes

$$I=\frac{1}{\mathrm{\Gamma }(D/2)}_0^{(pim_2)^2}\frac{dk^2(k^2)^{D/21}}{k^2+m_1^2}\frac{Z^1}{pk}\mathrm{\Phi }(Z^2)+\frac{1}{\mathrm{\Gamma }(D/2)}_{(pim_2)^2}^{\mathrm{}}\frac{dk^2(k^2)^{D/21}}{k^2+m_1^2}\frac{Z}{pk}\mathrm{\Phi }(Z^2).$$
(81)
This result will be used in section 7.2.

### 5.3 Large-$`x`$ behaviour of integrals: non-euclidean massive case below and at the deformation threshold

In the region below the deformation thresholds, in the case of non-euclidean massive integrals, the determination of the large-$`x`$ behaviour is quite similar to that described in section 5.1.

At the deformation threshold the determination of the large-$`x`$ behaviour presents some peculiarities, since some external momenta have on-mass-shell values. In the integral (61) denominators of the form $`k_1^22\overline{p}k_1`$ appear, where $`\overline{p}`$ is some on-mass-shell momentum. Considering the expansion in $`k_1`$ of $`g(k_1)`$, these denominators vanish at $`k_1=0`$ and cannot be expanded in a way as straightforward as for off-mass-shell denominators; they must be included in the angular integrals, complicating quite a lot the integration.

In order to explain the situation, let us consider here the case of a one-loop diagram with $`N+1`$ external lines, and the integral

$$\frac{[d^Dk]}{(k^2+m_0^2)((p_1k)^2+m_1^2)\mathrm{}((p_Nk)^2+m_N^2)}.$$
(82)
Setting $`k=0`$ we see that the deformation threshold is $`p_i^2=m_i^2`$, $`i=1`$, …, $`N`$. Now we set the external momenta to these on-mass-shell values, and we consider the integral

$$W(x)=\frac{[d^Dk]}{(k^2+m_0^2)^x(k^22p_1k)\mathrm{}(k^22p_Nk)},$$
(83)
whose we want to calculate the large-$`x`$ leading behaviour. Following the notation of section 5.1 we write

$$W(x)=_0^{\mathrm{}}\frac{dk^2(k^2)^{D/21}}{(k^2+m_0^2)^x}\frac{f(k^2)}{\mathrm{\Gamma }(D/2)},$$
(84)
$$f(k^2)=\frac{1}{\mathrm{\Omega }_D}\frac{d\mathrm{\Omega }_D(\widehat{k})}{(k^22p_1k)\mathrm{}(k^22p_Nk)}.$$
(85)
Now we extract the leading behaviour of $`f(k^2)`$ for $`k^20`$. Inserting Feynman parameters $`x_i`$, $`i=1`$, $`\mathrm{}`$, $`N1`$ one finds

$$f(k^2)=\frac{\mathrm{\Gamma }(N)}{\mathrm{\Omega }_D}𝑑x_1\mathrm{}𝑑x_{N1}\frac{d\mathrm{\Omega }_D(\widehat{k})}{(k^22P(x_i)k)^N},$$
(86)
where $`P(x_i)=_{i=1}^Nx_ip_i`$ and $`_{i=1}^Nx_i=1`$. The angular integral in $`D`$ dimensions can be expressed using a generalization of the formula (77)

$$\frac{1}{\mathrm{\Omega }_D}\frac{d\mathrm{\Omega }_D(\widehat{k})}{((pk)^2+m^2)^x}=\left(\frac{Z}{pk}\right)^xF(x,x+1D/2;D/2;Z^2);$$
(87)
using properties of hypergeometric function one finds for $`k^20`$

$$\frac{1}{\mathrm{\Omega }_D}\frac{d\mathrm{\Omega }_D(\widehat{k})}{(k^22P(x_i)k)^N}(P^2(x_i)k^2)^{N/2}\frac{\mathrm{\Gamma }(D/2)\mathrm{\Gamma }(N/2)}{2\mathrm{\Gamma }(N)\mathrm{\Gamma }((DN)/2)},$$
(88)
then

$$f(k^2)\frac{\mathrm{\Gamma }(D/2)\mathrm{\Gamma }(N/2)}{2\mathrm{\Gamma }((DN)/2)}(k^2)^{N/2}\frac{dx_1\mathrm{}dx_{N1}}{\left(P^2(x_i)\right)^{N/2}}.$$
(89)
We see that $`f(k^2)`$ is singular in $`k=0`$ and, if $`N>D`$, the integral over $`k^2`$ of Eq.(84) diverges for small $`k^2`$; Eq.(70) with arbitrary $`\alpha `$ takes care of the necessary analytical continuation.

Let us now consider a one-loop integral in $`E`$ dimensions with the same denominators as the angular integral (85)

$$L_E(p_i)\frac{1}{2}\frac{[d^Eq]}{(q^22p_1q)\mathrm{}(q^22p_Nq)};$$
(90)
introducing Feynman parameters in the same way as for Eq.(86) we obtain

$$L_E(p_i)=\frac{\mathrm{\Gamma }(NE/2)}{2}\frac{dx_1\mathrm{}dx_{N1}}{\left(P^2(x_i)\right)^{NE/2}}.$$
(91)
If we choose $`E=N`$ the integrals over $`x_i`$ of Eq.(89) and Eq.(91) become identical, so that we can rewrite Eq.(89) as

$$f(k^2)(k^2)^{N/2}\frac{\mathrm{\Gamma }(D/2)}{\mathrm{\Gamma }((DN)/2)}L_N(p_i).$$
(92)
This result shows that the angular integral (85) with $`N`$ on-mass-shell denominators in the limit $`k^20`$ is proportional to a one-loop integral with the same denominators (one less than Eq.(83)), calculated in a *number of dimensions equal to the number of on-mass-shell denominators*. Moreover, using Eq.(92) and Eq.(70) (with $`\alpha =N/2`$ according to Eq.(89)) one finds the large-$`x`$ leading behaviour of Eq.(83):

$$W(x)(m_0^2)^{D/2N/2x}x^{D/2+N/2}L_N(p_i).$$
(93)
This result is similar to euclidean leading behaviour for a one-loop vacuum integral (Eq.(71) with $`f(0)=1`$); the differences are the change of the exponents and the multiplication by the factor $`L_N(p_i)`$ which is *independent* of the number of dimensions $`D`$ of the integral (83). As a consequence of the divergence of the radial integral and of the necessary analytical continuation, we observe an apparently paradoxical result: if $`N>D`$ and $`m_0=1`$, from Eq.(93) we see that the integral $`W(x)`$ increases as the exponent $`x`$ increases, while in the euclidean case the integral $`U(x)`$ decreases as $`x`$ increases. The integral $`L_N(p_i)`$ may be calculated analytically, for example, extracting the $`k^20`$ behaviour from the known analytical expressions of the angular integral (85) for $`D=4`$ and $`N=1`$  and $`N=2,3`$ , or numerically, by inserting an exponent in one denominator of Eq.(90) and building and solving a difference equation, or by using the identities of the section 5.5 (for example see Eq.(152)). In section 9.2 we will need the following values $`\stackrel{~}{L}_N`$ of $`L_N(p_i)`$ in the on-mass-shell and equal masses case $`m_i=1`$, $`p_i^2=1`$ and $`(p_ip_j)^2=1`$ for every $`i`$ and $`j`$:

$$\begin{array}{cc}\hfill \stackrel{~}{L}_1=& \sqrt{\pi }/2,\stackrel{~}{L}_2=\pi \sqrt{3}/9,\hfill \\ \hfill \stackrel{~}{L}_3=& \left(\mathrm{arctan}\sqrt{8}\mathrm{arctan}\sqrt{3}\right)\sqrt{9\pi /8},\stackrel{~}{L}_4=0.172751462\mathrm{}.\hfill \end{array}$$
(94)
The generalization of Eq.(93) to multi-loop integrals is straightforward; considering

$$W^{}(x)=\frac{[d^Dk]}{(k^2+m_0^2)^x(k^22p_1k)\mathrm{}(k^22p_Nk)}h(k),$$
(95)
$$h(k)=\frac{[d^Dk_2]\mathrm{}[d^Dk_{N_k}]}{D_{N+2}\mathrm{}D_{N_d}},$$
(96)
the large-$`x`$ leading behaviour is

$$W^{}(x)(m_0^2)^{D/2N/2x}x^{D/2+N/2}L_N(p_i)h(0).$$
(97)
The condition (73) in this case must be modified to

$$\mu _j=1/m_1^2,K_j+D/2N/2=\text{integer or half-integer}0\eta _j0.$$
(98)
### 5.4  Non-euclidean case above the deformation threshold and massless case

The cases described above, euclidean and non-euclidean below the deformation thresholds, are both united by the fact that the large-$`x`$ behaviour of $`U(x)`$ depends only on the behaviour of $`f(k_1^2)`$ near the origin. In these cases it is easy to express this behaviour in terms of simpler diagrams obtained by putting $`k_1=0`$ in the main diagram. The situation gets more complicated if $`f(k_1^2)`$ has other singularities on the integration path: for the $`i`$th singularity placed in $`q_i^20`$, an additional contribution to the large-$`x`$ behaviour of $`U(x)`$ appears, of the form $`(q_i^2+m_1^2)^xV(x)`$, with $`V(x)x^K`$ for large $`x`$. The precise form of $`V(x)`$ depends on the behaviour of $`f`$ near the singularity; unfortunately its determination requires the calculation of angular integrals for $`k_1^20`$ (see Eq.(64))

$$_{k_1^2q_i^2}𝑑\mathrm{\Omega }_D(\widehat{k}_1)g(k_1),$$
(99)
which is difficult and case-dependent. In this work we have performed these calculations only in the one-loop example of section 7.2. In the case of more complicated diagrams we have preferred to avoid the calculation of these angular integrals and to look at the problem from a different point of view, using a method with broader applicability based on differential equations (see section 10).

The cases where singularities of $`f`$ other than the origin appear are:

* The non-euclidean massive case, above the deformation threshold: the path of integration must be deformed, turning around a number of singularities of $`f(k_1^2)`$. See section 7.2.1 for an example.
* The zero mass case, where $`m_1>0`$ and some of the masses $`m_2`$, …, $`m_{N_d}`$ are zero. If the external momenta are euclidean, singularities on the positive axis of $`k_1^2`$ may appear, but this fact does *not* require deformation of the path. See section 7.2.2 for an example.

The situation worsens if the mass $`m_1`$ of the denominator raised to $`x`$ is zero. In this case the behaviour of $`U(x)`$ for large $`x`$ depends on the behaviour of $`f(k_1^2)`$ on the whole real axis, and not in some isolated points. This can be easily understood by noting that in this case $`U(x+D/21)`$ is the Mellin transform of $`f(k_1^2)`$. See section 7.2.3 for an example.

### 5.5 Zero exponent condition

A very useful relation between the constants $`\eta _j`$ arises from Eq.(59) calculated at $`x=0`$:

$$U(0)U^{NH}(0)=\underset{j=1}{\overset{R}{}}\eta _jU_j^{HO}(0).$$
(100)
$`U(0)`$ is an integral without the denominator $`D_1`$; if the value of $`U(0)`$ can be found by some method, as an example, by adding an exponent to $`D_2`$ and by solving a difference equation, and if at least one $`U_j^{HO}(0)`$ has a non-zero value, Eq.(100) establishes one new relation between the constants $`\eta _j`$. The existence of this relation must be verified in each particular case. The advantage of this relation over the relations found using the asymptotic behaviour of $`U(x)`$, is that it is valid for every kind of integral, regardless of the value of masses or external momenta.

It is possible to construct identities analogous to Eq.(100), by inserting in Eq.(59) $`x=1,2,\mathrm{}`$ instead of $`x=0`$. Unfortunately in all the analyzed cases the identities so found turned out to be equivalent to Eq.(100).

## 6 Evaluation of factorial series

Once the constants $`\eta _j`$ of Eq.(59) have been determined using the methods described in section 5, in order to obtain the value of the master integral $`U(1)`$ we must calculate the values of the homogeneous solutions $`U_j^{HO}(1)`$ and of the nonhomogeneous solution $`U^{NH}(1)`$.

### 6.1 Convergence of factorial series and instabilities of recurrence relations

Let us unify the notation by considering the solution $`U^{(\alpha )}(x)`$, where $`(\alpha )`$ indicates one of the solutions of the homogeneous or nonhomogeneous equation, and expand it in factorial series:

$$U^{(\alpha )}(x)=\left(\mu ^{(\alpha )}\right)^x\underset{s=0}{\overset{\mathrm{}}{}}a_s^{(\alpha )}\frac{\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(x+1K^{(\alpha )}+s)}.$$
(101)
Convergence of the series depends on the value of the abscissa of convergence $`\lambda `$. Analyzing the large-$`s`$ behaviour of the coefficients $`a_s^{(\alpha )}`$ one finds that the series has $`\lambda <\mathrm{}`$ if none of the solutions $`\mu _j`$ of the characteristic equation (45) satisfies the condition

$$0<|\mu _j/\mu ^{(\alpha )}1|<1,j=1,\mathrm{},R.$$
(102)
If $`\lambda =\mathrm{}`$ the series is everywhere divergent and the expansion in factorial series is only formal; in this case another method must be used to calculate $`U^{(\alpha )}(1)`$ (see section 8 on the Laplace’s transformation).

If $`\lambda <\mathrm{}`$ the convergence of the factorial series is logarithmic, that is, if $`S_m(x)`$ is the sum of the first $`m`$ terms of Eq.(101), one finds that $`|S_m(x)U^{(\alpha )}(x)|m^{\lambda x}`$ for large $`m`$. The series converges if $`x>\lambda `$, more and more quickly as $`x`$ increases. As usually one finds $`\lambda 1`$, it is not convenient or possible to calculate $`U^{(\alpha )}(1)`$ directly by summing the factorial series. Therefore, chosen a number $`x_{max}`$ conveniently large, one calculates $`U^{(\alpha )}(x)`$ for $`R`$ contiguous values of $`x`$, $`U^{(\alpha )}(x_{max})`$, $`U^{(\alpha )}(x_{max}+1)`$, …, $`U^{(\alpha )}(x_{max}+R1)`$, where the series converges faster, and one uses repeatedly the corresponding recurrence relation, Eq.(29) or Eq.(31), in order to obtain the values of $`U^{(\alpha )}(x)`$ for $`x=x_{max}1,x_{max}2,\mathrm{}`$ up to $`x=1`$. A drawback of this procedure is that the recurrence relation may be unstable, so that each iteration causes a loss of precision.

Let $`A^1=\mathrm{min}_j|\mu _j/\mu ^{(\alpha )}|`$. The recurrence relation is unstable if

1. $`A>1`$: in this case each iteration increases the error on $`U^{(\alpha )}(x)`$ of a factor $`A`$.
2. $`A=1`$, and $`\mu ^{(\alpha )}`$ is a root of Eq.(45) of multiplicity $`m>1`$: in this case $`n`$ iterations of the recurrence relation increase the error on $`U^{(\alpha )}(x)`$ of a factor $`n^{m1}`$; this is a kind of instability weaker than the preceding one.

If the recurrence relation is stable $`x_{max}`$ can be chosen large at will, with no effect on the precision of $`U^{(\alpha )}(1)`$. In the case of instability with $`A>1`$, in order to obtain the result $`U^{(\alpha )}(1)`$ with a number $`E`$ of exact digits, the calculations of $`U^{(\alpha )}(x_{max}+i)`$ must be performed with a greater number of decimal digits $`C=E+x_{max}\mathrm{log}_{10}A`$. Supposing $`a_s^{(\alpha )}s!`$, a rough estimate of the number $`s_{max}`$ of terms of the series needed to obtain the sum with such a precision is $`s_{max}A^{\frac{C}{CE}}x_{max}`$. Performing the calculations with fixed precision arithmetic with $`C`$ digits, it is convenient to choose a value of $`x_{max}`$ as low as possible in order to obtain the greatest $`E`$, compatibly with the rapid increases of $`s_{max}`$ and of the computing time (see an example in section 7.2.4). Performing calculations with multiprecision arithmetic, $`C`$ can be chosen at will, and $`x_{max}`$ can be increased; a convenient choice which minimizes the estimate of $`s_{max}`$ is $`C/E1+\mathrm{ln}A`$ (see an example in section 9.4). Fixed $`C/E`$, by varying $`E`$ one sees that $`x_{max}E`$ and $`s_{max}E`$. Therefore the number of terms of the series (and even the computation time) is proportional to the number of digits of precision of the results; this is true even in the stable case provided that one chooses $`x_{max}E`$.

### 6.2 Truncated expansion in $`ϵ`$

We are interested in the results in the limit $`D4`$; therefore, defined $`ϵ=(4D)/2`$, we expand all the quantities in $`ϵ`$, truncating the expansions at the first $`n_ϵ`$ terms, and we perform all the calculations using truncated series; in this way the coefficients of all the powers of $`ϵ`$ are found numerically, including the *negative* powers. This technique is perhaps not the most efficient, but is very versatile. We have implemented in the program SYS (see section 11) the arithmetic of truncated series, so that the use of series becomes as simple as with ordinary numbers.

The time of computation of a factorial series grows approximately as $`n_ϵ^2`$, and it is mainly due to multiplications. The division between series is an operation less frequent than the multiplication as it occurs only once in the calculation of each $`a_s`$ with the recurrence relations (LABEL:systa2) or (LABEL:systb2) or in the calculation of each $`U(x)`$ using the recurrence relations (29) or (31). As effect of cancellations of terms when series are summed, and of divisions by series beginning with a non-zero power of $`ϵ`$, first or last terms of the expansions in $`ϵ`$ may be lost; in this case the calculations must be necessarily performed with an initial number $`n_ϵ`$ of terms of the expansions greater than the desired final number $`n_ϵ^{}`$ of terms of $`U(x)`$. The number $`n_ϵ`$ is chosen empirically. The loss of terms of the expansions in $`ϵ`$ is frequent when the recurrence relations (29) or (31) are used to obtain $`U(x)`$ for $`x\lambda `$; this is due to the fact that for such values of $`x`$ the expansion in $`ϵ`$ of $`p_0(x)`$ begins with a non-zero power of $`ϵ`$. Furthermore, we found that it is critical to recognize the cancellation of the first coefficient of a series used as divisor when, because of the numerical errors, the coefficient is a very small number instead of zero; in this case a numerical cutoff must be carefully used.

## 7 Applications to simple one-loop integrals

In this section we discuss in detail (as we assume the reader unfamiliar with these operator techniques) the solution of the difference equations for the one-loop vacuum and self-energy integrals.

### 7.1 One-loop vacuum integral

Defining

$$J(x)=\frac{[d^Dk]}{(k^2+m_1^2)^x},$$
(103)
we want to calculate the master integral $`J(1)`$. According the single identity of the set $`[1;\frac{[}{0}pt]00]`$ $`J(x)`$ satisfies the homogeneous difference equation

$$m_1^2(x1)J(x)(x1D/2)J(x1)=0.$$
(104)
We look for a solution121212The exact solution $`J(x)=a_0\mathrm{\Gamma }(xD/2)/\mathrm{\Gamma }(x)`$ does not have the form of Eq.(33). of this equation in the form of a factorial series (see section 4.3),

$$J(x)=\mu ^xV(x)=\mu ^x\underset{s=0}{\overset{\mathrm{}}{}}a_s𝝆^{Ks}.$$
(105)
Introducing the operators $`𝝅`$ and $`𝝆`$, multiplying by $`x`$ and using Eq.(39) we find the first canonical form of the difference equation

$$\left((\mu m_1^21)𝝆^2+((2\mu m_1^21)(𝝅1)+D/2)𝝆+\mu m_1^2𝝅(𝝅1)\right)V(x)=0.$$
(106)
The characteristic equation Eq.(45), $`f_2(\mu )=0`$, gives the value $`\mu =1/m_1^2`$. Substituting this value of $`\mu `$ in Eq.(106) one obtains

$$\left((𝝅1+D/2)𝝆+𝝅(𝝅1)\right)V(x)=0.$$
(107)
The indicial equation Eq.(49), $`f_1(𝝅=1+K)=0`$, gives the value $`K=D/2`$. Using the recurrence relation (LABEL:systa2) one finds

$$a_s=\frac{\mathrm{\Gamma }(s+D/2)\mathrm{\Gamma }(s+D/2+1)}{\mathrm{\Gamma }(D/2)\mathrm{\Gamma }(D/2+1)\mathrm{\Gamma }(s+1)}a_0.$$
(108)
From Eq.(71) one finds that $`a_0=(m_1^2)^{D/2}`$. The final result is

$$J(x)=(m_1^2)^{D/2x}\underset{s=0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(s+D/2)\mathrm{\Gamma }(s+D/2+1)}{\mathrm{\Gamma }(D/2)\mathrm{\Gamma }(D/2+1)\mathrm{\Gamma }(s+1)}\frac{\mathrm{\Gamma }(x+1)}{\mathrm{\Gamma }(x+1+D/2+s)}.$$
(109)
The coefficients $`a_s`$ behave for large $`s`$ as $`a_s/s!s^{D1}`$, therefore the term of the factorial series behaves as $`s^{D/21x}`$, so that the series has abscissa of convergence $`\lambda =D/2`$. This value signals the divergence of the integral (103). The recurrence relation (104), being of order one, is numerically stable. Therefore it is possible to obtain $`J(x)`$ for a large integer $`xD/2`$ by summing a few terms of the factorial series and to use the recurrence relation to obtain $`J(1)`$. For a numerical example in the limit $`D4`$ see section 7.2.4.

### 7.2 One-loop self-energy integral

Let us consider the one-loop self-energy integral

$$I(1)=\frac{[d^Dk]}{D_1D_2},$$
(110)
where $`D_1=k^2+m_1^2`$ and $`D_2=(pk)^2+m_2^2`$. A simple inspection of the system of identities $`[1\mathrm{}2;\frac{[}{0}pt]0\mathrm{}10]`$ shows that there are three master integrals, $`[d^Dk]/D_1D_2`$, $`[d^Dk]/D_1`$ and $`[d^Dk]/D_2`$. We must find three difference equations for

$$I(x)=\frac{[d^Dk]}{D_1^xD_2},J(x)=\frac{[d^Dk]}{D_1^x},K(x)=\frac{[d^Dk]}{D_2^x}.$$
(111)
The solution of the system made up of the set of identities $`[1\mathrm{}2;\frac{[}{0}pt]00\mathrm{}1]`$ gives the system of difference equations

$$\begin{array}{c}(xD)I(x2)+(p^2+m_2^2m_1^2)(2xD1)I(x1)+R^2(p^2,m_1^2,m_2^2)(x1)I(x)\hfill \\ \hfill +J(x1)((D/21)(p^2+m_2^2+m_1^2)(x2)(p^2+m_2^2))/m_1^2=0,\end{array}$$
(112)
$$m_1^2(x1)J(x)(x1D/2)J(x1)=0,$$
(113)
$$m_2^2(x1)K(x)(x1D/2)K(x1)=0,$$
(114)
where $`R`$ was defined in Eq.(79). The difference equation for $`J(x)`$ and $`K(x)`$ has been solved in the preceding section, so that we consider here only the second-order difference equation for $`I(x)`$. We discuss separately the massive case and the cases where one of the masses $`m_1`$ or $`m_2`$ is zero.

#### 7.2.1 Case 1: $`m_10`$, $`m_20`$

Following section 4.3 we substitute $`I(x)=\mu ^xV(x)`$ and $`J(x)=\mu ^xW(x)`$ in Eq.(112), multiply it by $`x(x1)`$ and insert the operators $`𝝅`$ and $`𝝆`$. We obtain the canonical form of the equation

$$\left(f_3(\mu )𝝆^3+f_2(𝝅,\mu )𝝆^2+f_1(𝝅,\mu )𝝆+f_0(𝝅,\mu )\right)V(x)=\mu \left(g_3(𝝅)𝝆^3+g_2(𝝅)𝝆^2+g_1(𝝅)𝝆\right)W(x);$$
(115)
$`f_i`$ and $`g_i`$ are polynomials in $`𝝅`$, $`\mu `$, $`D`$, $`p^2`$, $`m_1^2`$ and $`m_2^2`$, whose explicit expressions are not shown for brevity. The characteristic equation $`f_3(\mu )=0`$ has the two different roots

$$\mu =\mu _\pm =\frac{1}{(p\pm im_2)^2+m_1^2}$$
(116)
(we define $`p=\sqrt{p^2}`$, and if $`p^2<0`$ then $`p=i\sqrt{p^2}`$). The solution of the difference equation (112) is

$$I(x)=\eta _+I_+^{HO}(x)+\eta _{}I_{}^{HO}(x)+I^{NH}(x),$$
(117)
where $`I_\pm ^{HO}`$ are the solutions of the homogeneous equation corresponding to $`\mu _\pm `$ and $`I^{NH}`$ is one solution of the nonhomogeneous equation.

Let us now consider the homogeneous solution $`I_{}^{HO}`$. We write $`I_{}^{HO}(x)=\mu _{}^xV_{}^{HO}(x)`$ and we look for the coefficients of the expansion in factorial series of $`V_{}^{HO}(x)`$

$$V_{}^{HO}(x)=\underset{s=0}{\overset{\mathrm{}}{}}a_s^{}𝝆^{K_{}s}.$$
(118)
Substituting $`\mu =\mu _{}`$, the homogeneous part of the canonical form (115) becomes

$$\left(f_2^{}(𝝅)𝝆^2+f_1^{}(𝝅)𝝆+f_0^{}(𝝅)\right)V_{}^{HO}(x)=0,$$
(119)
$$\begin{array}{cc}\hfill f_2^{}(𝝅)=& 2im_2p(2𝝅+D5),\hfill \\ \hfill f_1^{}(𝝅)=& \left((p^2+m_1^2m_2^2)(𝝅+D3)+2im_2p(3𝝅4)\right)(𝝅1),\hfill \\ \hfill f_0^{}(𝝅)=& (p^2+m_1^2m_2^2+2im_2p)𝝅(𝝅1)^2;\hfill \end{array}$$
(120)
the indicial equation $`f_2^{}(K_{}+2)=0`$ gives $`K_{}=(1D)/2`$. Fixed $`a_0^{}=1`$, the other coefficients $`a_s^{}`$ can be found using the recurrence relation (LABEL:systa2); the behaviour of $`a_s^{}`$ for large $`s`$ can be determined by considering the recurrence relation as a difference equation in $`s`$ for $`a_s^{}`$ and solving it; one finds

$$a_s^{}/s!C_1s^{(3D7)/2}+C_2B^s/s,$$
(121)
where $`B=\mu _{}/(\mu _{}\mu _+)`$, while $`C_1`$ and $`C_2`$ are constants. If $`|B|>1`$ the series (118) never converges (in fact the condition (102) is satisfied); if $`|B|1`$, the series converges with abscissa of convergence $`\lambda =D2`$. The solution $`I_+^{HO}`$ can be obtained in analogous way with the replacement $`im_2pim_2p`$.

Let us now consider the nonhomogeneous solution $`I^{NH}(x)`$. We write $`I^{NH}(x)=\mu ^xV^{NH}(x)`$. The value of $`\mu `$ must be taken from the expansion in factorial series of $`J(x)`$ obtained in section 7.1: $`\mu =1/m_1^2`$. Substituting this value of $`\mu `$ in the canonical form, Eq.(115) becomes

$$\left(\widehat{f}_3(𝝅)𝝆^3+\widehat{f}_2(𝝅)𝝆^2+\widehat{f}_1(𝝅)𝝆+\widehat{f}_0(𝝅)\right)V^{NH}(x)=\left(\widehat{g}_3(𝝅)𝝆^3+\widehat{g}_2(𝝅)𝝆^2+\widehat{g}_1(𝝅)𝝆\right)W(x),$$
(122)
$$\begin{array}{cc}\hfill \widehat{f}_3(𝝅)=& (p^2+m_2^2)^2,\hfill \\ \hfill \widehat{f}_2(𝝅)=& \left(3(p^2+m_2^2)^2+2m_1^2(p^2m_2^2)\right)𝝅5(p^2+m_2^2)^2+(D5)m_1^2(p^2m_2^2),\hfill \\ \hfill \widehat{f}_1(𝝅)=& \left[\left((p^2+m_2^2)^2+m_1^2(p^2m_2^2)\right)(3𝝅4)+m_1^2(p^2+m_1^2m_2^2)(𝝅+D3)\right](𝝅1),\hfill \\ \hfill \widehat{f}_0(𝝅)=& R^2(p^2,m_1^2,m_2^2)𝝅(𝝅1)^2,\hfill \\ \hfill \widehat{g}_3(𝝅)=& p^2+m_2^2,\hfill \\ \hfill \widehat{g}_2(𝝅)=& (p^2+m_2^2)(2𝝅3)(D/2)(p^2+m_1^2+m_2^2)+m_1^2,\hfill \\ \hfill \widehat{g}_1(𝝅)=& (p^2+m_2^2)(𝝅1)^2\left((D/2)(p^2+m_1^2+m_2^2)m_1^2\right)(𝝅1),\hfill \end{array}$$
(123)
where $`J(x)=(m_1^2)^xW(x)`$, $`W(x)=_{s=0}^{\mathrm{}}b_s𝝆^{K_bs}`$, $`K_b=D/2`$, and the coefficients $`b_s`$ are the coefficients $`a_s`$ defined in Eq.(108).

The right-hand side of Eq.(122) can be written as $`_{s=0}^{\mathrm{}}c_s𝝆^{K_cs}`$, where $`K_c=K_b+3`$, and the expression of $`c_s`$ is given by

$$c_s=b_s\widehat{g}_3(K_cs)+b_{s1}\widehat{g}_2(K_cs)+b_{s2}\widehat{g}_1(K_cs).$$
(124)
Then, the coefficients $`\widehat{a}_s`$ of the expansion of $`V^{NH}(x)=_{s=0}^{\mathrm{}}\widehat{a}_s𝝆^{K_bs}`$ can be found by using the recurrence relation (LABEL:systb2); the large-$`s`$ behaviour is

$$\widehat{a}_s/s!C_1s^{(3D6)/2}+C_2s^{D2}+C_3B_+^s/\sqrt{s}+C_4B_{}^s/\sqrt{s},$$
(125)
where $`B_\pm =m_1^2/(m_1^2\mu _\pm )`$ and $`C_i`$ are constants. If $`|B_+|>1`$ or $`|B_{}|>1`$ the expansion of $`V^{NH}`$ never converges (the condition (102) is satisfied); if $`|B_+|<1`$ and $`|B_{}|<1`$, the series converges, and the abscissa of convergence is $`\lambda =\mathrm{max}(D2,D/21)`$ if $`C_10`$ and $`C_20`$.

Now we must determine the constants $`\eta _+`$ and $`\eta _{}`$; we compare the large-$`x`$ behaviours of the solutions $`I_\pm ^{HO}`$ and $`I^{NH}`$,

$`I_\pm ^{HO}(x)`$ $`\mu _\pm ^xx^{(1D)/2},`$ (126)
$`I^{NH}(x)`$ $`\{\begin{array}{cc}(p^2+m_2^2)^1(m_1^2)^{D/2x}x^{D/2}\hfill & p^2m_2^2,\hfill \\ (2D)(4m_2^2)^1(m_1^2)^{D/2x}x^{D/2}\hfill & p^2=m_2^2,\hfill \end{array}`$ (127)
obtained from the first term of factorial series, with the large-$`x`$ behaviour of $`I(x)`$. We point out that these behaviours and the values of $`\eta _\pm `$ later inferred are rigorously valid only if all the factorial series have finite abscissa of convergence; however the deduction may be extended to the cases where the series have $`\lambda =\mathrm{}`$ with the help of the integral representations of section 8.6.2.

If $`p^2m_2^2`$ we are below the deformation threshold, and we can write $`I(x)`$ as radial integral (see Eq.(80))

$$I(x)=\frac{1}{\mathrm{\Gamma }(D/2)}_0^{\mathrm{}}\frac{dk^2(k^2)^{D/21}}{(k^2+m_1^2)^x}\frac{Z}{pk}\mathrm{\Phi }(Z^2);$$
(128)
the large-$`x`$ behaviour is given by Eq.(71)

$$I(x)(m_1^2)^{D/2x}x^{D/2}(p^2+m_2^2)^1.$$
(129)
The comparison of Eq.(129) with Eqs.(126)-(127) shows that, as $`\mu _\pm 1/m_1^2`$, $`I_\pm ^{HO}`$ do not contribute to $`I`$, and that only $`I^{NH}`$ contributes, so that

$$\eta _+=\eta _{}=0(p^2>m_2^2,m_10,m_20).$$
(130)
At the deformation threshold $`p^2=m_2^2`$ one finds from Eq.(93) for large $`x`$

$$I(x)(m_1^2)^{D/21/2x}x^{(1D)/2}\frac{\sqrt{\pi }}{2m_2};$$
(131)
as $`\mu _{}=1/m_1^2`$ and $`K_{}=D/2+1/2`$, Eq.(98) is satisfied, so that $`I_{}^{HO}`$ contributes to $`I`$, and $`\eta _{}`$ is different from zero. Comparing Eq.(131) and Eqs.(126)-(127) one finds

$$\eta _{}=m_1^{D1}m_2^1\sqrt{\pi }/2,\eta _+=0(p^2=m_2^2,m_10,m_20).$$
(132)
If $`p^2<m_2^2`$ we are above the deformation threshold, and the radial path of integration turns around the point $`k^2=(pim_2)^2`$. We split $`I(x)=I_a(x)+I_b(x)`$ (see Eq.(81)), where $`I_a`$ is given by Eq.(128) and $`I_b`$ is

$$I_b(x)=\frac{1}{\mathrm{\Gamma }(D/2)}_{(pim_2)^2}^0\frac{dk^2(k^2)^{D/21}}{(k^2+m_1^2)^x}\left(\frac{Z}{pk}\mathrm{\Phi }(Z^2)\frac{Z^1}{pk}\mathrm{\Phi }(Z^2)\right).$$
(133)
Noting that for large $`x`$ the contribution to this integral comes from the neighbourhood of the singularity, one finds

$$I_b(x)\mu _{}^xx^{(1D)/2}\frac{\sqrt{\pi }}{m_2}\left(\frac{im_2/p}{\mu _{}}\right)^{(D1)/2},$$
(134)
so that

$$\eta _{}=\frac{\sqrt{\pi }}{m_2}\left(\frac{im_2/p}{\mu _{}}\right)^{(D1)/2},\eta _+=0(p^2<m_2^2,m_10,m_20).$$
(135)
#### 7.2.2 Case 2: $`m_10`$, $`m_2=0`$

In this case we have $`\mu _+=\mu _{}=1/(p^2+m_1^2)`$. The second-order difference equation (112) simplifies to the first-order equation

$$(Dx1)I(x1)+(x1)(p^2+m_1^2)I(x)(x1)J(x)=0,$$
(136)
whose solution is of the kind

$$I(x)=\eta I^{HO}(x)+I^{NH}(x).$$
(137)
In a manner analogous to the previous case, one obtains the expansion in factorial series

$$I^{HO}(x)=(p^2+m_1^2)^x\underset{s=0}{\overset{\mathrm{}}{}}a_s𝝆^{2Ds};$$
(138)
the large-$`s`$ behaviour of the coefficients is $`a_s/s!s^{2D5}`$, so that the series has abscissa of convergence $`\lambda =D2`$. Considering the nonhomogeneous solution,

$$I^{NH}(x)=(m_1^2)^x\underset{s=0}{\overset{\mathrm{}}{}}\widehat{a}_s𝝆^{D/2s},$$
(139)
the coefficients have the large-$`s`$ behaviour

$$\widehat{a}_s/s!C_1s^{(3D6)/2}+C_2s^{D2}+B_0^ss^{(D2)/2},$$
(140)
where $`B_0=1+m_1^2/p^2`$ and $`C_i`$ are constants. If $`p^2>m_1^2/2`$, $`|B_0|>1`$ and the series does not converge; if $`p^2m_1^2/2`$ the series converges with abscissa of convergence $`\lambda =\mathrm{max}(D2,D/21)`$, if $`C_10`$ and $`C_20`$.

Now we determine the constant $`\eta `$. The large-$`x`$ behaviour of the solutions is given by the first term of the factorial series

$`I^{HO}(x)`$ $`(p^2+m_1^2)^xx^{2D},`$ (141)
$`I^{NH}(x)`$ $`\{\begin{array}{cc}p^2(m_1^2)^{D/2x}x^{D/2}+\mathrm{?}\hfill & p^2>0,\hfill \\ p^2(m_1^2)^{D/2x}x^{D/2}\hfill & p^2<0.\hfill \end{array}`$ (142)
In the whole region $`p^2>0`$ the factorial series expansion of $`I^{NH}`$ never converges so that the large-$`x`$ behaviour shown above is valid only in asymptotic sense. Using the integral representations of section 8.6.2 one can show that if $`p^2>0`$ the large-$`x`$ behaviour of $`I^{NH}`$ contains an additional contribution, exponentially small in comparison with the main contribution

$$I^{NH}(x)p^2(m_1^2)^{D/2x}x^{D/2}+\mathrm{\Gamma }(D/21)(p^2)^{1D/2}(p^2+m_1^2)^{D2x}x^{2D}.$$
(143)
The large-$`x`$ behaviour of $`I(x)`$ can be deduced from Eq.(128) and Eq.(133), as in the massive case. But if $`p^2>0`$ there is a difference: the integrand of Eq.(128) is discontinuous in the point $`k^2=p^2`$, in fact

$$Z=\frac{p^2+k^2|p^2k^2|}{2pk}=\{\begin{array}{cc}k/pk<p,\hfill & \\ p/kkp.\hfill & \end{array}$$
(144)
The presence of this discontinuity gives rise to an exponentially small additional contribution proportional to $`(p^2+m_1^2)^x`$ in the large-$`x`$ behaviour of $`I(x)`$. So we split the integral over $`k^2`$ into two parts,

$$\begin{array}{c}I(x)=\frac{1}{\mathrm{\Gamma }(D/2)}_0^{\mathrm{}}\frac{dk^2(k^2)^{D/21}}{(k^2+m_1^2)^x}\frac{1}{p^2}\mathrm{\Phi }\left(\frac{k^2}{p^2}\right)\hfill \\ \hfill +\frac{1}{\mathrm{\Gamma }(D/2)}_{p^2}^{\mathrm{}}\frac{dk^2(k^2)^{D/21}}{(k^2+m_1^2)^x}\left(\frac{1}{k^2}\mathrm{\Phi }\left(\frac{p^2}{k^2}\right)\frac{1}{p^2}\mathrm{\Phi }\left(\frac{k^2}{p^2}\right)\right),\end{array}$$
(145)
each one having a different large-$`x`$ behaviour. One finds that the large-$`x`$ behaviour of $`I(x)`$ is identical to that of $`I^{NH}(x)`$ for $`p^2>0`$, Eq.(143), so that

$$\eta =0(p^2>0,m_10,m_2=0).$$
(146)
If $`p^2<0`$ the large-$`x`$ behaviour of $`I(x)`$ can obtained from Eq.(133) (with $`m_2=0`$). One finds the same result found for $`p^2>0`$. Therefore for $`p^2<0`$ the exponentially small contribution is present in the large-$`x`$ behaviour of $`I(x)`$ but not in that of $`I^{NH}(x)`$: it must come from $`\eta I^{HO}(x)`$. Therefore the constant $`\eta `$ is

$$\eta =\mathrm{\Gamma }(D/21)\left((p^2+m_1^2)/(ip)\right)^{D2}(p^2<0,m_10,m_2=0).$$
(147)
#### 7.2.3 Case 3: $`m_1=0`$, $`m_20`$

In this case the denominator raised to $`x`$ has zero mass, and the term containing $`J`$ in Eq.(112) disappears, so that the difference equation becomes homogeneous; the solution is of the kind

$$I(x)=\eta _+I_+^{HO}(x)+\eta _{}I_{}^{HO}(x).$$
(148)
$`I_\pm ^{HO}(x)`$ are the same functions considered in section 7.2.1 with $`m_1=0`$, and their large-$`x`$ behaviour was given in Eq.(126). But the large-$`x`$ behaviour of $`I(x)`$ cannot be found as in the massive case; in fact if $`m_1=0`$ the integrals (128) and (133) are strongly divergent in $`k^2=0`$ for large $`x`$ and only dimensional regularization can give a finite value to $`I(x)`$. The large-$`x`$ behaviour of $`I(x)`$ must be determined by other methods, for example by writing $`I(x)`$ as integral over one Feynman parameter and using the saddle-point method; one finds

$$\eta _\pm =\frac{\sqrt{\pi }}{2m_2\mathrm{sin}(\pi D/2)}\left(ip\mu _\pm /m_2\right)^{(1D)/2}(m_1=0,m_20).$$
(149)
#### 7.2.4 Numerical example

Let us describe in some detail the calculation of $`I(1)`$, $`J(1)`$ and $`K(1)`$ in the case $`p^2=m_1^2=m_2^2=1`$. As $`m_1=m_2`$ then $`K(x)=J(x)`$, and we consider only $`J(x)`$ and $`I(x)`$. The root of characteristic equation of the difference equation for $`J(x)`$ (Eq.(113)) is $`\mu =1`$; the roots of the characteristic equation of the homogeneous part of the equation for $`I(x)`$ (Eq.(112)) are $`\mu _{}=1`$ and $`\mu _+=1/3`$. $`I(x)`$ is written as sum of the homogeneous and nonhomogeneous part, Eq.(117); from Eq.(132) one finds $`\eta _+=0`$, so that $`I_+^{HO}`$ does not contribute. Therefore we have to evaluate three functions: $`J(x)`$, $`I_{}^{HO}(x)`$ and $`I^{NH}(x)`$ for $`x=1`$ and $`D4`$. The factorial series expansions of all the three functions have abscissa of convergence $`\lambda 2`$ for $`D4`$ and do not converge for $`x=1`$. Following the procedure of section 6 we evaluate the series for a large value $`x=x_{max}`$ and, as the equation for $`I(x)`$ is of second order, even for $`x=x_{max}+1`$. The convergence becomes faster by increasing $`x_{max}`$, but, as $`A=|\mu _{}/\mu _+|=3>1`$, the recurrence relation for $`I(x)`$ Eq.(112), rewritten in order to obtain $`I(x2)`$ from $`I(x1)`$ and $`I(x)`$, is unstable. Each application of the recurrence relation increases the error of the value of $`I(x)`$ of about a factor $`3`$. Therefore we must choose a value of $`x_{max}`$ which is a compromise between speed of computation and loss of precision in the result. We choose $`x_{max}=8`$. Calculations are performed by setting $`D=42ϵ`$ and expanding in $`ϵ`$, truncating the expansions at the first three terms. Doing the calculations with 19 digits, the convergence is attained for all the factorial series for $`x=8`$ (and consequently for $`x=9`$) in about 4000 terms. The values for $`x<8`$ are calculated by using repeatedly the recurrence relations (112) and (113). The application of the unstable recurrence relation enlarges the error of about a factor $`3^72000`$, corresponding to a loss of about 3 significant digits in $`I(1)`$. Values of $`J(x)`$, $`I_{}^{HO}(x)`$ and $`I^{NH}(x)`$ are shown in Table 1, with only the first two terms of the expansion in $`ϵ`$, and the coefficients with only 6 digits to save space. The value of $`J(1)`$ calculated agrees with the exact result $`ϵ^1+\gamma 1+O(ϵ)`$ within the precision of the calculation ($`\gamma `$ is the Euler’s constant). Inserting $`I_{}^{HO}(1)`$, $`I^{NH}(1)`$ and the value $`\eta _{}=\sqrt{\pi }/2`$ (from Eq.(132)) in Eq.(117) one finds

$$I(1)=(12\times 10^{17})ϵ^10.3910150291357503+O(ϵ),$$
(150)
with an error of $`4\times 10^{16}`$ on the constant term in comparison with the exact result $`ϵ^1+2\pi /\sqrt{3}\gamma +O(ϵ)`$.

A numerical value of $`\eta _{}`$ can be obtained independently, using the identity (100)

$$I(0)=\eta _{}I_{}^{HO}(0)+I^{NH}(0);$$
(151)
inserting the values of $`I(0)=K(1)=J(1)`$, $`I_{}^{HO}(0)`$ and $`I^{NH}(0)`$ one finds

$$\eta _{}=0.8862269254527578+5\times 10^{15}ϵ+O(ϵ^2),$$
(152)
with an error of $`2\times 10^{16}`$ on the constant term. Using this numerical value in the place of the analytical value one obtains a value of $`I(1)`$ with the same precision as Eq.(150).

In Table 2 we show for different choices of $`x_{max}`$ the number of terms of series necessary to evaluate $`J(8)`$, $`I_{}^{HO}(8)`$ and $`I^{NH}(8)`$ with 19 digits of precision and the obtained values of the finite part of $`I(1)`$; it is evident that by increasing $`x_{max}`$ the series converges faster, but the precision degrades because of the increasing number of applications of the unstable recurrence relation.

## 8 Solutions of difference equations by means of Laplace’s transformation

The expansion in factorial series is certainly the most direct method of solution of difference equations with polynomial coefficients; however, for some values of masses and external momenta of the diagram, as we have seen in the above example, the abscissa of convergence of factorial series may become infinite. In this case the factorial series become divergent for every value of $`x`$ and therefore useless, so that another method of solution must be used: the Laplace’s transformation method.

This method is described in section 8.1, and applied to the systems of difference equations in section 8.2, 8.3 and 8.4. Techniques used for integrating the differential equations obtained from the application of the method are described in section 8.5. The application to simple one-loop integrals is shown in section 8.6.

### 8.1 Transformation of a difference equation

Let us consider the difference equation

$$p_0(x)U(x)+p_1(x)U(x+1)+\mathrm{}+p_N(x)U(x+N)=0,$$
(153)
where $`p_i(x)`$ are polynomials in $`x`$ of maximum degree $`P`$. The Laplace’s transformation method consists in the substitution

$$U(x)=_l𝑑tt^{x1}v(t),$$
(154)
where $`l`$ is a line of integration suitably determined and where $`v(t)`$ is found from a certain differential equation. Writing the coefficients as

$$p_k(x)=A_{k0}+\underset{i=1}{\overset{P}{}}A_{ki}\underset{j=0}{\overset{i1}{}}(x+k+j),$$
(155)
substituting Eq.(154) in Eq.(153) and integrating by parts one finds

$$\underset{k=0}{\overset{N}{}}p_k(x)U(x+k)=_l𝑑tt^{x1}\underset{i=0}{\overset{P}{}}\mathrm{\Phi }_i(t)(t)^iv^{(i)}(t)+[I(x,t)]_l,$$
(156)
where

$$\mathrm{\Phi }_i(t)=\underset{k=0}{\overset{N}{}}A_{ki}t^k,$$
(157)
$$I(x,t)=\underset{i=0}{\overset{P1}{}}(1)^iv^{(i)}(t)\underset{m=0}{\overset{P1i}{}}\left(\frac{d}{dt}\right)^m\left(\mathrm{\Phi }_{m+i+1}(t)t^{x+m+i}\right).$$
(158)
Eq.(154) provides a solution of the difference equation (153) if $`v(t)`$ is a solution of the differential equation

$$\underset{i=0}{\overset{P}{}}\mathrm{\Phi }_i(t)(t)^iv^{(i)}(t)=0,$$
(159)
provided that the line of integration $`l`$ be chosen so that $`I(x,t)`$ has the same value at each endpoint of the line, if the line is open. Note that the difference equation (153) has order $`N`$ with coefficients of degree $`P`$, while the differential equation (159) has order $`P`$ with coefficients of degree $`N+P`$.

The singular points of this differential equation are $`0`$, $`\mathrm{}`$ and the zeros $`t_i`$ (of multiplicity $`m_i`$) of the *characteristic equation*

$$\mathrm{\Phi }_P(t)=0;$$
(160)
these points turn to be always regular singular points in the case of the differential equations encountered in this work.

We choose as lines of integration the lines which begin in the origin and end in one of the singular points $`t_i`$. This is a convenient choice. One can show that $`I(x,t)=0`$ at each endpoint of such lines, if the integral over $`t`$ of Eq.(154) is finite. Under these conditions, $`U(x)`$ is a solution of the difference equation.

We can construct a set of $`_im_i=N`$ functions $`U_{ij}(x)`$ which form a fundamental system of solutions of the difference equation (153) by defining

$$U_{ij}(x)=_{l_i}𝑑tt^{x1}v_{ij}(t)j=1,\mathrm{},m_i,$$
(161)
where $`v_{ij}(t)`$ is one of the $`m_i`$ solutions of the differential equation singular in $`t_i`$, and $`l_i`$ is the line which begins in $`t=0`$ and ends in $`t=t_i`$.

It is important to note that the characteristic equation (160) of the differential equation and the characteristic equation of the difference equation (153) turn out to be identical, so that we can readily identify the singular points $`t_i`$ of the differential equation with the solutions $`\mu _i`$ of Eq.(45).

Now we consider the nonhomogeneous equation

$$\underset{k=0}{\overset{N}{}}p_k(x)U(x+k)=\underset{k=0}{\overset{N^{}}{}}q_k(x)T(x+k),$$
(162)
where $`q_k(x)`$ are polynomials and $`T(x)`$ is a solution of some difference equation. A particular solution $`U^{NH}(x)`$ is found by substituting into the equation

$$T(x)=_{l_T}𝑑tt^{x1}w(t),U^{NH}(x)=_{l_T}𝑑tt^{x1}v_{NH}(t),$$
(163)
where $`l_T`$ is a known line of integration and $`w(t)`$ is a known function, solution of a differential equation analogous to Eq.(159), obtained from the difference equation satisfied by $`T(x)`$. Provided that $`I(x,t)=0`$ at the endpoints of $`l_T`$, one finds the nonhomogeneous differential equation

$$\underset{i=0}{\overset{P}{}}\mathrm{\Phi }_i(t)(t)^iv_{NH}^{(i)}(t)=\underset{i=0}{\overset{P^{}}{}}\mathrm{\Psi }_i(t)(t)^iw^{(i)}(t),$$
(164)
whose solution gives $`v_{NH}`$. Then $`U^{NH}`$ is found using Eq.(163).

### 8.2 Transformation of the system of difference equations

In section 3.2 we described the construction of the system of difference equations. The algorithm used yields a system in triangular form; applying to this system the Laplace’s transformation one obtains a system of differential equations with the same triangular structure and the same ease of solution. There is, however, a complication: the $`l`$th difference equation

$$\underset{k=0}{\overset{N}{}}p_k(x)U_l(x+k)=\underset{j=0}{\overset{l1}{}}\underset{k=0}{\overset{N_j^{}}{}}q_{jk}(x)U_j(x+k)$$
(165)
of order $`N`$ and with polynomial coefficients of degree $`P`$ is transformed, using the Laplace’s transformation, into a differential equation of order $`P`$

$$\underset{i=0}{\overset{P}{}}\mathrm{\Phi }_i(t)(t)^iv_l^{(i)}(t)=\underset{j=0}{\overset{l1}{}}\underset{i=0}{\overset{P_j^{}}{}}\mathrm{\Psi }_{ij}(t)(t)^iv_j^{(i)}(t)$$
(166)
and coefficients of degree $`N+P`$. As effect of the particular algorithm of construction and solution of the system of identities, $`N`$ is usually small (typically $`14`$) while $`P`$ may be large ($`330`$). Therefore the differential equation may have a high order. Calculations in some test cases have shown that solution of high order equations slows down the calculations and is source of undesired numerical errors, so that, if possible, it is better to avoid it.

We have discovered that this difficulty can be overcome by applying the Laplace’s transformation to the identities obtained by integration-by-parts *before* the insertion in the system of identities, building a system of “transformed” identities between the “transformed” integrals $`v(t)`$ instead of the real integrals $`U(x)`$. The solution of the system of transformed identities will provide a system of differential equations of smaller order; a small price to pay is the (possible) appearance of spurious singular points in the equations so obtained (see section 8.5.2). More in detail, a generic integral is transformed into

$$U_{ni\alpha \beta }(x)=[d^Dk_1]\mathrm{}[d^Dk_{N_d}]\frac{_{j=1}^{N_{sp}n}(pk\text{ irred.})_j^{\beta _j}}{D_{i_1}^{x+\alpha _1}D_{i_2}^{\alpha _2}\mathrm{}D_{i_n}^{\alpha _n}}=_l𝑑tt^{x1}v_{ni\alpha \beta }(t),$$

where the line $`l`$ is unspecified; the values of the functions $`I(x,t)`$ at the endpoints of the line are always zero because of dimensional regularization. Therefore, a generic integration-by-parts identity,

$$\underset{ni\alpha \beta }{}(xr_{ni\alpha \beta }+s_{ni\alpha \beta })U_{ni\alpha \beta }(x)=0,$$

($`r_{ni\alpha \beta }`$ and $`s_{ni\alpha \beta }`$ are independent of $`x`$) becomes the transformed identity

$$\underset{ni\alpha \beta }{}\left((s_{ni\alpha \beta }r_{ni\alpha \beta }\alpha _1)v_{ni\alpha \beta }(t)tr_{ni\alpha \beta }\frac{dv_{ni\alpha \beta }}{dt}(t)\right)t^{\alpha _1}=0$$
(167)
which is a differential equation between the functions $`v_{ni\alpha \beta }(t)`$.

### 8.3 Construction of the system of differential equations

We want to build a triangular system of differential equations

$$\underset{i=0}{\overset{P_l}{}}p_{il}(t)v_{ml}^{(i)}(t)=\underset{k=1}{\overset{l1}{}}\underset{j=0}{\overset{P_{lk}}{}}q_{jkl}(t)v_{mk}^{(j)}(t),l=1,\mathrm{},L_m^{},$$
(168)
between a set of “master transformed” functions $`v_{ml}(t)`$ ($`m=1,\mathrm{},N_dN_k+1`$) analogous to the system of difference equations (25) between the master functions $`U_{ml}(x)`$ discussed in section 3.2. For this reason we devise an algorithm of construction and solution of the system of identities similar to the algorithm 2:

###### Algorithm 3

Consider the algorithm 2 with the following modifications:

1. After the step 8 of algorithm 1, transform each integration-by-parts identity using Eq.(167); in the subsequent steps replace everywhere the integrals $`[d^Dk_1]\mathrm{}[d^Dk_{N_k}]W_{ni\alpha \beta }`$ with the functions $`v_{ni\alpha \beta }(t)`$.
2. The transformed identity obtained does not depend on the index $`\alpha _1`$ of $`W`$ in the step 7 of algorithm 1 because it appears as an overall factor $`t^{\alpha _1}`$; in order to restore the total number of different identities we have decided (somewhat arbitrarily) to differentiate $`\alpha _1`$ times with respect to $`t`$ each transformed identity.
3. Ignore step 4 of algorithm 2.
4. Derivatives of master transformed functions have priority of extraction lower than other generic transformed integrals.
5. Add the new entry *the greatest derivative* to the list of priorities after the entry 9(b)vi of algorithm 1.

With a suitable choice of the parameters $`a_i`$ and $`b_i`$ (see the end of section 3.2), by means of this algorithm we can identify the master transformed functions $`v_{ml}(t)`$ as the functions which satisfy equations of non-zero order, and we can work out a set of differential equations among them. If each function $`v_{ml}(t)`$ corresponds to an integral with a different combination of denominators, the system is obtained directly in the triangular form (168); if, on the contrary, there are different master transformed functions $`v_{m,l+1}(t)`$,…,$`v_{m,l+G}(t)`$ corresponding to integrals containing the same combination of denominators, the algorithm provides a set of $`G`$ simultaneous differential equations containing all these $`G`$ functions, which are conveniently transformed into triangular form using a procedure quite analogous to that described in section 3.2. As final result one obtains a system of differential equations with the triangular form (168). It is important to note that the functions

$$F_{ml}(x)=_l𝑑tt^{x1}v_{ml}(t)$$
(169)
are not necessarily identical to the master functions $`U_{ml}(x)`$ defined in Eq.(22). While $`v_{ml}(t)`$ is a master transformed function and satisfies a differential equation of non-zero order, $`F_{ml}(x)`$ satisfies a difference equation of different order, which may be zero; if so, $`F_{ml}(x)`$ is not a *master* function. This fact frequently occurs when many master integrals have the same combination of denominators. For these reasons the number and the structure of the master transformed functions $`v_{ml}(t)`$ must be found independently of the master functions $`U_{ml}(x)`$. See section 9.4 for an example.

### 8.4 Correspondence with factorial series and initial conditions

Let $`U(x)=_\alpha U_\alpha (x)`$ be the solution of a generic difference equation; the initial conditions for the integration of the differential equation can be determined by comparing the integral representation of $`U_\alpha (x)`$ with the factorial series expansion

$$U_\alpha (x)=_0^{\mu _\alpha }𝑑tt^{x1}v_\alpha (t)=\mu _\alpha ^x\underset{s=0}{\overset{\mathrm{}}{}}a_{s\alpha }𝝆^{K_\alpha ^Fs}.$$
(170)
The behaviour of $`v_\alpha (t)`$ near the singular point $`t=\mu _\alpha `$

$$v_\alpha (t)A_{0\alpha }(\mu _\alpha t)^{K_\alpha ^L},t\mu _\alpha ,$$
(171)
can be deduced from the known behaviour of $`U_\alpha (x)`$ for large $`x`$. Substituting Eq.(171) in Eq.(170), integrating over $`t`$ and comparing the result with the large-$`x`$ leading behaviour of the first term of the series $`U_\alpha (x)\mu _\alpha ^xa_{0\alpha }x^{K_\alpha ^F}`$ one finds $`K_\alpha ^L`$ and $`A_{0\alpha }`$:

$$K_\alpha ^L=K_\alpha ^F1,A_{0\alpha }=a_{0\alpha }/(\mu _\alpha ^{K_\alpha ^L}\mathrm{\Gamma }(K_\alpha ^L+1)).$$
(172)
If necessary, the subsequent coefficients of the series expansion of $`v_\alpha (t)`$ can be deduced in the same way from the coefficients $`a_{s\alpha }`$.

### 8.5 Integrating differential equations

The choice of a effective numerical method of integration of the differential equations obtained by applying Laplace’s transformation is not simple. In fact we must consider that:

* The initial and final point of the path of integration $`l`$ are singular points of the solutions; the numerical methods usually used to integrate differential equations (for example, the Runge-Kutta method) cannot be used in singular points.
* Master functions, coefficients of differential equations and (sometimes) singular points depend on $`D`$ and are represented by truncated series in $`ϵ`$.
* The method must be able to provide very high-precision values, and the time of computation must grow linearly with the number of exact digits of the result; by using fixed-order methods (like the Runge-Kutta method) the time grows exponentially.

Therefore we have decided to solve the differential equations by using power series, expanding the general solution around a number of selected points near or on the path $`l`$ and equating the sums of the series in some intermediate points. Each power series will be evaluated inside the respective circle of convergence, so that the number of terms of the series necessary to attain a precision of $`E`$ digits in the results (and also the computation time) will be proportional to $`E`$.

#### 8.5.1 Integration over the path

The line of integration can have any shape, but for convenience it is assumed here to be the segment $`[0,\mu ]`$. The segment is subdivided into $`M`$ intervals $`[t_{M+1},t_M]`$, $`[t_M,t_{M1}]`$, …, $`[t_2,t_1]`$, where $`t_1=\mu `$ and $`t_{M+1}=0`$; the (possible) singular points placed on the segment are $`0=t_P^{si}<\mathrm{}<t_2^{si}<t_1^{si}<\mu `$. The choice of a line of integration which passes through singular points of the differential equation, instead of avoiding them, may be convenient: it allows one to check if $`v(t)`$ is regular or singular in these points and it speeds up the calculations, avoiding the use of complex numbers. Let us consider the first interval $`[t_2,t_1]`$. The solution is expanded around the point $`t_1=\mu `$ (which may be a singular point). The first coefficients of the expansions of the $`P`$ solutions corresponding to the roots of the indicial equation are obtained using Eq.(172). All the subsequent coefficients are obtained using recurrence relations obtained by substituting the expansions into the differential equation. Then the power series are evaluated in the point $`t_2=t_1r_1/2`$, where $`r_1`$ is the radius of convergence of the series and the factor $`1/2`$ has been chosen in order to minimize the total time of calculations. In the next interval the solution is expanded around the regular point $`t_2`$, the first $`P`$ coefficients of the expansion are obtained from the already known values of $`v`$ and its derivatives in $`t_2`$, the subsequent coefficients are obtained using the recurrence relations, and the power series is evaluated in a point $`t_3=t_2r_2/2`$. This process involving expansions around regular points continues for $`m`$ steps until a point $`t_m`$ is reached, placed inside the circle of convergence of the series with center the singular point $`t_1^{si}`$, at a distance less than one half of the radius $`r_1^{si}`$ of the circle. Now we expand the solution around the singular point $`t_1^{si}`$. We write the general solution as $`v(t)=_{j=1}^Pc_jv_j(t)`$, where each $`v_j`$ corresponds to one of the $`P`$ solutions of the indicial equation. The values of $`c_j`$ are found by solving the linear system131313 It is important to note that if $`m+1`$ solutions of the indicial equation coincide for $`D4`$, the Wronskian determinant of the linear system is proportional to $`(D4)^m`$; this causes a loss of $`m`$ terms in the expansions in $`D4`$ of $`c_j`$ and, consequently, in the expansions of $`v`$ in the rest of integration, and therefore in $`U(x)`$. This mishap sometimes occurred in the final singular point $`t=0`$ in the calculations of section 9.2.  $`_{j=1}^Pc_jv_j^{(i)}(t_m)=v^{(i)}(t_m)`$, $`i=1,\mathrm{},P`$. The right-hand sides contain the already known values of $`v`$ and its derivatives in $`t_m`$; the necessary values of partial solutions $`v_j(t_m)`$ are worked out by summing the corresponding expansions about $`t_1^{si}`$, whose coefficients are found using recurrence relations. Then the singular point is got over by evaluating $`v`$ (and its derivatives) in the regular point $`t_{m+1}=t_1^{si}r_1^{si}/2`$. The process is repeated until the next singular point $`t_2^{si}`$ is reached, then $`t_3^{si}`$, $`t_4^{si}`$, etc., up to the final point $`t_P^{si}=t_{M+1}=0`$.

The integration of $`v(t)`$ over $`t`$ needed to obtain $`U(x)`$ can be easily carried out by integrating the expansions in series in the corresponding intervals:

$$U(x)=_0^\mu 𝑑tt^{x1}v(t)=\underset{i=1}{\overset{M}{}}_{t_i}^{t_{i+1}}𝑑tt^{x1}\underset{j=1}{\overset{P}{}}\underset{s=0}{\overset{\mathrm{}}{}}a_s^{(i,j)}(t\overline{t}_i)^{K_{ij}+s}.$$
(173)
The integrals

$$I(j,s)=_{\overline{t}+a}^{\overline{t}+b}𝑑tt^j(t\overline{t})^{k+s}=_a^b𝑑y(y+\overline{t})^jy^{k+s}$$
(174)
which appear in Eq.(173) can be expressed in terms of incomplete Beta function (note that $`k`$ is not an integer). If $`j=0`$ the integral is immediate; if $`j`$ is positive integer the value can be efficiently computed using the recurrence relation

$$I(j,s)=I(j1,s+1)+\overline{t}I(j1,s).$$
(175)
#### 8.5.2 Singular points depending on $`D`$

The coefficients $`a_s^{(i,j)}`$ and the exponents $`K_{ij}`$ in Eq.(173) depend on $`D`$; therefore all quantities are expanded around $`D=4`$, and truncated series are used in the calculation, as described in section 6.2. A new feature, characteristic of the differential equations obtained by solving the system of transformed identities, is the appearance of spurious apparent singular points $`\overline{t}`$, not corresponding to any solution $`\mu `$ of the characteristic equations of the difference equations obtained by solving the system of ‘original’ identities. These apparent singular points are depending on $`D`$, and correspond to regular points of the solution of differential equation, in contrast with regular singular points, which are independent of $`D`$. In general the line of integration can be deformed in order to avoid these spurious singular points; but if for $`D4`$ one or more of these *mobile* points tends to one of the endpoints of the line, $`t=\mu `$ or $`t=0`$, we cannot avoid it, and we encounter difficulty in working out the solution near these coalescing points. Let us explain this fact, by considering a homogeneous differential equation with polynomial coefficients

$$\underset{i=0}{\overset{P}{}}p_i(t)v^{(i)}(t)=0,p_i(t)=t^i\underset{j=0}{\overset{g_i}{}}p_{ij}t^j,$$
(176)
and supposing for simplicity that the equation has only one apparent singular point $`t_0(D)`$ such that

$$t_0(D)=O(D4)\text{for}D4.$$
(177)
The coefficients of the expansion $`v(t)=_{s=0}^{\mathrm{}}a_st^{s+K}`$ can be found using the recurrence relation

$$a_s=\frac{_{j=1}^ma_{sj}f_j(K+sj)}{f_0(K+s)},$$
(178)
where $`a_i0`$ if $`i<0`$, $`m=\mathrm{max}_ig_i`$, $`f_j(k)=_{i=0}^Pk(k1)\mathrm{}(ki+1)p_{ij},`$ and $`K`$ is one of the roots of indicial equation $`f_0(K)=0`$ (note the analogy with the solution of a difference equation with expansions in factorial series).

As the point $`t_0`$ is a regular point, the coefficients $`a_s`$, which are functions of $`D`$, have values in $`D=4`$ finite and, in general, different from zero; unfortunately, there are problems for calculating them. As a consequence of Eq.(177), the function $`f_0`$ vanishes if $`D=4`$, and the other $`f_j`$ do not vanish (but the sum in the numerator of Eq.(178) always vanishes for $`D=4`$), so that the recurrence relation (178) turns out to be very unstable, with a degree of instability proportional to $`1/(D4)`$. Performing the calculations of $`a_s`$ using truncated expansions in $`D4`$, each iteration of Eq.(178) (with increasing $`s`$) yields one new coefficient $`a_s`$, whose expansion in $`D4`$ has a number of terms reduced by one in comparison with $`a_{s1}`$; after a few iterations, the given number of terms of the expansion in $`D4`$ is exhausted. We found that a solution of the problem is to modify the recurrence relation to

$$a_{s1}^{(n+1)}=\frac{_{j=2}^ma_{sj}^{(n)}f_j(K+sj)+a_s^{(n)}f_0(K+s)}{f_1(K+s1)},$$
(179)
so that the denominator is $`f_1`$, which does not vanish for $`D=4`$. The new recurrence relation requires as input the value of the $`s`$th coefficient before that its value is obtained; therefore Eq.(179) must be seen as part of an iterative process:

1. set $`a_s^{(0)}=0`$ for $`s=1,2,\mathrm{}`$, $`s_{max}`$;
2. apply Eq.(179) for $`s=1,2,\mathrm{}`$, $`s_{max}`$ obtaining the coefficients $`a_s^{(1)}`$;
3. repeat $`n`$ times the step 2 until $`|a_s^{(n+1)}a_s^{(n)}|=O((D4)^m)`$ for every $`s`$, where $`m`$ is the desired number of term of the expansion in $`D4`$; the convergence is guaranteed in about $`m`$ steps by the fact that $`f_0=O(D4)`$.

Analogous modifications must be made to the recurrence relation in the case of nonhomogeneous equations, or in case of two or more spurious singular points coalescing to the same endpoint. An example of equation with one mobile singular point is the third-order equation141414 The equation (180) was found by solving the system of transformed identities. Applying the Laplace’s transformation directly to the difference equation (181) one gets a higher order differential equation which does not have the mobile singular point. This equation can be also derived from Eq.(180) by writing $`(3td\mathrm{\Phi }(t)/dt+4(D4)\mathrm{\Phi }(t))/(4(D1)t(D4))=0`$ where $`\mathrm{\Phi }(t)`$ is the left-hand side of Eq.(180).

$$\begin{array}{c}2(t1)(8t+1)(4(D1)tD+4))t^3v_{\text{2}\text{h}}^{\prime \prime \prime }+(32(D1)(5D7)t^3\hfill \\ \hfill +8(D1)(25D63)t^2+(4D^2+100D192)t(D4)(9D22))t^2v_{\text{2}\text{h}}^{\prime \prime }\\ \hfill +(32(D1)(D2)(3D5)t^3+32(D1)(D3)(5D12)t^2\\ \hfill +(18D^370D^244D+288)t+(D4)(13D^2+72D100))tv_{\text{2}\text{h}}^{}\\ \hfill +2(D4)(D3)^2\left(16(D1)t^2+8Dt3D+10\right)v_{\text{2}\text{h}}=0,\end{array}$$
(180)
which is the homogeneous part of the equation satisfied by the function $`v_{\text{2}\text{h}}(t)`$, corresponding to the integral of Fig. 2h of section 9.2, used in the calculation of Eq.(208) (Eq.(180) remains the same for both choices of the line $`D_1`$). The singular points are $`t=1`$, $`1/8`$, $`0`$ and $`(D4)/(4D4)`$. The last singular point, which satisfies the condition (177), is a regular point with exponents $`0`$, $`1`$ and $`3`$. The characteristic equation of the homogeneous part of the corresponding difference equation,

$$\begin{array}{c}(xD+1)(x2D+4)(2x3D+6)(3x4D+10)U_{\text{2}\text{h}}(x1)\hfill \\ \hfill +2(xD+2)(21x^3+(13667D)x^2+(243+64D^2253D)x\\ \hfill 8(D1)(D2)(2D5))U_{\text{2}\text{h}}(x)\\ \hfill 8x(xD+3)(2x3D+7)(3x4D+7)U_{\text{2}\text{h}}(x+1)=0,\end{array}$$
(181)
has only the roots $`1`$ and $`1/8`$.

### 8.6 Applications to simple one-loop integrals

Now we consider the solution with Laplace’s transformation of the difference equations analyzed in the examples of sections 7.1 and 7.2.

#### 8.6.1 One-loop vacuum integral

Considering the integral $`J(x)`$ of Eq.(103), the endpoints of the line of integration are the origin and the root of the characteristic equation $`\mu =1/m_1^2`$, so that we write

$$J(x)=_0^{1/m_1^2}𝑑tt^{x1}v_J(t),$$
(182)
where $`v_J`$ satisfies the differential equation

$$t(m_1^2t1)v_J^{}(t)+(D/2m_1^2t)v_J(t)=0.$$
(183)
The solution is

$$v_J(t)=C(1/m_1^2t)^{D/21}t^{D/2};$$
(184)
the constant $`C`$ can be deduced from the value of $`a_0`$ (see section 108), the value of $`K=D/2`$ and the relation (172). One finds

$$C=(m_1^2)^{D/21}/\mathrm{\Gamma }(D/2).$$
(185)
#### 8.6.2 One-loop self-energy integral

Here we consider $`I(x)`$ of Eq.(111), and in particular the case of non-zero masses. The homogeneous solution can be written using the Laplace’s transformation as

$$I_\pm ^{HO}(x)=_0^{\mu _\pm }𝑑tt^{x1}v_{HO}(t).$$
(186)
The function $`v_{HO}(t)`$ satisfies the differential equation

$$t\mathrm{\Phi }_1(t)v_{HO}^{}(t)+\mathrm{\Phi }_0(t)v_{HO}(t)=0,$$
(187)
where

$$\begin{array}{cc}\hfill \mathrm{\Phi }_1(t)& =R^2(p^2,m_1^2,m_2^2)(t\mu _+)(t\mu _{}),\hfill \\ \hfill \mathrm{\Phi }_0(t)& =t^2R^2(p^2,m_1^2,m_2^2)+(D1)(p^2+m_1^2m_2^2)t+2D.\hfill \end{array}$$
(188)
The solution of this equation is

$$v_{HO}(t)=C\left((\mu _+t)(\mu _{}t)\right)^{(D3)/2}t^{2D}.$$
(189)
Values of $`C`$ such that $`I_\pm ^{HO}`$ are the same functions considered in section 7.2.1 are

$$C_\pm =\mu _\pm ^{(D1)/2}\left(\mu _{}\mu _\pm \right)^{(3D)/2}/\mathrm{\Gamma }((D1)/2).$$
(190)
Considering now the nonhomogeneous solution

$$I^{NH}(x)=_0^{1/m_1^2}𝑑tt^{x1}v_{NH}(t),$$
(191)
the function $`v_{NH}`$ satisfies the differential equation

$$t\mathrm{\Phi }_1(t)v_{NH}^{}(t)+\mathrm{\Phi }_0(t)v_{NH}(t)=t\varphi _1(t)v_J^{}(t)\varphi _0(t)v_J(t),$$
(192)
where $`\mathrm{\Phi }_0`$ and $`\mathrm{\Phi }_1`$ are given in Eq.(188), $`\varphi _0`$ and $`\varphi _1`$ are

$$\begin{array}{cc}\hfill \varphi _1(t)& =t(p^2+m_2^2)/m_1^2,\hfill \\ \hfill \varphi _0(t)& =t(D(p^2+m_1^2+m_2^2)/(2m_1^2)1),\hfill \end{array}$$
(193)
and $`v_J`$ is given in Eq.(184).

#### 8.6.3 Numerical example

We consider the calculation using Laplace’s transformation of $`J(1)`$, $`I_{}^{HO}(1)`$ and $`I^{NH}(1)`$, numerically calculated in section 7.2.4. The corresponding differential equations are Eq.(183), Eq.(187) and Eq.(192). The singular points of the system are $`t=1`$, $`1/3`$ and $`0`$. The equations are integrated using the method described in section 8.5. The integral over $`t`$ is divided into 4 intervals, with endpoints $`0`$, $`1/8`$, $`1/4`$, $`1/2`$, $`1`$; a cutoff $`\lambda 1`$ is conveniently introduced in the first interval $`[\lambda ,1/8]`$ because the solutions are not regular in $`t=0`$. Doing the calculations with 19 digits, the convergence of the expansions in $`t`$ is attained in about 80 terms; the finite values of $`J(x)`$, $`I_{}^{HO}(x)`$ and $`I^{NH}(x)`$ for $`x=3`$ and $`x=4`$ are obtained by calculating the integrals with Eq.(173); the unstable recurrence relations are used only to calculate the values for $`x2`$, where the integrals are divergent, so that the error on $`I(1)`$ is reduced by two orders of magnitude with respect to Eq.(150).

## 9 Application to multi-loop diagrams

After the self-energy diagram discussed in section 7.2, now we consider more complicated diagrams: the vacuum and the self-energy diagrams up to three loop, shown in Figs. 1 and 2, and the vertex and box diagrams up to two loops, shown in Figs. 3 and 4 (we have considered all diagrams such that the scalar integral (18) is always a master integral which does not factorize in a product of simpler master integrals). A complete discussion of all these diagrams would be rather long and we postpone it to future papers. However, to give an idea of the kind and complexities of the equations involved in the calculations, we will show some results.

First of all, we distinct in each diagram $`g`$ the topologically different lines which are indicated in the figures with a number; of course in diagrams without numbers all lines are topologically equivalent. For each topologically different line $`l`$ we set $`D_1`$ equal to the denominator of such line and we consider the difference equation satisfied by the scalar master function $`U_{gl}(x)`$ and the differential equation satisfied by the Laplace-transformed function $`v_{gl}(t)`$, both corresponding to the scalar master integral,

$$U_{gl}(x)=_l𝑑tt^{x1}v_{gl}(t)=\frac{[d^Dk_1]\mathrm{}[d^Dk_{N_k}]}{D_1^xD_2\mathrm{}D_{N_d}}.$$
(194)
### 9.1 Arbitrary case

In the left half of Table 3, considering arbitrary (non exceptional) values of masses and momenta, we list151515 Due to the present limitations of the program used for the calculations, the values listed in this part of the table were calculated by giving arbitrary rational values to square masses and scalar products of external momenta. It is possible, but very unlikely, that the chosen values correspond to some particular case so that the results obtained do not correspond to the real arbitrary case.

* The diagram considered and, if present, the indication of the possible values of the index of the topologically different lines.
* The number $`n_b`$ of master integrals containing all the $`N_d`$ denominators, determined with the procedure of section 2.5, subdivided according to the number of scalar products in the numerator; for example, 1,4,3 means that we have found 8 master integrals altogether, of which 1 with no scalar product, 4 with 1 scalar product and 3 with a product of 2 scalar products.
* In the $`R_C`$ column we list the values of the order $`R`$ of the difference equations in $`x`$ satisfied by $`U_{gl}(x)`$, for each possible choice of $`D_1`$ as one of the topologically distinct lines; the index $`C`$, where present, indicates (assuming values of external momenta below the deformation threshold) the number of constants $`\eta _j`$ which are different from zero because the corresponding partial solutions of the homogeneous equation satisfy the condition (73).
* The order $`S`$ of the differential equations in $`t`$ satisfied by the function $`v_{gl}(t)`$, for each possible choice of $`D_1`$.

The orders shown within parenthesis are estimated from the subsystems of simultaneous equations (26), avoiding the transformation into triangular form; the index $`C`$ is not shown in these cases. Analyzing the data of the table we observe that

1. The differential equation has order less than the order of the difference equation; this is expected as $`v_{gl}(t)`$ is an object simpler than $`U_{gl}(x)`$.
2. The order of the differential equation in $`t`$ is equal to the number of master integrals, with the curious exception of the diagrams 2f3 and 2g1 where $`S<n_b`$.
3. The choice of the line heavily affects the order and the complexity of the equations, as in the case of the diagram 2f, where the difference equation may have order 5,15 or 2.
4. Some similarities appear between vertex diagrams, self-energy diagrams and vacuum diagrams with different number of loops:161616 The equivalence of recurrence relations very recently described in  may probably throw light on this.

1. between a vertex diagram and the self-energy diagrams obtained by connecting two external vertices of the vertex diagram with a line and inserting an external momentum in this new line;
    2. between a self-energy diagram and the vacuum diagram obtained by connecting its two external lines.

The total number of master integrals and the number of master integrals subdivided according to the number of scalar products turn out to be the same, as well as the orders of the equations for the lines present in both diagrams. For example, compare the diagrams (3a, 2d, 1e), (2a, 1b), (2b, 1c), (2c, 1d), (3b, 2j, 2l), (3c, 2o, 2p), (3d, 2r, 2s), (3e, 2t, 2u) and (3f, 2v). The similarities presumably exist even between three-loop self-energy diagrams and four-loop vacuum diagrams (some preliminary results seem to confirm it). Perhaps there is a relation with the heuristic “rule of the mapping” described in  for massless self-energy diagrams.
5. The number of master integrals grows probably exponentially171717We consider here the master integrals containing *one particular* combination of denominators. The total number of combinations (8) itself grows exponentially with $`N_d`$, clearly along with the total number of master integrals with every possible combination of denominators.  with the number of loops and external vertices, and may be large. Consider for example the class of $`L`$-loop “sunset” self-energy diagrams with $`L+1`$ denominators, shown for one, two and three loops respectively in Fig. 2a, 2b and 2e. The number $`n_b(L)`$ of master integrals for $`L=1`$ to $`L=5`$ is $`1,4,11,26,57`$, respectively (the last two values come from preliminary analysis). These values seem to follow the law $`n_b(L)=2^{L+1}L2=_{i=2}^{L+1}\left(\genfrac{}{}{0pt}{}{L+1}{i}\right)`$, corresponding to the (alternative) choice of all the master integrals with numerator equal to one and with one or two as exponents of the denominators, with at least two exponents one.

### 9.2 The test case

As first test of our approach, we have considered in particular detail the case where all masses are equal, $`m_1=\mathrm{}=m_{N_d}=1`$, and all the external lines are on-mass-shell. Denoting the incoming external momenta $`p_1`$, $`p_1p_2`$,…, $`p_{N_p1}p_{N_p}`$, $`p_{N_p}`$, we choose $`p_i^2=1`$ and $`(p_ip_j)^2=1`$ for every $`i`$ and $`j`$. In the case of box diagrams this corresponds to setting the Mandelstam variables $`s=t=1`$ and $`u=2`$. These values of masses and momenta have been chosen because they introduce symmetries which allow several consistency checks on equations and results. Similarly to the arbitrary case, in the right half of Table 3, for each diagram we list the number $`n_b^{}`$ of master integrals with all the denominators, the order $`R^{}`$ of the difference equation satisfied by the function $`U_{gl}(x)`$, the order $`S^{}`$ of the differential equation satisfied by the corresponding function $`v_{gl}(t)`$ and, if present, the number $`C^{}`$ of non-zero constants needed to determine the solution.

For these particular values of masses and momenta the equations turn out to be, as expected, simpler than in the arbitrary case. The homogeneous parts of the equations which presented some similarities in the arbitrary case (see the observation 4 of section 9.1) here turn out to be *identical*; of course the nonhomogeneous parts of these equations are quite different.

The number of constants $`\eta _j`$ to find turns out to be greater than or equal to that of the arbitrary case (except for diagram 2e); in this connection the test case is more complicated than the arbitrary case. With the exception of the diagram 4d6, the calculation of the scalar master integrals requires no more than two constants, easily determined using the identities of section 5.5 (which turn out always to provide one useful relation involving the constants) and the large-$`x`$ leading behaviours. One can prove that the values of external momenta are always below or at the deformation threshold, with the exception of some subdiagrams of diagrams 4f and 4i; in the cases at deformation threshold the values of the constants $`\stackrel{~}{L}_1,\stackrel{~}{L}_2`$ and $`\stackrel{~}{L}_3`$ given in Eq.(94) are needed (for example, for the diagrams 2a, 3a and 4a).

The characteristic equations of all difference equations have always the solution $`\mu =1`$. We list the values of the other roots for some diagrams: (Fig. 1b: $`\mu =1/3`$), (1c: $`1/8`$), (1e: $`1/2`$), (2e: $`1/3`$, $`1/15`$), (2i1: $`1/3`$, $`1/2`$), (2i2: $`1/3`$), (2i3, 2j2, 2l1 and 3b2 : $`1\pm \sqrt{4/3}`$), (2v1-2 and 3f: $`1/3`$, $`1/9`$), (3a: $`1/2`$), (4a: $`3/5`$), (4b1-2: $`1/2`$, $`1/4\pm i\sqrt{1/8}`$), (4b3: $`3/4\pm \sqrt{27/32}`$), (4f1: $`1/2`$, $`1/3`$, $`1`$, $`(1\pm \sqrt{3})/2`$, $`0.098\pm 0.050i`$, $`0.934`$, $`3.632`$), (4f2: $`1/2`$, $`1/3`$, $`1`$, $`(1\pm \sqrt{3})/2`$, $`0.049\pm 0.099i`$, $`2.607`$, $`0.819`$, $`1.756`$, $`12.07`$), (4i1: $`1/9`$, $`3/2`$, $`1/2`$, $`1/3`$, $`1`$) and (4i2: $`1/9`$, $`3/2`$, $`1/3`$, $`3`$). Most of these values, written in the form $`1/(q^2+1)`$, are due to presence of singularities at $`k^2=q^2`$ in the function $`f(k^2)`$ (see section 5.4). All characteristic equations turn out to have one negative rational solution $`1<\mu <0`$, so that, according to section 6.1, all recurrence relations are unstable; the diagram 2e shows the greatest instability ($`A=15`$). One can ask whether the difference equations of all the analyzed diagrams admit as solutions *convergent* factorial series expansions. The answer is negative. In fact the root $`\mu =1/9`$ of the diagrams 2v1-2 and 3f, the root $`(1+\sqrt{3})/2`$ and the complex roots of diagrams 4f1-2, the roots $`1/9`$ and $`3/2`$ of diagrams 4i1-2 (note, all the diagrams with lines crossed) and the root $`3/4+\sqrt{27/32}`$ of the diagram 4b3 satisfy the condition (102) (with $`\mu ^{(\alpha )}=1`$) so that the factorial series expansions of the solutions never converge. Therefore we are forced to use Laplace’s transformation for these diagrams and these choices of the line $`D_1`$, and for all the diagrams which become these diagrams by deleting lines. In all the other cases factorial series expansions can be used.

### 9.3 Numerical results

For each diagram shown in Figs. 1, 2, 3 and 4, except for the diagrams 4f and 4i (see section 9.3.4), we have calculated the values of the master integrals for $`D=42ϵ`$, using the values of masses and momenta shown in the previous section.

Calculations were carried out by using the program SYS described in section 11. The master integrals of the diagrams 1b, 1e, 2a, 2d, 2i, 2j, 2l, 3a, 3e, 4a and 4c were first calculated using expansions in factorial series; these integrals were also recalculated using Laplace’s transformation in order to provide important checks of the calculations. The master integrals of all the remaining diagrams were calculated using Laplace’s transformation. Excluding the simplest diagrams, calculations with Laplace’s transformation turned out faster than calculations with factorial series; this is due to the instabilities of the recurrence relations, which become deeper increasing the number of loops, and which force calculations with factorial series to be performed with a larger number of digits. To give an idea of the size of calculations, the systems of difference equations between the master integrals of the diagrams 1e, 2d, 2t, 2u, 2v, 3f, 4g, and 4h are formed with 44, 28, 245, 304, 362, 81, 139 and 158 equations, respectively; note that in each system, in the right-hand side of the equation corresponding to the more complicated master integral almost all the other master functions appear. We made no use of the symmetries due to the particular values of masses and momenta in order to simplify or reduce the number of the equations, as the aim of program SYS is to deal with calculations of multi-scale integrals, lacking in such symmetries; we used them only to check the final results. In order to guarantee results with at least 20 digits of precision, calculations with factorial series were performed with precision up to 77 digits (depending on the degree of instability), while calculations with Laplace’s transformation were all performed with 38 digits of precision.

For example, the calculation from scratch of the integral $`I(\text{2}\text{t})`$, Eq.(220), with Laplace’s transformation, requested about 128 hours of CPU time on a 133 MHz Pentium PC; 16 hours were used for the determination of the systems of difference and differential equations, obtained by solving systems up to 43000 identities. The solution of the systems (245 equations) yields, as a byproduct, also the values of all simpler master integrals, including Eqs.(205)-(208), Eqs.(210)-(214), Eq.(216) and Eq.(219). We stress that at this preliminary stage of development we directed our efforts to devise tests and cross checks rather than to speed up the program.

For brevity, for each diagram we list here only the values of the master scalar integrals

$$I(diagram)=\frac{[d^Dk_1]\mathrm{}[d^Dk_{N_k}]}{D_1D_2\mathrm{}D_{N_d}}$$
(195)
without scalar products and containing all the $`N_d`$ denominators. As usual, the results have been normalized with the division by $`\mathrm{\Gamma }_ϵ\mathrm{\Gamma }(1+ϵ)`$ raised to the number of loops of the diagram. Coefficients are shown with only 13 digits to save space. Values for $`ϵ=0`$ of all finite integrals have been checked by comparing them with numerical values obtained by performing Monte-Carlo integrations over Feynman parameters, or by performing low dimensional Gaussian integrations on integrands obtained using dispersion relations and hyperspherical variables. As additional consistency check we have repeated the calculation of some diagrams with different choices of the line $`D_1`$, and we have checked that the results obtained are the same.

#### 9.3.1 Vacuum diagrams

$$\begin{array}{c}I(\text{1}\text{a})\mathrm{\Gamma }_ϵ^1=ϵ^11ϵϵ^2ϵ^3ϵ^4+O(ϵ^5),\hfill \end{array}$$
(196)
$$\begin{array}{c}I(\text{1}\text{b})\mathrm{\Gamma }_ϵ^2=1.5ϵ^24.5ϵ^16.98413914196618.00878162355ϵ\hfill \\ \hfill 27.99422356368ϵ^272.00378659799ϵ^3111.9974983355ϵ^4+O(ϵ^5),\end{array}$$
(197)
$$\begin{array}{c}I(\text{1}\text{c})\mathrm{\Gamma }_ϵ^3=2ϵ^3+7.666666666667ϵ^2+17.5ϵ^1+22.91666666667\hfill \\ \hfill +21.25179105129ϵ184.2300051053ϵ^2\\ \hfill 661.1105861534ϵ^33685.054779382ϵ^4+O(ϵ^5),\end{array}$$
(198)
$$\begin{array}{c}I(\text{1}\text{d})\mathrm{\Gamma }_ϵ^3=ϵ^35.666666666667ϵ^215.30161161726ϵ^1\hfill \\ \hfill 46.07511172933148.3508545129ϵ394.1378145809ϵ^2\\ \hfill 1375.669435211ϵ^33466.9749998996ϵ^4+O(ϵ^5),\end{array}$$
(199)
$$\begin{array}{c}I(\text{1}\text{e})\mathrm{\Gamma }_ϵ^3=2.404113806319ϵ^110.03527847977+35.94478903214ϵ\hfill \\ \hfill 119.1503507802ϵ^2+379.7433345095ϵ^31183.320931551ϵ^4+O(ϵ^5).\end{array}$$
(200)
Eq.(197), first six terms of Eq.(198) and first two terms of Eq.(200) agree with the analytical expressions given in ,  and , respectively; remaining terms and other results are new.

#### 9.3.2 Self-energy diagrams

$$\begin{array}{c}I(\text{2}\text{a})\mathrm{\Gamma }_ϵ^1=ϵ^1+0.186200635766+0.021156303568ϵ+0.001726745353ϵ^2\hfill \\ \hfill +0.000109897792ϵ^3+0.000005730593ϵ^4+O(ϵ^5),\end{array}$$
(201)
$$\begin{array}{c}I(\text{2}\text{b})\mathrm{\Gamma }_ϵ^2=1.5ϵ^24.25ϵ^17.37517.22197253479ϵ\hfill \\ \hfill 29.55920705372ϵ^268.87789517038ϵ^3118.2464846454ϵ^4+O(ϵ^5),\end{array}$$
(202)
$$\begin{array}{c}I(\text{2}\text{c})\mathrm{\Gamma }_ϵ^2=0.5ϵ^2+0.6862006357658ϵ^10.6868398873414\hfill \\ \hfill +1.486398391913ϵ2.938796587745ϵ^2+5.871086365958ϵ^3\\ \hfill 11.73616571449ϵ^4+O(ϵ^5),\end{array}$$
(203)
$$\begin{array}{c}I(\text{2}\text{d})\mathrm{\Gamma }_ϵ^2=0.92363182651991.284921671848ϵ+2.689507626490ϵ^2\hfill \\ \hfill 5.338399227511ϵ^3+10.67136736912ϵ^4+O(ϵ^5),\end{array}$$
(204)
$$\begin{array}{c}I(\text{2}\text{e})\mathrm{\Gamma }_ϵ^3=2ϵ^3+7.333333333333ϵ^2+16.02777777778ϵ^1\hfill \\ \hfill +21.92956264368+3.605127475161ϵ184.1413665431ϵ^2\\ \hfill 838.2364324178ϵ^33647.102197031ϵ^4+O(ϵ^5),\end{array}$$
(205)
$$\begin{array}{c}I(\text{2}\text{f})\mathrm{\Gamma }_ϵ^3=ϵ^32.612634286982ϵ^23.906420490690ϵ^1\hfill \\ \hfill +0.5840769314959+1.76460041453ϵ+107.0072031435ϵ^2\\ \hfill +163.2855293783ϵ^3+1372.241466189ϵ^4+O(ϵ^5),\end{array}$$
(206)
$$\begin{array}{c}I(\text{2}\text{g})\mathrm{\Gamma }_ϵ^3=ϵ^35.5ϵ^215.48413914197ϵ^145.68793012675\hfill \\ \hfill 149.1607636537ϵ392.298867227ϵ^2\\ \hfill 1380.125833167ϵ^33455.548194007ϵ^4+O(ϵ^5),\end{array}$$
(207)
$$\begin{array}{c}I(\text{2}\text{h})\mathrm{\Gamma }_ϵ^3=ϵ^35.333333333333ϵ^216ϵ^143.91483126325\hfill \\ \hfill 154.918028663ϵ374.0941853334ϵ^2\\ \hfill 1436.672712535ϵ^33281.940436319ϵ^4+O(ϵ^5),\end{array}$$
(208)
$$\begin{array}{c}I(\text{2}\text{i})\mathrm{\Gamma }_ϵ^3=2.404113806319ϵ^19.763424447585+34.99888165588ϵ\hfill \\ \hfill 116.0420477564ϵ^2+370.0407274069ϵ^31153.646312515ϵ^4+O(ϵ^5),\end{array}$$
(209)
$$\begin{array}{c}I(\text{2}\text{j})\mathrm{\Gamma }_ϵ^3=0.3333333333333ϵ^3+0.5195339690991ϵ^2\hfill \\ \hfill +0.5753609494269ϵ^12.981838135558+12.56596204108ϵ\\ \hfill 44.85302351538ϵ^2+149.1742811721ϵ^3477.1440886129ϵ^4+O(ϵ^5),\end{array}$$
(210)
$$\begin{array}{c}I(\text{2}\text{k})\mathrm{\Gamma }_ϵ^3=0.1666666666667ϵ^3+0.5931003178829ϵ^2\hfill \\ \hfill +0.06234894542402ϵ^11.364667486582+7.062482427894ϵ\\ \hfill 26.87419915573ϵ^2+91.91641284417ϵ^3298.196943613ϵ^4+O(ϵ^5),\end{array}$$
(211)
$$\begin{array}{c}I(\text{2}\text{l})\mathrm{\Gamma }_ϵ^3=0.1666666666667ϵ^3+0.5931003178829ϵ^2\hfill \\ \hfill +0.06234894542402ϵ^11.158877567105+6.268660583427ϵ\\ \hfill 24.1193749759ϵ^2+82.9872059343ϵ^3270.1688760103ϵ^4+O(ϵ^5),\end{array}$$
(212)
$$\begin{array}{c}I(\text{2}\text{m})\mathrm{\Gamma }_ϵ^3=0.3333333333333ϵ^3+0.8528673024324ϵ^2\hfill \\ \hfill 1.728169411584ϵ^1+6.07014140974719.48651365516ϵ\\ \hfill +61.38828756627ϵ^2190.3302695306ϵ^3+583.8045381529ϵ^4+O(ϵ^5),\end{array}$$
(213)
$$\begin{array}{c}I(\text{2}\text{n})\mathrm{\Gamma }_ϵ^3=0.3333333333333ϵ^3+0.8528673024324ϵ^2\hfill \\ \hfill 1.728169411584ϵ^1+6.12035970837519.67063042467ϵ\\ \hfill +62.00178253235ϵ^2192.2586253184ϵ^3+589.7212544716ϵ^4+O(ϵ^5),\end{array}$$
(214)
$$\begin{array}{c}I(\text{2}\text{o})\mathrm{\Gamma }_ϵ^3=0.9236318265199ϵ^12.423491634417+8.381349710069ϵ\hfill \\ \hfill 26.99362121677ϵ^2+85.10096322999ϵ^3263.903318629ϵ^4+O(ϵ^5),\end{array}$$
(215)
$$\begin{array}{c}I(\text{2}\text{p})\mathrm{\Gamma }_ϵ^3=0.9236318265199ϵ^12.116169718457+6.929544685259ϵ\hfill \\ \hfill 21.50327837738ϵ^2+66.32213380401ϵ^3202.887025717ϵ^4+O(ϵ^5),\end{array}$$
(216)
$$\begin{array}{c}I(\text{2}\text{q})\mathrm{\Gamma }_ϵ^3=1.3264482082725.196648136965ϵ+18.37758387804ϵ^2\hfill \\ \hfill 60.41191503661ϵ^3+191.5963941ϵ^4+O(ϵ^5),\end{array}$$
(217)
$$\begin{array}{c}I(\text{2}\text{r})\mathrm{\Gamma }_ϵ^3=1.3413992414475.197752955896ϵ+18.38704656407ϵ^2\hfill \\ \hfill 60.4233521301ϵ^3+191.614009625ϵ^4+O(ϵ^5),\end{array}$$
(218)
$$\begin{array}{c}I(\text{2}\text{s})\mathrm{\Gamma }_ϵ^3=2.0025000411058.162562835907ϵ+29.46716463085ϵ^2\hfill \\ \hfill 98.13080591819ϵ^3+313.871881187ϵ^4+O(ϵ^5),\end{array}$$
(219)
$$\begin{array}{c}I(\text{2}\text{t})\mathrm{\Gamma }_ϵ^3=0.27960892328260.1380294113932ϵ+0.3194688268113ϵ^2\hfill \\ \hfill 0.4399664109267ϵ^3+0.6650515012166ϵ^4+O(ϵ^5),\end{array}$$
(220)
$$\begin{array}{c}I(\text{2}\text{u})\mathrm{\Gamma }_ϵ^3=0.18262723753920.06746690965803ϵ+0.1865462420623ϵ^2\hfill \\ \hfill 0.2498713405447ϵ^3+0.3796187113121ϵ^4+O(ϵ^5),\end{array}$$
(221)
$$\begin{array}{c}I(\text{2}\text{v})\mathrm{\Gamma }_ϵ^3=0.14801330395840.009263002043238ϵ+0.1053308537397ϵ^2\hfill \\ \hfill 0.1224292041846ϵ^3+0.1898480457555ϵ^4+O(ϵ^5).\end{array}$$
(222)
Eq.(202), first four terms of Eq.(203), first term of Eq.(204) and first six terms of Eq.(208) agree with the analytical expressions given in , ,  and , respectively; remaining terms and other results are new.

#### 9.3.3 Vertex diagrams

$$\begin{array}{c}I(\text{3}\text{a})\mathrm{\Gamma }_ϵ^1=0.671253105748+0.1998957762816ϵ+0.03189366853371ϵ^2\hfill \\ \hfill +0.003532937320333ϵ^3+0.0003018185047825ϵ^4+O(ϵ^5),\end{array}$$
(223)
$$\begin{array}{c}I(\text{3}\text{b})\mathrm{\Gamma }_ϵ^2=0.5ϵ^2+0.6862006357658ϵ^10.5916667014024\hfill \\ \hfill +1.356196533114ϵ2.669112118814ϵ^2+5.336651358516ϵ^3\\ \hfill 10.66866283741ϵ^4+O(ϵ^5),\end{array}$$
(224)
$$\begin{array}{c}I(\text{3}\text{c})\mathrm{\Gamma }_ϵ^2=0.671253105748ϵ^10.08774519609257+0.7262375626947ϵ\hfill \\ \hfill 1.32112948587ϵ^2+2.667431469376ϵ^35.332675337091ϵ^4+O(ϵ^5),\end{array}$$
(225)
$$\begin{array}{c}I(\text{3}\text{d})\mathrm{\Gamma }_ϵ^2=0.9371395273151.27184968708ϵ+2.69185047506ϵ^2\hfill \\ \hfill 5.336932134961ϵ^3+10.67100342934ϵ^4+O(ϵ^5),\end{array}$$
(226)
$$\begin{array}{c}I(\text{3}\text{e})\mathrm{\Gamma }_ϵ^2=0.2711563494022+0.1833941077514ϵ+0.05375101058769ϵ^2\hfill \\ \hfill +0.01446103368419ϵ^3+0.000746187372276ϵ^4+O(ϵ^5),\end{array}$$
(227)
$$\begin{array}{c}I(\text{3}\text{f})\mathrm{\Gamma }_ϵ^2=0.173896742268+0.1816664876962ϵ+0.04440899181832ϵ^2\hfill \\ \hfill +0.02231547385785ϵ^30.003079810479797ϵ^4+O(ϵ^5).\end{array}$$
(228)
#### 9.3.4 Box diagrams

$$\begin{array}{c}I(\text{4}\text{a})\mathrm{\Gamma }_ϵ^1=0.3455029252972+0.4731008318818ϵ+0.1519459537543ϵ^2\hfill \\ \hfill +0.0275179284554ϵ^3+0.00348492177519ϵ^4+O(ϵ^5),\end{array}$$
(229)
$$\begin{array}{c}I(\text{4}\text{b})\mathrm{\Gamma }_ϵ^2=0.671253105748ϵ^10.06425178040942+0.7393966927045ϵ\hfill \\ \hfill 1.317187699112ϵ^2+2.668107343399ϵ^35.332500882459ϵ^4+O(ϵ^5),\end{array}$$
(230)
$$\begin{array}{c}I(\text{4}\text{c})\mathrm{\Gamma }_ϵ^2=0.95092356231711.258189955235ϵ+2.694588643167ϵ^2\hfill \\ \hfill 5.335347124508ϵ^3+10.67067120761ϵ^4+O(ϵ^5),\end{array}$$
(231)
$$\begin{array}{c}I(\text{4}\text{d})\mathrm{\Gamma }_ϵ^2=0.276209225359+0.1937422320842ϵ+0.06034849310181ϵ^2\hfill \\ \hfill +0.01640828588681ϵ^3+0.001301642516765ϵ^4+O(ϵ^5),\end{array}$$
(232)
$$\begin{array}{c}I(\text{4}\text{e})\mathrm{\Gamma }_ϵ^2=0.3455029252972ϵ^1+0.4347670080988+0.17885718095363ϵ\hfill \\ \hfill +0.05005382113385ϵ^2+0.006936292250698ϵ^3+0.002511559375421ϵ^4+O(ϵ^5),\end{array}$$
(233)
$$\begin{array}{c}I(\text{4}\text{g})\mathrm{\Gamma }_ϵ^2=0.1723367907503+0.2679578491711ϵ+0.13552112755141ϵ^2\hfill \\ \hfill +0.04468531532833ϵ^3+0.008430602827459ϵ^4+O(ϵ^5),\end{array}$$
(234)
$$\begin{array}{c}I(\text{4}\text{h})\mathrm{\Gamma }_ϵ^2=0.1036407209893+0.2142416932987ϵ+0.14046068671363ϵ^2\hfill \\ \hfill +0.04437197236388ϵ^3+O(ϵ^4).\end{array}$$
(235)
In the case of the diagrams 4f and 4i (from which the diagram 4f is derived with the deletion of a line) the equations (27), obtained after the transformation of the subsystems of equations into triangular form, have orders relatively high ($`10`$) and expressions of size greater than the limit of the computer used, so that our program failed to work out them (the coefficients of the equations would be polynomials in two variables of degree $`100`$); we were able to find only the characteristic and indicial equations. The calculation of these integrals without the transformation into triangular form will be considered in the next paper.

### 9.4 Some examples of equations

Here we discuss the equations satisfied by the functions $`I_g(x)=U_{gl}(x)`$ for some simple diagrams with all lines topologically equivalent, using the values of masses and momenta of section 9.2. For the diagrams 2a and 1b the difference equations are

$`𝕎_\mathrm{𝟏}I_{\text{2}\text{a}}(x)`$ $`={\displaystyle \frac{1}{2}}(D2)I_{\text{1}\text{a}}(x),`$ (236)
$`𝕎_\mathrm{𝟏}I_{\text{1}\text{b}}(x)`$ $`=(D2)I_{\text{1}\text{a}}(1)I_{\text{1}\text{a}}(x),`$ (237)
where $`𝕎_\mathrm{𝟏}`$ is the operator

$$𝕎_\mathrm{𝟏}I(x)=3xI(x+1)+(2xD+1)I(x)+(xD+1)I(x1),$$
(238)
and $`I_{\text{1}\text{a}}(x)`$ is the integral (103) with $`m=1`$. The identity of the homogeneous parts of the equations is an example of the similarities between equations described in section 9.1. The solutions of Eqs.(236)-(237) can be written as

$`I_{\text{2}\text{a}}(x)`$ $`=\eta _{\text{2}\text{a}}I_{\text{2}\text{a}}^{HO}(x)+I_{\text{2}\text{a}}^{NH}(x),`$ (239)
$`I_{\text{1}\text{b}}(x)`$ $`=\eta _{\text{1}\text{b}}I_{\text{2}\text{a}}^{HO}(x)+2I_{\text{1}\text{a}}(1)I_{\text{2}\text{a}}^{NH}(x).`$ (240)
Numerical calculation of $`I_{\text{2}\text{a}}^{HO}(1)`$ and $`I_{\text{2}\text{a}}^{NH}(1)`$ was described in section 7.2.4. The constant $`\eta _{\text{1}\text{b}}`$ may be found by comparing the large-$`x`$ behaviours of the solutions $`I_{\text{2}\text{a}}^{HO}(x)x^{D/2+1/2}`$, $`I_{\text{2}\text{a}}^{NH}(x)x^{D/2}`$ (see Eqs.(126)-(127)) with

$$I_{\text{1}\text{b}}(x)x^{D/2}\frac{[d^Dk_2]}{(k_2^2+1)^2}=x^{D/2}I_{\text{1}\text{a}}(2)$$
(241)
(see Eq.(71)). One finds $`\eta _{\text{1}\text{b}}=0`$. Therefore $`I_{\text{1}\text{b}}(1)=2I_{\text{1}\text{a}}(1)I_{\text{2}\text{a}}^{NH}(1)`$; we see that the two-loop integral $`I_{\text{1}\text{b}}(1)`$ factorizes into a product of a vacuum one-loop integral and *a part* of a one-loop self-energy integral.

Now let us consider the diagrams of Figs. 2b and 1c. Both diagrams have only one master integral with all the denominators, the scalar integral (compared with the four master integrals of the arbitrary case); the corresponding master functions $`I_{\text{2}\text{b}}(x)`$ and $`I_{\text{1}\text{c}}(x)`$ satisfy the difference equations

$`𝕎_\mathrm{𝟐}I_{\text{2}\text{b}}(x)`$ $`=(D2)^2I_{\text{1}\text{a}}(1)I_{\text{1}\text{a}}(x),`$ (242)
$`𝕎_\mathrm{𝟐}I_{\text{1}\text{c}}(x)`$ $`=\frac{3}{2}(D2)^2I_{\text{1}\text{a}}^2(1)I_{\text{1}\text{a}}(x),`$ (243)
where $`𝕎_\mathrm{𝟐}`$ is the operator

$$\begin{array}{c}𝕎_\mathrm{𝟐}I(x)=8x(xD+2)I(x+1)+(7x^2+(1310D)x\hfill \\ \hfill +(3D4)(D1))I(x)+(xD+1)(x3D/2+2)I(x1).\end{array}$$
(244)
Note again the identity of the homogeneous parts of the equations. The characteristic equation has the solutions $`1`$ and $`1/8`$. The index associated to $`\mu =1`$ is $`K=D/2`$; therefore the condition (73) is satisfied and the corresponding solution of the homogeneous equation $`I_{\text{2}\text{b}}^{HO}(x)`$ contributes to $`I_{\text{2}\text{b}}(x)`$ and $`I_{\text{1}\text{c}}(x)`$. The solutions of Eqs.(242)-(243) can be written as

$`I_{\text{2}\text{b}}(x)`$ $`=\eta _{\text{2}\text{b}}I_{\text{2}\text{b}}^{HO}(x)+I_{\text{2}\text{b}}^{NH}(x),`$ (245)
$`I_{\text{1}\text{c}}(x)`$ $`=\eta _{\text{1}\text{c}}I_{\text{2}\text{b}}^{HO}(x)+\frac{3}{2}I_{\text{1}\text{a}}(1)I_{\text{2}\text{b}}^{NH}(x).`$ (246)
The constants may be obtained by comparing the large-$`x`$ behaviours $`I_{\text{2}\text{b}}^{HO}(x)`$$``$ $`x^{D/2}`$ and $`I_{\text{2}\text{b}}^{NH}(x)x^{D/21}`$ with $`I_{\text{2}\text{b}}(x)x^{D/2}I_{\text{2}\text{a}}(1)`$ and $`I_{\text{1}\text{c}}(x)x^{D/2}`$ $`I_{\text{1}\text{b}}(1)`$. One finds $`\eta _{\text{2}\text{b}}=I_{\text{2}\text{a}}(1)`$ and $`\eta _{\text{1}\text{c}}=I_{\text{1}\text{b}}(1)`$. We consider the calculation of $`I_{\text{2}\text{b}}^{HO}(1)`$ and $`I_{\text{2}\text{b}}^{NH}(1)`$ with factorial series. The recurrence relation (244) is unstable with $`A=8`$. We fix a precision of the coefficients of the powers of $`ϵ`$ of the results of $`E=13`$ digits, with $`n_ϵ^{}=7`$ terms of the expansions in $`ϵ`$; following section 6 we choose $`x_{max}=25`$, and we perform the calculations with $`C=38`$ digits. Expansions in factorial series of $`I_{\text{2}\text{b}}^{HO}(x)`$ and $`I_{\text{2}\text{b}}^{NH}(x)`$ converge for $`x=25`$ in about 1800 terms. Values for $`x=1`$ of the solutions are calculated using repeatedly the recurrence relations (242). Two terms of the expansions in $`ϵ`$ are lost going beyond the abscissa of convergence $`\lambda =3`$, so that we must retain the first $`n_e=9`$ terms of the expansions. One obtains

$$\begin{array}{c}\eta _{\text{2}\text{b}}I_{\text{2}\text{b}}^{HO}(1)\mathrm{\Gamma }_ϵ^2=0.09188814923697ϵ^3+0.194632539439ϵ^2\hfill \\ \hfill 0.045490472375ϵ^1+6.50912255436+10.43240978278ϵ\\ \hfill +76.7023111407ϵ^2+118.6149739413ϵ^3+732.0021187015ϵ^4+O(ϵ^5),\end{array}$$
(247)
$$\begin{array}{c}I_{\text{2}\text{b}}^{NH}(1)\mathrm{\Gamma }_ϵ^2=0.09188814923697ϵ^31.694632539439ϵ^2\hfill \\ \hfill 4.204509527625ϵ^113.8841225543627.65438231756ϵ\\ \hfill 106.2615181944ϵ^2187.4928691117ϵ^3850.2486033469ϵ^4+O(ϵ^5).\end{array}$$
(248)
Summing the results one finds Eq.(202); note the cancellation of the $`ϵ^3`$ term, and the partial cancellations of digits in the $`ϵ^2`$, $`ϵ^1`$ and constant terms. Considering now the solution of Eq.(242) with Laplace’s transformation, there are *three* master transformed functions: $`v_{\text{2}\text{b}}(t)`$, $`w_1(t)`$ and $`w_2(t)`$ defined by

$$I_{\text{2}\text{b}}(x)=\frac{[d^Dk_1][d^Dk_2]}{(k_1^2+1)^x(k_2^2+1)((pk_1k_2)^2+1)}=_0^1𝑑tt^{x1}v_{\text{2}\text{b}}(t),$$
(249)
$$K_\alpha (x)=\frac{[d^Dk_1][d^Dk_2](pk_\alpha )}{(k_1^2+1)^x(k_2^2+1)((pk_1k_2)^2+1)}=_0^1𝑑tt^{x1}w_\alpha (t),\alpha =1,2.$$
(250)
Clearly $`w_1(t)w_2(t)`$ and $`K_1(x)K_2(x)`$ as the insertion of $`x`$ breaks the symmetry $`k_1k_2`$. The function $`v_{\text{2}\text{b}}(t)`$ satisfies the differential equation

$$\begin{array}{c}t^2(1t)(1+8t)v_{\text{2}\text{b}}^{\prime \prime }+t(8(1D)t^2+10(D2)t+(5D/26))v_{\text{2}\text{b}}^{}\hfill \\ \hfill +(3D8)(D3)(t+1/2)v_{\text{2}\text{b}}=(D2)^2tI_{\text{2}\text{a}}(1)v_{\text{2}\text{a}},\end{array}$$
(251)
where $`v_{\text{2}\text{a}}`$ is $`v_J`$ of Eq.(184) with $`m=1`$. $`K_1(1)`$ and $`K_2(1)`$ are not master integrals, as $`K_1(1)=K_2(1)=\frac{1}{3}I_{\text{2}\text{b}}(1)`$; therefore $`K_1(x)`$ and $`K_2(x)`$ are not master functions and satisfy difference equations of order *zero*, for example

$$(x+2D)(3K_1(x+1)I_{\text{2}\text{b}}(x+1)+I_{\text{2}\text{b}}(x))=(2D)I_{\text{2}\text{a}}(x+1)I_{\text{2}\text{a}}(1).$$
(252)
On the contrary, $`w_1(t)`$ and $`w_2(t)`$ satisfy *first* order differential equations, for example

$$6t(tw_1^{}+(D2)w_1)+2t(1t)v_{\text{2}\text{b}}^{}+2((2D)t+D3)v_{\text{2}\text{b}}=2(D2)I_{\text{2}\text{a}}(1)tv_{\text{2}\text{a}},$$
(253)
and therefore they are master functions.

## 10 Differential equations in masses and momenta

As shown in the one-loop example, the calculation of the constants $`\eta _j`$, the factors multiplying the solutions of the homogeneous equation, presents different degrees of difficulty according to the values of masses and momenta. For values below the deformation threshold, these factors (when non-zero) are easily expressed in form of integrals with one loop missing, but, above the deformation threshold or when some masses are zero, as discussed in section 5.4, the calculation of the factors becomes more complicated. An alternative to this calculation is to calculate the master integrals for values of masses and momenta such that the calculation of the factors is simple, and subsequently to build and integrate a system of differential equations in the masses and momenta in order to reconstruct the integrals for the requested values of masses and momenta. The approaches to calculation of diagrams based on construction of differential equations in masses or momenta were introduced in ; up to now they have been used in the analysis of particular diagrams (e.g. massless box diagrams in ) with a small number of master integrals by building and solving small systems of differential equations. Here we describe how to build and solve systems of differential equations in masses and momenta in way very similar to that of the differential equations obtained by Laplace’s transformation. In this way even large systems containing hundreds of equations can be built and solved, obtaining results expanded at will in $`ϵ`$.

### 10.1 Differentiation of master integrals

A generic master integral is a function $`B(m_i^2,p_ip_j)`$ of scalar products and masses. We parametrize the trajectory from the initial point to the final point in the space of scalar products and masses with a parameter $`y`$; the scalar products and masses of the points of the trajectory will be in general functions of $`y`$. In analogy with the differential equations obtained by using the Laplace’s transformation, we introduce the parametrization in such a way that the initial point corresponds to $`y=1`$ and the final point corresponds to $`y=0`$. Derivatives with respect to $`y`$ of master integrals will be written as derivatives with respect to scalar products and square masses:

$$\frac{dB}{dy}=\underset{i=1}{\overset{N_p}{}}\frac{dp_i^2}{dy}\frac{B}{p_i^2}+\underset{i<j}{}\frac{d(p_ip_j)}{dy}\frac{B}{(p_ip_j)}+\underset{i=1}{\overset{N_d}{}}\frac{dm_i^2}{dy}\frac{B}{m_i^2}.$$
(254)
Differentiation with respect to masses of the integral is trivial; differentiation with respect to scalar products is less immediate. Let us consider a diagram with $`N_p`$ independent external momenta $`p_i`$; the scalar products $`p_\alpha p_\beta `$ are $`N_{pr}=N_p(N_p+1)/2`$. The derivative with respect to $`p_\alpha p_\beta `$ has the general form

$$\frac{}{(p_\alpha p_\beta )}=\underset{i=1}{\overset{N_p}{}}\underset{j=1}{\overset{N_p}{}}a_{ij\alpha \beta }p_i\frac{}{p_j},$$
(255)
for $`(\alpha ,\beta )`$ equal to each one of the $`N_{pr}`$ possible different pairs of integer. The coefficients $`a_{ij\alpha \beta }`$ can be determined by imposing the $`N_{pr}^{}{}_{}{}^{2}`$ conditions

$$\frac{(p_\gamma p_\delta )}{(p_\alpha p_\beta )}=\{\begin{array}{cc}1\text{ if }p_\gamma p_\delta =p_\alpha p_\beta ,\hfill & \\ 0\text{ otherwise }.\hfill & \end{array}$$
(256)
But the number of the coefficients $`a_{ij\alpha \beta }`$ is $`N_{p}^{}{}_{}{}^{2}N_{pr}`$ which is greater than $`N_{pr}^{}{}_{}{}^{2}`$; therefore some of these coefficients are not constrained by these conditions and may be chosen arbitrarily. We conveniently rewrite Eq.(255) as

$$\frac{}{(p_\alpha p_\beta )}=\underset{i=1}{\overset{N_p}{}}a_{ii\alpha \beta }p_i\frac{}{p_i}+\underset{i<j}{}a_{ij\alpha \beta }\left(p_i\frac{}{p_j}+M_{ij}p_j\frac{}{p_i}\right),$$
(257)
where the coefficients $`a_{ij\alpha \beta }`$ to be found are $`N_{pr}^{}{}_{}{}^{2}`$ and there are $`N_p(N_p1)/2`$ arbitrary constants $`M_{ij}`$. In order to calculate the coefficients it is convenient, fixed $`\alpha `$ and $`\beta `$, to apply Eq.(257) to a generic scalar product $`p_kp_l`$, and to construct the tensor $`A_{ijkl}`$ formed by the factors multiplying $`a_{ij\alpha \beta }`$ in the obtained expression. Now, considering the various pairs of indices $`\{i,j\}`$ as one single index $`ij`$, $`A`$ becomes a square matrix $`A_{ij,kl}`$ of dimension $`N_{pr}`$. Inverting it we find the desired coefficients

$$a_{kl,ij}=\left(A^1\right)_{ij,kl}.$$
(258)
As an example, for a vertex diagram $`N_p=2`$, and the matrix $`A`$ is

$$A=\left(\begin{array}{ccc}2p_1^2& 2M_{12}p_1p_2& 0\\ p_1p_2& p_1^2+M_{12}p_2^2& p_1p_2\\ 0& 2p_1p_2& 2p_2^2\end{array}\right),$$
(259)
where the indices 1,2,3 correspond to $`ij=11,12,22`$.

### 10.2 Construction of the system of differential equations

The algorithm used for the generation and solution of the system of identities is the following:

###### Algorithm 4

Follow the algorithm 1 with the following modifications:

1. New identities which transform the derivatives with respect to $`y`$ of master integrals into combinations of integrals obtained using Eq.(254) are added to the system.
2. Master integrals and their derivatives have a priority of extraction lower than other integrals.
3. Add the new entry “the greatest derivative” to the list of priorities after the entry 9(b)vi.

Analogously to sections 2.5 and 3.2, the previous algorithm allows one to determine the list of the master integrals $`B_n(y)`$, and to obtain a system of differential equations between them. Applying a procedure of transformation into triangular form of the subsystems of equations between master integrals with the same denominators, the whole system of differential equations takes the triangular form. The arbitrary constants $`M_{ij}`$ turn out to be particularly useful as a check, because the expression of derivatives (254) contains $`M_{ij}`$ (through Eq.(257)), while the differential equations obtained from the system of identities are obviously independent of them.

The integration of the system of differential equations requires the knowledge of the values of the master integrals and of some derivatives in the initial point $`y=1`$. The values of the master integrals are obtained by solving a system of difference equations. The derivatives of the master integrals in the initial point are expressed in terms of the master integrals of the system of difference equations using the identities provided by the algorithm 1.

Note that the master integrals of the system of difference equations may be different in number and structure from the master integrals of the system of differential equations. The system of differential equation is solved with the same procedure used for the Laplace’s transformation method. We stress that mobile apparent singular points, as described in section 8.5.2, may appear in the system of differential equations.

### 10.3 Expansion in $`ϵ`$ and infrared divergences

All quantities, coefficients of differential equations, solutions, $`\mathrm{\Gamma }`$ functions, etc. are expanded in $`ϵ=(4D)/2`$. Solution of differential equations using power series in $`y`$ with coefficients expanded in $`ϵ`$ does not cause difficulties, unless we deal with IR divergences. If a master integral $`B(y,ϵ)`$ is IR finite in the initial point $`y=1`$, but is IR divergent in the final point $`y=0`$, it is not possible to obtain the correct dimensionally regularized value, which contains an additional IR pole in $`ϵ`$, by integrating the differential equation from the initial point, where the integral is IR finite: the coefficients of the expansion in $`ϵ`$ of $`B(y,ϵ)`$ tend to infinity as $`y0`$, so that $`lim_{y0}B(y,ϵ)B(0,ϵ)`$. In this unfortunate case the only remedy is to substitute the master integral IR divergent with another integral (or combination of integrals) IR finite. Our preference for master integrals with scalar products in the numerator and all exponents 1 in the denominator rather than integrals with exponents greater than 1 in the denominator was devised to limit the appearance of master integrals IR divergent.

In the case where $`B(0,ϵ)`$ is IR finite, but some derivatives with respect to $`y`$ of $`B(y,ϵ)`$ are IR divergent in $`y=0`$, the problem is overcome by integrating the differential equation up to a small value $`y0`$, suitable chosen so that the value of the master integral will be calculated with the requested precision.

### 10.4 Examples of calculation of integrals with zero masses

For illustrative purpose, we consider here only two simple diagrams; applications to more complicated diagrams will be shown in future papers. Let us consider the diagram of Fig. 5a, and the integral

$$\overline{I}(x)=\frac{[dk]}{(k_1^2+1)^xk_2^2(pk_1k_2)^2},p^2=1,$$
(260)
where $`[dk]=[d^Dk_1][d^Dk_2]`$. It satisfies the homogeneous difference equation

$$x(x2D+5)\overline{I}(x+1)(x3D/2+3)(xD+2)\overline{I}(x)=0;$$
(261)
the solution of this equation is

$$\overline{I}(x)=\eta \frac{\mathrm{\Gamma }(x3D/2+3)\mathrm{\Gamma }(xD+2)}{\mathrm{\Gamma }(x)\mathrm{\Gamma }(x2D+5)}.$$
(262)
The calculation of the constant $`\eta `$ is not simple; in fact, the value of $`p^2=1`$ is above the deformation threshold $`p^2=0`$, so that the large-$`x`$ behaviour of $`\overline{I}(x)`$ contains an additional contribution due to the turning point of the integration path; moreover, the identity (100) is useless because $`\overline{I}(0)=0`$.

We replace the zero mass with a square mass $`y`$ in the denominators

$$D_1=k_1^2+1,D_2=k_2^2+y,D_3=(pk_1k_2)^2+y,$$
(263)
then we build the system of identities following the algorithm 4. The master integrals are $`I_1(y)=[dk]/D_1D_3`$, $`I_2(y)=[dk]/D_1D_2`$, $`I_3(y)=[dk]/D_2D_3`$, $`I_4(y)=[dk]/D_1D_2D_3`$ and $`I_5(y)=[dk](pk_2)/D_1D_2D_3`$, and satisfy the system of differential equations

$$\begin{array}{cc}\hfill yI_1^{}& =(D2)I_1,2yI_2^{}=(D2)I_2,2yI_3^{}=(D2)I_3,\hfill \\ \hfill 8y(y1)I_4^{\prime \prime }& =4\left((4D13)y+72D\right)I_4^{}2(3D8)(D3)I_4+(D2)^2(I_1I_2I_3),\hfill \\ \hfill 6(D2)I_5& =4y(y1)I_4^{}+2\left((D3)y+52D\right)I_4+(D2)(2I_2I_1I_3).\hfill \end{array}$$
(264)
The values of the integrals in the initial point $`y=1`$ can be expressed using the values of integrals with masses equal to one given in section 9.3: $`I_1(1)=I_2(1)=I_3(1)=I^2(\text{2}\text{a})`$, $`I_4(1)=I(\text{2}\text{b})`$, $`I_5(1)=I(\text{2}\text{b})/3`$. The singular points of the system are $`y=1`$ and $`y=0`$. Considering the equation for $`I_4`$, the exponents in these points are

$$\rho _{1,2}^{(y=1)}=0,22ϵ,\rho _{1,2}^{(y=0)}=0,3/22ϵ.$$
(265)
In the initial point the integral $`I_4(1)`$ and its derivatives are IR finite, so that the solution must be regular; therefore the singular solution associated with $`\rho _2^{(y=1)}`$ does not contribute, and the second-order equation for $`I_4`$ needs only $`I_4(1)`$ as initial condition. In the final point the integral $`I_4(0)`$ is IR finite but $`I_4^{}(0)`$ is IR divergent, so that the singular solution associated with $`\rho _2^{(y=0)}`$ must contribute to $`I_4(y)`$.

The integration of the system is divided into two parts: from $`y=1`$ to $`y=1/2`$ expanding in $`y=1`$, and from $`y=1/2`$ to the value $`y=\lambda ^21`$ expanding in $`y=0`$. The presence of the cutoff $`\lambda `$ is needed because the solution is not regular in the origin. Perusing the numerical results, the effect of the cutoff on the coefficient of $`ϵ^s`$ of the expansion in $`ϵ`$ of the final value turns out to be approximately proportional to $`\lambda ^2\mathrm{log}^{m+s}\lambda `$, where $`m`$ is some small integer. A value $`\lambda =10^{15}`$ suffices to obtain the first coefficients with 20 exact digits. The normalized result is

$$\begin{array}{c}I_4(0)\mathrm{\Gamma }_ϵ^2=0.5ϵ^21.25ϵ^14.6648681336964528729\hfill \\ \hfill 9.595397946879509324ϵ26.045799878017610383ϵ^2\\ \hfill 49.501934187562851546ϵ^3120.38235865133474218ϵ^4+O(ϵ^5).\end{array}$$
(266)
Another more complicated example is the integral of Fig. 5b

$$L=\frac{[d^Dk_1][d^Dk_2][d^Dk_3]}{(k_1^2+1)\left((k_1k_2)^2+1\right)\left((k_2k_3)^2+1\right)k_2^2k_3^2(k_1k_3)^2}.$$
(267)
As before, we give a square mass $`y`$ to massless denominators. In this case there are 43 master integrals from 2 to 6 denominators, of which only 1 with 6 denominators. The system of differential equations is too long to be shown here. The singular points are $`y=1,0,1/9,1/4,(3\pm \sqrt{5})/2,1,4,9`$; the points $`1/9`$, $`1/4`$, and $`(3\sqrt{5})/2`$ on the interval $`[0,1]`$ turn out to be regular points of the solutions. The effect of the cutoff is different from the previous case, being proportional to $`\lambda \mathrm{log}^{m+s}\lambda `$. A value $`\lambda =10^{30}`$ suffices to obtain the first coefficients with 20 exact digits. The result is

$$\begin{array}{c}L\mathrm{\Gamma }_ϵ^3=2.4041138063191885708ϵ^13.0270094939876520198\hfill \\ \hfill +22.804522068631748454ϵ53.102275435449702689ϵ^2\\ \hfill +201.84333994219396694ϵ^3577.74024368094326834ϵ^4+O(ϵ^5).\end{array}$$
(268)
The first two terms agree with the results ; subsequent terms are new.

These calculations were carried out by using the program SYS; the calculation from scratch of the integral $`L`$ and all the other master integrals, including the calculation of the master integrals in the initial point, required about 1.5 hours on a 133 MHz Pentium PC.

## 11 The program SYS

Here we report in brief some information concerning the program SYS used to calculate the values of the integrals of the sections 9.2 and 10.4. Further details will be given in a separate paper.

* C program, about 23000 lines, 1Mb executable.
* The program SYS allows one to calculate the value of integrals in almost completely automatic way; only needed input is a description of the diagram, the constants $`a_i`$ and $`b_i`$ of section 2.6 which limit the identities generated, the number of dimensions $`D_0`$ about which to expand the integrals, and the number of terms of the expansions in $`DD_0`$.
* The program contains a simplified algebraic manipulator, used for the solution of the systems of identities.
* Coefficients of the integrals in the identities can be unlimited precision integers, rationals, ratios of polynomials in one and two variables (for example $`D`$ and $`x`$) with integer coefficients. At present the values of square masses and products of external momenta must be rational numbers.
* Efficient management of systems of identities of size up to the limit of disk space (tested up to half million of identities).
* Numerical solution of systems of difference/differential equations up to 500 equations.
* Numerical variables used in the solution are arbitrary precision floating point complex numbers and truncated series in $`ϵ`$ with this kind of coefficients.
* Arithmetic libraries which deal with operations on integers, polynomials, rationals, floating point numbers and truncated series in $`ϵ`$ were written on purpose.

Two versions of the program exist, one using factorial series and one using Laplace’s transformation. The Laplace’s transformation version is much more complicated than the factorial series version but it turns to be faster in many cases: both systems of difference and Laplace-transformed differential equations are generated, then the system of differential equations is solved, and the functions $`U(x)`$ are obtained by integrating over $`t`$ the solutions $`v(t)`$, also checking that $`U(x)`$ are solutions of the system of difference equations. Both versions of program were used for calculating the single-scale integrals of section 9.2.

## 12 Conclusions

Most part of this paper has been devoted to the description of a new method of calculation of master integrals, based on numerical solution of systems of difference equations, obtained by solving systems of identities obtained by integration-by-parts. Important features of the method are the applicability to arbitrary diagrams, inherited from integration-by-parts method, and the ability to obtain very high precision results (even 100-200 digits) expanded at will in $`ϵ`$, due to fact that the calculation of integrals is reduced to sums of *convergent* factorial or power series in one variable.

We also have described two important complements to this method: an algorithm for the reduction of generic Feynman integrals to master integrals, and a procedure for construction and numerical solution of differential equations is masses and momenta; among other things, at present the latter is needed to calculate generic master integrals with massless denominators or with external momenta ‘deep’ in the non-euclidean region (over the deformation thresholds), at least until a working automated general procedure for the calculation of the arbitrary constants of difference equations in these cases will be found.

The implementation of these methods and algorithms in the program $`\mathrm{𝚂𝚈𝚂}`$ turned out to be essential to test and prove the validity of the approach.

In this first exploratory work, mostly devoted to test our approach, we have focused our attention to the calculation of single-scale master integrals, in particular vacuum and self-energy diagrams up to three loops and vertex and box diagrams up to two loops. The calculated values of these master integrals are useful, as they may appear expanding Feynman integrals with respect to ratio of different scale parameters.

But the final targets of our approach are the calculation of four-loop $`g`$-$`2`$ contribution in QED, and the calculation of multi-scale multi-loop master integrals especially in cases where there are no hierarchies between scales, where asymptotic expansions in large masses or momenta seem to be not useful. No insurmountable difficulty seems to exist for applying our approach to these problems. Clearly, that means to deal with a larger number of master integrals, or with more complicated equations, and that will imply modifications or improvements of various parts of the algorithms implemented in the program $`\mathrm{𝚂𝚈𝚂}`$ which, at this stage of development, is far to be optimal. The experience gained by performing the calculations of this work has given many suggestions on the changes which should be made and which will be discussed in future papers.

## Acknowledgement

The author wants to thank E. Remiddi and M. Caffo for useful discussions and encouragement in the very early stage of this work.

.

## Figure Captions

Figure 1: Vacuum diagrams up to three loops.

Figure 2: Self-energy diagrams up to three loops.

Figure 3: Vertex diagrams up to two loops.

Figure 4: Box diagrams up to two loops.

Figure 5: Diagrams with massless lines (dashed). .

(S. Laporta, High-precision calculation of multi-loop Feynman integrals by difference equations)

.

## Condensed paper title

Calculation of Feynman integrals by difference equations .

(S. Laporta, High-precision calculation of multi-loop Feynman integrals by difference equations)

.

.

(S. Laporta, High-precision calculation of multi-loop Feynman integrals by difference equations)

.

.

(S. Laporta, High-precision calculation of multi-loop Feynman integrals by difference equations)