DESY 16-003, DO-TH 16/01A toolbox to solve coupled systems of differential and difference equationsThis work was supported in part by the Austrian Science Fund (FWF) grant SFB F50 (F5009-N15) and the European Commission through contract PITN-GA-2012-316704 (HIGGSTOOLS).

# Desy 16-003, Do-Th 16/01 A toolbox to solve coupled systems of differential and difference equations1

## Abstract

We present algorithms to solve coupled systems of linear differential equations, arising in the calculation of massive Feynman diagrams with local operator insertions at 3-loop order, which do not request special choices of bases. Here we assume that the desired solution has a power series representation and we seek for the coefficients in closed form. In particular, if the coefficients depend on a small parameter (the dimensional parameter), we assume that the coefficients themselves can be expanded in formal Laurent series w.r.t.  and we try to compute the first terms in closed form. More precisely, we have a decision algorithm which solves the following problem: if the terms can be represented by an indefinite nested hypergeometric sum expression (covering as special cases the harmonic sums, cyclotomic sums, generalized harmonic sums or nested binomial sums), then we can calculate them. If the algorithm fails, we obtain a proof that the terms cannot be represented by the class of indefinite nested hypergeometric sum expressions. Internally, this problem is reduced by holonomic closure properties to solving a coupled system of linear difference equations. The underlying method in this setting relies on decoupling algorithms, difference ring algorithms and recurrence solving. We demonstrate by a concrete example how this algorithm can be applied with the new Mathematica package SolveCoupledSystem which is based on the packages Sigma, HarmonicSums and OreSys. In all applications the representation in -space is obtained as an iterated integral representation over general alphabets, generalizing Poincaré iterated integrals.

\xyoption

all \ShortTitleA toolbox to solve coupled systems of differential and difference equations \FullConference12th International Symposium on Radiative Corrections (Radcor 2015) and LoopFest XIV (Radiative Corrections for the LHC and Future Colliders)
15-19 June, 2015
UCLA Department of Physics & Astronomy Los Angeles, USA

## 1 Introduction

A massive three loop Feynman diagram with a local operator insertion can be written in terms of multiple integrals or multiple sums, which depend on a discrete variable and the dimensional parameter , where denotes the space-time dimension, see e.g. see [1, 2]. Then one is interested in the first coefficients of its formal Laurent series w.r.t.  (in short -expansion)

 I(N)=εoIo(N)+εo+1Io+1(N)+εo+2Io+2(N)+… (1)

with order . In this article we present tools to decide algorithmically if the coefficients up to a certain order can be written in terms of indefinite nested hypergeometric sums (in short, nested hypergeometric sums) which can be defined as follows. Let be an expression that evaluates at non-negative integers (from a certain point on) to elements of a field containing the rational numbers . Then is called a nested hypergeometric sum expression w.r.t.  if it is composed by elements from the rational function field , the three operations (), hypergeometric expressions of the form with and being a rational function in and being free of , and sums of the form with and with being a nested hypergeometric sum expression w.r.t.  and being free of . This class of special functions covers as special cases harmonic sums [3],

 Sa1,…,ak(N)=∑N≥i1≥i2≥⋯≥ik≥1sign(a1)i1i|a1|1⋯sign(ak)iki|ak|k (2)

for non-negative integers and nonzero integers and more generally, generalized harmonic sums [4, 5], cyclotomic harmonic sums [6] or nested binomial sums [7, 8].2

In order to calculate the coefficients in (1), the summation package Sigma enhanced by the package EvaluateMultiSums [10, 11], the integration package MultiIntegrate [12], and the package HarmonicSums [13] have been applied successfully in many applications. However, in the course of recent calculations, these tools turned out to be not sufficient and we extended them significantly by using the integration by parts (IBP) identities [14]. Namely, encoding by a generating function (formal power series)

 ^I(x)=∞∑N=0I(N)xN,

we can activate the powerful C++ program Reduze 2 [15] based on Laporta’s algorithm [16] to reduce to a linear combination of master integrals. In many calculations, see e.g., [17, 18, 19] these remaining integrals are now suitable for symbolic summation and integration. However, some of the master integrals are rather hard too handle or are not in the proper form for symbolic summation and integration. Here we utilize the fact that Reduze 2 can produce recursively defined coupled systems of linear differential equations in terms of these master integrals. Following the tactics in [20] the main task is to extract the required information from these coupled systems and to reassemble information of the coefficients in (1) for further processing. In this article we are interested in computing the in closed form. If a first-order coupled system has a specific form, one could use, e.g., the methodology described in [21]. In the following we introduce a very general and efficient approach relying on decoupling algorithms [22, 23, 24], recurrence solvers [25, 26] and difference ring algorithms [27, 28, 29]: we obtain a complete algorithm that extracts the first coefficients of the Laurent series and computes simultaneously the representation of the coefficients in terms of nested hypergeometric sum expressions, whenever this is possible. The first ideas of this new algorithm have been introduced in [30] and the main features are worked out in [19] by concrete examples coming from massive -loop ladder and -diagrams. In the following we will complement these achievements by precise input-output specifications and further details of the algorithms. Moreover, we will illustrate how the differential equation algorithm can be executed within the new package SolveCoupledSystem that relies on the packages Sigma and OreSys. In addition, the package HarmonicSums is used to gain significant speed-ups.

We will use the following notations. Let be a computable field containing the rational numbers as sub-field (e.g., ). In the following (or ) denotes the ring of polynomials in the variable (or in the variables and ). Moreover, (or ) denotes the field of rational functions in the variable (or in the variables and ). We denote by the field of formal Laurent series, i.e., elements are of the form with and . Furthermore, we denote by the ring of power series whose elements are of the form with . Furthermore, (or ) denotes the ring of sequences with entries form (or from ).

## 2 The algorithmic machinery

First, we will address the problem how one can decide algorithmically, if a sequence can be calculated by a nested hypergeometric sum expression provided that the sequence is described by a linear recurrence (linear difference equation) in terms of nested hypergeometric sum expressions (see Theorem 1). Given this technology, we can extract the first coefficients of a Laurent series expansion in terms of nested hypergeometric sum expressions provided that the Laurent series is a solution of a recurrence of certain kind (see Theorem 2). Using uncoupling algorithms, this result can be generalized further to coupled systems (see Theorem 3). Finally, we can carry over this result to coupled systems of linear differential equations (see Theorem 4).

### 2.1 Finding nested hypergeometric solutions of linear recurrences

The main engine relies on the following algorithmic result [11, 29, 25, 26].

Theorem 1. Suppose that a sequence is a solution of the difference equation

with for given rational functions , not all zero, and a nested hypergeometric sum expression . Then one can determine an with the following property.
If one is given the values for all , then one can decide algorithmically if there exists a nested hypergeometric sum expression that calculates the values for all (or at least from a certain point on).

Proof. For this result we refer to Section 4.3 in [11] and Section 2.4 [29]. It has been implemented within Sigma as follows. The recurrence operator is factorized as much as possible into linear factors. Then each linear factor leads to one additional linearly independent solution of the homogeneous version of the recurrence by introducing one extra indefinite summation quantifier and introducing one hypergeometric expression [25]. In this way one finds a basis of all solutions of the homogeneous recurrence that can be expressed in terms of nested hypergeometric sum expressions. If the recurrence factorizes completely, the particular solution can be obtained straightforwardly. However, if the recurrence does not fully factorize, one has to activate algorithms from [26] in order to calculate a particular solution in terms of nested hypergeometric sum expressions or to prove that such a representation is not possible. If there is not such a solution, cannot be represented by a nested hypergeometric sum expression.
In the process of this calculation one can determine a such that the solutions in terms of nested hypergeometric sum expressions can be evaluated and are a solution of the recurrence. Now compute the finite set of all non-negative integer roots of . If , set , else set . Thus for all with we have that . Then using initial values, namely , one can check if the found solutions can be combined to an expression which produces the same initial values. If this is possible, this expression agrees with for all : since the leading coefficient of (3) does not evaluate to zero, there is exactly one sequence which has these initial values and which is a solution of the recurrence. Note that this construction is always successful if one computes linearly independent solutions of the input recurrence (the particular solution is then just a by-product). Otherwise, if this construction fails, it follows that cannot be expressed by a nested hypergeometric sum expression (i.e., that at least one solution of the homogeneous recurrence is of different nature). Summarizing, if we can compute the first initial values of , we can execute the decision procedure described above.

Example. Consider the sequence that is determined by the linear recurrence

In[1]:=

(loaded into Mathematica) and the initial values , , . Loading in the summation package Sigma into Mathematica, one can solve the recurrence in terms of nested hypergeometric sum expressions as follows.

In[2]:=

Sigma - A summation package by Carsten Schneider © RISC-Linz

In[3]:=

Out[3]=

Here the first three entries are linearly independent solutions of the homogeneous version of the recurrence and the last entry is a particular solution of the recurrence itself.
Within Sigma a strong toolbox has been developed to simplify these solutions by flattening the sums optimally and by finding denominators with minimal degrees within the setting of difference rings [10, 11]. In particular, the simplified expressions are built by sums that are algebraically independent [28]. These features can be activated by dropping IndefiniteSummationFalse in In[2.1]. For the class of harmonic sums, generalized harmonic sums, cyclotomic harmonic sums and binomial nested sums these features are also available within the package

In[4]:=

HarmonicSums by Jakob Ablinger © RISC-Linz

For such sums the rather involved difference ring theory can be avoided, and one can calculate very efficiently the simplified representation as follows:

In[5]:=

Out[5]=

Remark. Internally, only the so-called basis-sums remain that cannot be eliminated by relations induced by the underlying quasi-shuffle algebra. For harmonic sums these ideas are worked out in [31] and have been extended for cyclotomic sums, generalized harmonic sums and nested binomial sums [13, 5, 6, 8]. We remark further that the basis sums produce sequences which are algebraically independent [32].

Summarizing, the solution set of Out[2.1] is given by the set

 {c1−1N+1+c2−NN+1+c3(1(N+1)2+S1(N)N+1)+2(N2+N−1)3(N+1)2−2(N+2)3(N+1)S1(N)|c1,c2,c3∈K}

of nested hypergeometric sum expressions. The initial values can be fulfilled with which yields

 I(N)=59N2+120N+499(N+1)2−2(N+3)S1(N)3(N+1). (3)

Hence we have shown that can be calculated for by a nested hypergeometric expression.

### 2.2 Finding Laurent series solutions of linear difference equations

During the calculation of a Feynman integral one often obtains linear recurrences for depending on a dimensional parameter with . In lucky situations one finds a nested hypergeometric sum representation of using the recurrence solver of Subsection 2.1 with . However, in most cases one fails to find any solution of the given recurrence in terms of nested hypergeometric sum expressions, but one finds an -expansion whose coefficients can be represented by nested hypergeometric expressions using the following algorithmic machinery [1].

Theorem 2. Suppose that the sequence with (1) is a solution of the difference equation

for explicitly given and for a sequence with

 r(x)=εoro(N)+εo+1ro+1(N)+εo+2ro+2(N)+… (5)

Here we assume3 that the do not introduce poles for and that not all are zero. Then for any one can determine an with the following property.
If one is given the values for all and and one is given for all nested hypergeometric sum expressions that calculate the values for all (or at least from a certain point on), then one can decide algorithmically if for all there exist nested hypergeometric sum expressions that calculate the values for all (or at least from a certain point on).

Proof. Let and suppose that the can be represented in terms of nested hypergeometric sum expressions. We make the Ansatz (1) with unknown coefficients and plug them into (4). Then the left and right hand sides are both Laurent series which are equal if and only if the coefficients agree. In particular, the lowest term must agree, i.e., we obtain the following constraint

which is a linear recurrence of order . Activating Theorem 1 to this recurrence with appropriately chosen initial values (given by Theorem 1) one can decide algorithmically if is expressible by a nested hypergeometric sum expression. If such a representation is not possible, the theorem is proven. Otherwise, we take this representation and make the Ansatz (1) with the known coefficient and the unknown coefficients () and plug them into (4). Then by construction the coefficients of on the left and right hand sides agree and the term of can be eliminated by subtracting it on both sides. Now one repeats this process for the next lowest term. In this way one can decide algorithmically if all can be represented by nested hypergeometric sum expressions. In the process of this construction we choose such that the used initial values are covered by with .

Example. Take with the -expansion (1) of order where the first two coefficients are determined by the initial values

 I(1)=5ε3−16312ε2+O(ε−1),I(2)=13027ε3−69554ε2+O(ε−1),I(3)=16936ε3−39532ε2+O(ε−1) (7)

and the linear recurrence

In[6]:=

We seek to calculate a nested hypergeometric sum representation for , and , i.e., we set . Note that the expansion on the right hand side of In[2.2] is sufficiently high expanded. The recurrence (6) in our concrete instance is precisely In[2.1] with . In addition, the initial values agree with (7). Hence the found nested hypergeometric sum expression (3) represents . Continuing this process we can calculate the coefficient . Within Sigma this calculation can be carried out automatically with the function call

In[7]:=

Out[7]=

Remark. For further details and speed-ups we refer to [1].

### 2.3 Finding Laurent series solutions of coupled systems of linear difference equations

We generalize Theorem 2 to solve coupled systems as follows [19].

Theorem 3. Suppose that the sequences with

 Ii(N)=εoIi,o(N)+εo+1Ii,o+1(N)+εo+2Ii,o+2(N)+…

are solutions of the coupled system of difference equations

 A0⎛⎜ ⎜⎝I1(N)⋮In(N)⎞⎟ ⎟⎠+A1⎛⎜ ⎜⎝I1(N+1)⋮In(N+1)⎞⎟ ⎟⎠⋯+Ad⎛⎜ ⎜⎝I1(N+d)⋮In(N+d)⎞⎟ ⎟⎠=⎛⎜ ⎜⎝r1(N)⋮rn(N)⎞⎟ ⎟⎠ (8)

for explicitly given matrices with entries from and for some sequences with

 ri(N)=εoro,i(N)+εo+1ro+1,i(N)+εo+2ro+2,i(N)+… (9)

Then for any one can determine and with the following property.
If one is given the values for all , and and one is given for all and nested hypergeometric sum expressions that calculate the values for all (or at least from a certain point on), then one can decide algorithmically if for all and there are nested hypergeometric sum expressions that calculate the values for all (or at least from a certain point on).

Proof. The algorithmic steps can be summarized as follows. First, we transform the coupled system (8) to a first order system as explained in [19]. Then we can apply any decoupling algorithm from [22] to uncouple the system, e.g., w.r.t. . In our implementation we chose Zürcher’s algorithm [23] implemented in the package OreSys [24]. In the generic case one obtains one linear recurrence in which is of the form (5). We can assume that the evaluations are possible and that not all are zero (see Footnote 2). In addition, the decoupling algorithm expresses the remaining sequences by a linear combination of the shifted versions of and shifted versions of the (9).
A subtle point is to which order the expressions in (9) should be expanded. To extract this knowledge, we decouple the system by considering first as unspecified sequences. Then analysing the corresponding output gives an upper bound for the ; details on these aspects can be found in [19]. Since the description of the is given in terms of a linear combination of the shifted versions of , it might be necessary to expand higher than . E.g., if with is one of the components. Analysing these combinations yields the required order . Now we are ready to apply Theorem 2 to calculate the values with by means of a linear recurrence. In this process we determine that the first initial values are needed4. If one fails to get the representation of in terms of nested hypergeometric sum expressions, the theorem is proven. Otherwise, by the properly chosen , the and , this yields also a nested hypergeometric sum representation of the .
In the degenerated case, the decoupling algorithm provides several scalar linear recurrences, say in the , and the remaining are expressed by a linear combination of the shifted versions of the and the . Applying Theorem 2 with the corresponding (as described for the generic case)  times leads to the desired result.

Example. Consider the sequences which are solutions of the coupled system (8) with and where

 A0=(N+100ε(3ε+2)−2(3ε+1)−2(−1+ε−2N)−ε(3ε+2)2(3+3ε+2N)2(ε+1)), A1=(−2−ε−N20−2ε(3ε+2)2(5ε+2)4(−1+ε−N)0−2(4+ε+2N)0)

and

 r1(N)=(−4(N+3)3(N+2)ε−3+(236N3+29N2+45N+21(N+1)(N+2)2−2(2N+3)S1(N)3(N+2))ε−2+O(ε−1),r2(N)=−83ε−3+(4(3N+1)3(N+1)−8S1(N)3)ε−2+O(ε−1),r3(N)=83ε−3+(−4(3N+1)3(N+1)+8S1(N)3)ε−2+O(ε−1)). (10)

In particular, we are given the initial values (7) for We want to derive the -expansions for the up to the order for . Uncoupling the system (first with generic right hand sides) shows that , i.e., the -expansions in (10) are sufficiently high expanded. In particular, we obtain the linear recurrence In[2.2] with and can express and by the shifted versions of and . Summarizing, the coefficients and are computed in Out[2.2]. This yields the needed information to calculate the -expansions of and up to .

In[8]:=

OreSys by Stefan Gerhold (optimized by C. Schneider) © RISC-Linz

In[9]:=

SolveCoupledSystem by Carsten Schneider © RISC-Linz

First, we execute the following command from the package SolveCoupledSystem:

In[10]:=

Out[10]=

This means that one can solve the system by providing the three consecutive initial values of up to order (the starting point depends usually on the physical problem) and that one needs the -expansions of up to the orders . Providing the required information, we execute

In[11]:=

In[12]:=

Out[12]=

and obtain the -expansions of and up to the orders and , respectively.

### 2.4 Finding power series solutions of coupled systems of linear differential equations

Finally, we are ready to present our differential equation solver for coupled systems.

Theorem 4. Suppose that the power series with

 ^Ii(x)=∞∑N=0(Ii,o(N)εo+Ii,o+1(N)εo+1+Ii,o+2(N)εo+2+Ii,o+3(N)εo+3+…)xN (11)

for some common are solutions of the coupled system of differential equations5

 (12)

for explicitly given matrices with entries from and for some with

 ^ri(x)=∞∑N=0(ri,o(N)εo+ri,o+1(N)εo+1+ri,o+2(N)εo+2+ri,o+3(N)εo+3+…)xN. (13)

Then for any one can determine and with the property as stated in Theorem 3.

Proof. This result follows straightforwardly by holonomic closure properties. Namely, take the partial linear differential equations in the and clear denominators by multiplying them with an appropriate polynomial from . Take one of the terms of the equations, say with , and . With the Ansatz (11) we get

 a(ε,x)Dk^Ii(x)= a(ε,x)Dk∞∑N=0(Ii,o(N)εo+Ii,o+1(N)εo+1+Ii,o+2(N)εo+2+…)xN = a(ε,x)∞∑N=0(Ii,o(N)εo+Ii,o+1(N)εo+1+Ii,o+2(N)εo+2+…)k−1∏j=0(N−j)xN−k.

Now we plug all these terms and (13) into the equations. Doing coefficient comparison w.r.t.  leads to a coupled system of linear difference equations in the and the . Performing an appropriate shift in yields a coupled system of the form (8) with the right hand sides which depend linearly on the and their shifted versions from (13). Note that the coefficients in can be expressed in terms of nested hypergeometric sum expressions since the coefficients in can be expressed in terms of nested hypergeometric sum expressions. Thus we can apply Theorem 3 and obtain the claimed result.

Example. Consider the integrals given by , and where

 J6(ν2,ν4;x)=∫dDk1(2π)DdDk2(2π)DdDk3(2π)D1Pν22Pν44P5P7P8P10 (14)

with the propagators , , , , and . Note that can be given in the power series representation (11) and the main task is to determine the coefficients . The IBP algorithm of Reduze 2 [15] delivers e.g. the coupled system (12) with and where is the identity matrix and

 A0=−⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1+ε−x(x−1)x−2(x−1)x0−ε(3ε+2)(x−2)4(x−1)x−2−5ε+x+3εx2(x−1)x(−2ε−x+εx)2(x−1)xε(3ε+2)4(x−1)2+ε−3x−3εx2(x−1)x−ε+12(x−1)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

Furthermore, the are given by a linear combination of master integrals where each one has a power series representation of the form (11). In particular, one can use symbolic summation tools [10, 11] to calculate the coefficients of the coefficients (up to a certain modest order in ) in terms of nested hypergeometric sum expressions. For other situations multiple integration methods [12] are the appropriate tool. Thus the have a power series representation of the form (13) where the can be given (up to a certain modest -order) explicitly in terms of nested hypergeometric sum expressions.
Now we activate our machinery. By holonomic closure properties we obtain precisely the coupled difference system from Example 2.3. Taking the data from Out[2.3] we calculate the required by means of symbolic summation. Furthermore, we calculate the initial values (7) by exploiting the -parametrization of the integrals; for further details on this method we refer to [17]. Note that in other situations we also used our summation tools, provided a reasonable sum representation has been derived. Finally, we activate the function call In[2.3] to get the final result.

Remark. The integrals , and with (14) themselves are master integrals produced by Reduze 2 [15] in order to calculate the diagram6

 (15)

in [19]. There we calculated the -expansion up to order (and not just to order ) in terms of 40 harmonic sums up to weight 7. The total calculation time was 229 seconds. The most complicated coupled system for diagram (15) had dimension . Interestingly enough, the right hand sides of (8) can be given in terms of generalized harmonic sums only, but the solution is given in terms of nested binomial sums. This strongly indicates that straightforward tactics, like transforming the system to a particular shape and reading off the solutions, are not sufficient for such systems.
Note that the obtained results in -space presented in [19] and also in [18, 17] can be transformed to -space by using iterated integral representations over general alphabets, generalizing in parts Poincaré iterated integrals [8].

## 3 Conclusion

We worked out a solver for coupled systems of linear differential equations where in one stroke

1. the coefficients of the formal power series solution are expanded in its