Fuchsia: a tool for reducing differential equations for Feynman master integrals to epsilon form

TL;DR

Fuchsia is a software tool implementing the Lee algorithm to transform differential equations for Feynman master integrals into epsilon form, simplifying their solutions in terms of polylogarithms.

Contribution

It provides an automated implementation of the Lee algorithm specifically tailored for reducing Feynman integral differential equations to epsilon form.

Findings

Successfully reduces differential equations to epsilon form

Ensures solutions contain only regular singularities

Facilitates solving Feynman integrals analytically

Abstract

We present $\text{Fuchsia}$ $-$ an implementation of the Lee algorithm, which for a given system of ordinary differential equations with rational coefficients $\partial_x\,\mathbf{f}(x,\epsilon) = \mathbb{A}(x,\epsilon)\,\mathbf{f}(x,\epsilon)$ finds a basis transformation $\mathbb{T}(x,\epsilon)$, i.e., $\mathbf{f}(x,\epsilon) = \mathbb{T}(x,\epsilon)\,\mathbf{g}(x,\epsilon)$, such that the system turns into the epsilon form: $\partial_x\, \mathbf{g}(x,\epsilon) = \epsilon\,\mathbb{S}(x)\,\mathbf{g}(x,\epsilon)$, where $\mathbb{S}(x)$ is a Fuchsian matrix. A system of this form can be trivially solved in terms of polylogarithms as a Laurent series in the dimensional regulator $\epsilon$. That makes the construction of the transformation $\mathbb{T}(x,\epsilon)$ crucial for obtaining solutions of the initial equations. In principle, $\text{Fuchsia}$ can deal with any regular systems, however its primary task is to reduce differential equations for Feynman master integrals. It ensures that solutions contain only regular singularities due to the properties of Feynman integrals.

Equations

$${\partial }_x\mathbf{J}(x,ϵ)=\mathbb{A}(x,ϵ)\mathbf{J},$$

$$\mathbb{A}(x\to x_k,ϵ)=\{\begin{array}{ll}({\mathbb{A}}_{k0}(ϵ)+{\mathbb{A}}_{k1}(ϵ)(x-x_k)+\dots )/(x-x_k{)}^{1+p_k}& \text{if xk=∞},\\ -({\mathbb{A}}_{k0}(ϵ)+{\mathbb{A}}_{k1}(ϵ)x^{-1}+\dots )x^{-1+p_k}& \text{if xk=∞},\end{array}$$

$$\mathbf{J}=\mathbb{T}(x,ϵ){\mathbf{J}}^{\prime }.$$

$${\partial }_x{\mathbf{J}}^{\prime }={\mathbb{A}}^{\prime }(x,ϵ){\mathbf{J}}^{\prime },$$

$${\mathbb{A}}^{\prime }(x,ϵ)={\mathbb{T}}^{-1}\left(\mathbb{AT}-{\partial }_x\mathbb{T}\right).$$

$$\mathcal{B}(\mathbb{P}(ϵ),x_1,x_2;x)=\mathbb{I}-\mathbb{P}(ϵ)+c\frac{x-x_1}{x-x_2}\mathbb{P}(ϵ),$$

$$c\equiv \{\begin{array}{ll}1/x_1& \text{if x1=∞},\\ x_2& \text{if x2=∞},\\ 1& \text{otherwise},\end{array}$$

$$m_k(\mathbb{A})=p_k+\frac{rank({\mathbb{A}}_{k0})}{n}$$

$${\mu }_k(\mathbb{A})=\min m_k({\mathbb{A}}^{\prime }),\text{ for }\forall \mathbb{T}.$$

$${\Delta }^{r_k}\det \left(\frac{{\mathbb{A}}_{k0}}{\Delta }+{\mathbb{A}}_{k1}-\lambda \mathbb{I}\right),$$

$${\mathbb{A}}^{\prime }(x,ϵ)\equiv ϵ\sum _i\frac{{\mathbb{S}}_i}{x-x_i}=\sum _i{\mathbb{T}}^{-1}(ϵ)\frac{{\mathbb{A}}_{i0}(ϵ)}{x-x_i}\mathbb{T}(ϵ),$$

$${\mathbb{S}}_i={\mathbb{T}}^{-1}(ϵ)\frac{{\mathbb{A}}_{i0}(ϵ)}{ϵ}\mathbb{T}(ϵ).$$

$${\mathbb{S}}_i={\mathbb{T}}^{-1}(ϵ)\frac{{\mathbb{A}}_{i0}(ϵ)}{ϵ}\mathbb{T}(ϵ)={\mathbb{T}}^{-1}(\mu )\frac{{\mathbb{A}}_{i0}(\mu )}{ϵ}\mathbb{T}(\mu ),$$

$$\frac{{\mathbb{A}}_{i0}(ϵ)}{ϵ}\mathbb{T}(ϵ,\mu )=\mathbb{T}(ϵ,\mu )\frac{{\mathbb{A}}_{i0}(\mu )}{\mu },\text{for all i.}$$

$$\mathbb{A}(x,ϵ)=\left(\begin{array}{cccc}{\mathbb{A}}^{(11)}& 0& 0& 0\\ {\mathbb{A}}^{(21)}& {\mathbb{A}}^{(22)}& 0& 0\\ ⋮& \cdots & \ddots & 0\\ {\mathbb{A}}^{(m1)}& {\mathbb{A}}^{(m2)}& \cdots & {\mathbb{A}}^{(mm)}\end{array}\right),$$

$$\det (\mathbb{A}(x,ϵ)-\lambda \mathbb{I})=\det ({\mathbb{A}}^{(11)}(x,ϵ)-\lambda \mathbb{I})\cdot \dots \cdot \det ({\mathbb{A}}^{(mm)}(x,ϵ)-\lambda \mathbb{I}),$$

$$\{\begin{array}{l}{\partial }_x{\mathbf{J}}^{(a)}={\mathbb{A}}^{(aa)}{\mathbf{J}}^{(a)}+\dots \\ {\partial }_x{\mathbf{J}}^{(b)}={\mathbb{A}}^{(ba)}{\mathbf{J}}^{(a)}+{\mathbb{A}}^{(bb)}{\mathbf{J}}^{(b)}+\dots \end{array}$$

$${\mathbf{J}}^{(b)}={\mathbf{J}}^{\prime (b)}+{\Delta }^{-r}\mathbb{D}{\mathbf{J}}^{(a)},$$

$${\mathbb{A}}_{k0}^{\prime (ba)}(ϵ)={\mathbb{A}}_{k0}^{(ba)}(ϵ)+r\mathbb{D}+{\mathbb{A}}_{k0}^{(bb)}(ϵ)\mathbb{D}-\mathbb{D}{\mathbb{A}}_{k0}^{(aa)}(ϵ)$$

Full text

In principle, Fuchsia can deal with any regular systems, however its primary task is to reduce differential equations for Feynman master integrals. It ensures that solutions contain only regular singularities due to the properties of Feynman integrals.

PROGRAM SUMMARY

Program Title: Fuchsia

Authors: O. Gituliar and V. Magerya

Program obtainable from: https://github.com/gituliar/fuchsia/

Journal Reference:

Catalog identifier:

Licensing provisions: ISC license

Programming language: Python 2.7

Operating system: Linux, Unix-like

RAM: Dependent upon the input data. Expect hundreds of megabytes.

Keywords: Computer algebra, Feynman integrals, differential equations, epsilon form, Fuchsian form, Moser reduction

Classification: 5 Computer Algebra, 11.1 High Energy Physics and Computing

External routines/libraries: SageMath (7.0 or higher), Maple (optional)

Nature of problem: Feynman master integrals may be calculated from solutions of a linear system of differential equations with rational coefficients. Such a system can be easily solved as an $\epsilon$-series when its epsilon form is known. Hence, a tool which is able to find the epsilon form transformations can be used to evaluate Feynman master integrals.

Solution method: Methods of Moser [Mos59] and Lee [Lee15].

Restrictions: Systems of single-variable differential equations are considered, however unknown functions may also depend on other symbolic arguments. A system needs to be reducible to Fuchsian form and eigenvalues of its residues must be of the form $n+m\,\epsilon$, where $n$ is integer.

Running time: Depends upon the input data, its size and complexity. Around an hour in total for an example $74\times 74$ matrix with 20 singular points on a PC with a 1.7GHz Intel Core i5 CPU.

1. 1 Introduction

2. 2 Overview of the Lee method

3. 2.1 Notation and definitions

4. 2.2 Reduction to epsilon form

5. 3 Using Fuchsia

6. 3.1 Installation

7. 3.2 Usage from the command line

8. 3.3 Usage from SageMath or Python

9. 4 Summary

1 Introduction

More than 60 years have passed since Richard Feynman proposed a diagrammatic approach for calculating perturbative processes in quantum field theories. Since then Feynman integrals calculus has grown to a separate branch of mathematical physics with a big community of scientists conducting research in this exciting field. With no doubt we can say that none of the recent discoveries in the high-energy particle physics could happen without precise theoretical calculations, which are based on the Feynman integrals calculation techniques. It is also clear that such techniques will play a key role for discoveries at the present and future high-energy colliders, hence their development and further improvement are a very important task.

Recent progress in computational techniques has made it possible to automate the calculation of Feynman integrals; problems which seemed impossible 10 years ago now are successfully solved with state-of-the-art computer algorithms. Among the most popular are integration-by-parts (IBP) reduction [CT81] and the method of differential equations (DE) [Kot91c, Kot91b, Kot91a]; for a detailed overview of these and other methods see [Smi06].

In this paper we discuss the method of differential equations. In particular, we focus on the fact that a solution to the system of DEs may be easily found as an $\epsilon$-series when an epsilon form of this system is known [Hen13]. We consider a general algorithm to find an epsilon form of a given system of differential equations in one variable developed by Lee [Lee15], for the review of the subject see [Hen15] and [Pap14, Tan15, ABB*+*16]. The Lee algorithm allows

To find a Fuchsian form of the system using improved Moser reduction algorithm [Mos59]; and then 2. 2.

To normalize eigenvalues of the Fuchsian system in all singular points.

If these two steps are successfully completed111In principle, the first step can always be done because Feynman integrals contain only logarithmic singularities. For the discussion of potential complications in the second step see Section 2. then a transformation which puts the initial system into the epsilon form may be easily found. Although this method is focused on systems for Feynman integrals, it also may be successfully used for different problems provided that the requirements on particular properties of the initial system are satisfied. For another method for finding the epsilon form of a given system with multiple scales see [Mey16a, Mey16b].

Until recently, no implementation of the Lee method was made publicly available. This fact motivated us to develop Fuchsia — the first public implementation of the Lee algorithm [Lee15] which was presented in [GM16]. Another implementation of this method, called Epsilon, was recently presented in [Pra17]. In this paper we provide a description of some implementation and usage details of Fuchsia which together with algorithms and tools for the integration-by-parts reduction [Lap00, Smi08, vMS12, Lee12, Lee14, SS13, Smi15, GLZ16] form a powerful tandem for evaluating Feynman integrals.

This paper is organized as follows: in Section 2 we introduce notation and definitions followed by a brief review of the Lee method and related algorithms implemented in Fuchsia. In Section 3 we describe how to install and how to use Fuchsia from different environments, depending on your goal and programming experience.

2 Overview of the Lee method

2.1 Notation and definitions

Let us consider a system of ordinary differential equations (ODEs) of this form:

[TABLE]

where $\mathbf{J}(x,\epsilon)$ is a column-vector of $n$ unknown functions (e.g., master integrals); $x$ is a free variable; $\epsilon$ is an infinitesimally small parameter (e.g., a dimensional regulator in $d=4-2\epsilon$ dimensions); $\mathbb{A}(x,\epsilon)$ is an $n\times n$ matrix, rational in both $x$ and $\epsilon$.

In the general case $\mathbb{A}(x,\epsilon)$ may have a finite number of poles in $x$ at $x\in\{x_{k}\}$, including a pole at infinity. The asymptotic behavior of $\mathbb{A}(x,\epsilon)$ around these poles can be described as:

[TABLE]

where $p_{k}$ is the Poincaré rank of $\mathbb{A}(x,\epsilon)$ at the singular point $x=x_{k}$. If $p_{k}=0$, we call $\mathbb{A}(x,\epsilon)$ Fuchsian in $x=x_{k}$, and $\mathbb{A}_{k0}$ matrix residue of $\mathbb{A}(x,\epsilon)$ at $x=x_{k}$. If all $p_{k}=0$, we call $\mathbb{A}(x,\epsilon)$ Fuchsian.

The behavior of the system at $x=\infty$ is a bit of a special case, but it is essential in the overall reduction process. One must always keep in mind the $x=\infty$ point, and treat it on the same footing as other singular points.

Equivalent systems.

We are interested in transforming system (1) into a simpler form. For this purpose let us consider a change of basis from $\mathbf{J}$ into $\mathbf{J}^{\prime}$ using the linear transformation $\mathbb{T}(x,\epsilon)$:

[TABLE]

This leads to the equivalent system of ODEs

[TABLE]

with the new matrix being

[TABLE]

Generally speaking, the transformations $\mathbb{T}(x,\epsilon)$ may have an arbitrary form. However in the scope of this paper we will require the transformation matrices to be rational in both $x$ and $\epsilon$. This restriction guarantees that the equivalent matrix $\mathbb{A}^{\prime}(x,\epsilon)$ and hence all equivalent systems are in rational form, thus making the expansion (2) possible.

In particular, we will be using the transformation constructed by stepwise application of a $\mathbb{P}$-balance between $x=x_{1}$ and $x=x_{2}$, defined as:

[TABLE]

[TABLE]

where $\mathbb{P}(\epsilon)$ is a projector matrix (that is, $\mathbb{P}^{2}=\mathbb{P}$).

Classification of singularities.

Following [Mos59], for a system (1) and it’s Laurent expansion (2) we define a rational number

[TABLE]

as the Moser order of $\mathbb{A}(x,\epsilon)$ at point $x=x_{k}$.

Equivalent systems do not necessarily have identical Moser orders. In fact [Mos59] introduces an algorithm that constructs a transformation decreasing $rank(\mathbb{A}_{k0})$ by at least one (thus reducing $m_{k}(\mathbb{A})$),222Note, that if $rank(\mathbb{A}_{k0})$ reaches zero, this means that $\mathbb{A}_{k0}$ is now zero itself, and thus the Poincaré rank $p_{k}$ was decreased by at least one. or certifies that no further order reduction can be achieved.

Let us then denote the minimal order of $\mathbb{A}(x,\epsilon)$ at $x=x_{k}$ as:

[TABLE]

If $\mu_{k}(\mathbb{A})<m_{k}(\mathbb{A})$ we say that the matrix $\mathbb{A}(x,\epsilon)$ is Moser-reducible at $x=x_{k}$.

With this in mind we can classify a point $x=x_{k}$ of $\mathbb{A}(x,\epsilon)$ as:

• •

regular point, if $m_{k}(\mathbb{A})=0$;

• •

apparent singularity, if $m_{k}(\mathbb{A})>0$ and $\mu_{k}(\mathbb{A})=0$;

• •

regular singularity, if $0<\mu_{k}(\mathbb{A})\leq 1$;

• •

irregular singularity, if $\mu_{k}(\mathbb{A})>1$.

The matrix $\mathbb{A}(x,\epsilon)$ is called Fuchsian if it does not contain irregular singularities at any value of $x$ including $\infty$.

2.2 Reduction to epsilon form

In the previous section we have introduced the notation and key definitions related to the Fuchsian theory of ODEs. Now we are ready to review the reduction method proposed by Lee in [Lee15]. With its help we can construct a rational transformation $\mathbb{T}(x,\epsilon)$ which converts a system of ordinary differential equations with rational coefficients given by the matrix $\mathbb{A}(x,\epsilon)$ to an equivalent system given by the matrix $\mathbb{M}(x,\epsilon)$ which is Fuchsian and has an epsilon form (also called canonical in [Hen13]), i.e., $\mathbb{M}(x,\epsilon)=\epsilon\,\mathbb{S}(x)$. When the epsilon form $\mathbb{M}(x,\epsilon)$ of the initial system $\mathbb{A}(x,\epsilon)$ is found we can easily solve it as a Laurent series in $\epsilon$ and restore solutions for the initial system $\mathbb{A}(x,\epsilon)$ — which is our ultimate goal — by solving a linear system of equations.

The whole method is performed in the following three steps:

Given a matrix $\mathbb{A}(x,\epsilon)$, find an equivalent system $\mathbb{A}^{\prime}(x,\epsilon)$ and a corresponding transformation $\mathbb{T}(x,\epsilon)$, such that $\mathbb{A}^{\prime}(x,\epsilon)$ is Fuchsian. We call this step fuchsification. 2. 2.

Given a Fuchsian matrix $\mathbb{A}(x,\epsilon)$ with eigenvalues of all its residues of the form $n+m\,\epsilon$, where $n$ is integer, find an equivalent system (along with the transformation) which is still Fuchsian, but with residue eigenvalues being all of the form $k\,\epsilon$. We call this step normalization. 3. 3.

Given a normalized matrix $\mathbb{A}(x,\epsilon)$, find an equivalent matrix $\mathbb{A}^{\prime}(x)$ in epsilon form, i.e. such that $\mathbb{A}^{\prime}(x,\epsilon)=\epsilon\,\mathbb{S}^{\prime}(x)$. We call this step factorization.

2.2.1 Fuchsification

To fuchsify a system (1) means to find an equivalent Fuchsian system. This, of course, is only possible if $\mathbb{A}(x,\epsilon)$ has no irregular singularities, or in other words $\mu_{k}(\mathbb{A})\leq 1$ for all $k$.333 We expect this to be often the case in practice, in particular for ODEs corresponding to Feynman integrals which are known to have only logarithmic singularities, hence be solutions of some Fuchsian ODEs.

In the case of a single ODE of order $n$, the minimal Moser order can be computed explicitly from power counting analysis of its coefficients (see the generalization of Fuchs’ theorem in [Mos59]). This is not possible for ODE systems like (1). Instead, we have a criterion for Moser-reducibility of the form:

Theorem 1. If $m_{k}(\mathbb{A})>1$ then the system (1) is Moser-reducible at $x=x_{k}$ if and only if the polynomial

[TABLE]

vanishes identically in $\lambda$ at $x=x_{k}$, where $r_{k}=\mathrm{rank}(\mathbb{A}_{k0})$ and $\Delta=x-x_{k}$ if $x_{k}\neq\infty$, or $\Delta=1/x$ if $x_{k}=\infty$.

When this condition fails the singularity in $x=x_{k}$ is irregular.

In addition to this criterion we have a method for constructing a transformation that lowers the Moser order at $x_{i}$, provided that the system is Moser-reducible at that point (possibly at the expense of increasing the Poincaré rank $p_{j}$ at another point $x=x_{j}$ by one). This is done by selecting a projector matrix $\mathbb{P}$ equal to a sum of products of a particular subset of (generalized) eigenvectors of $\mathbb{A}_{i0}$ and $\mathbb{A}_{j0}$, and constructing a $\mathbb{P}$-balance between $x_{i}$ and $x_{j}$.

The reader can find the details of this construction in [Lee15], but it is important to note that even if the system is Moser-reducible at $x_{i}$, it is only sometimes possible to construct a transformation that lowers $m_{i}(\mathbb{A})$ without increasing Poincaré rank at $x_{j}$. Sometimes the best we can do is to choose $x_{j}$ to be some regular point (with $p_{j}=-1$), and use a transformation that decreases $m_{i}(\mathbb{A})$ at the expense at increasing $p_{j}$ to [math], effectively introducing an apparent singularity where there was none before.444 In practice these apparent singularities are not a major problem, since they are subsequently removed during the normalization step. Still, we try not to introduce them if possible in order to decrease the intermediate expression sizes and to increase the overall performance.

With this in mind, to reduce a system to Fuchsian form we need to combine the reducibility check with the Moser rank-lowering transformation in stepwise fashion, as follows:

Select some point $x_{k}$ with $m_{k}(\mathbb{A})>1$. If none exist, reduction is complete. 2. 2.

Check if $\mathbb{A}(x,\epsilon)$ is reducible at $x=x_{k}$. If not, fail. 3. 3.

Find a transformation $\mathbb{T}(x,\epsilon)=\mathcal{B}(\mathbb{P}(\epsilon),x_{k},x_{j};x)$ that lowers $m_{k}(\mathbb{A})$. 4. 4.

Apply $\mathbb{T}(x,\epsilon)$ and repeat from Step 1.

In Fuchsia, this process is implemented by function fuchsify; see Section 3.3 for it’s usage.

Finally, let us mention that a similar problem of reducing Poincaré ranks and Moser orders of rational matrices was actively studied by Barkatou and co-authors, e.g., see [BP99]. They developed algorithms implementing their method [BP99] which is available in the standard Maple package DEtools as moser_reduce and super_reduce routines.

2.2.2 Normalization

To normalize a Fuchsian system (1) with matrix residue eigenvalues of the form $n+m\,\epsilon$ (where $n$ is integer), implies to find an equivalent Fuchsian system with residue eigenvalues of the form $m\,\epsilon$.

Just like in the previous step, the normalizing transformations are found by stepwise application of balance transformations (6). We refer the reader to [Lee15, p. 11] for the description of how such balances are constructed, but we will note that this transformation is possible due to these two facts:

• •

Given a properly selected projector matrix $\mathbb{P}(\epsilon)$, the balance $\mathcal{B}(\mathbb{P}(\epsilon),x_{i},x_{j};x)$ shifts one of the eigenvalues of $\mathbb{A}_{i0}$ by $\pm 1$, shifts one of the eigenvalues of $\mathbb{A}_{j0}$ by $\mp 1$, and does not change the Poincaré ranks at any point.

• •

Since the system is Fuchsian, the sum of all its matrix residues is zero, and thus, the sum of all residue eigenvalues is zero as well.

Combining these two facts, we perform the normalization by shifting the residue eigenvalues by $1$ at each step, until a state is reached where the integer parts of all the eigenvalues are zero.

In Fuchsia, the normalization step is implemented by the function normalize. See Section 3.3 for its usage.

It may happen that after fuchsification we obtain a system for which the residue eigenvalues do not fit to the $n+m\,\epsilon$ form neatly. In this case it is sometimes possible to rectify the problem by using some non-linear change of variables. Unfortunately we do not have an automated solution for such cases, and users are expected to find transformations appropriate for their system manually.

2.2.3 Factorization

After the normalization we have obtained a matrix $\mathbb{A}(x,\epsilon)$ with all residue eigenvalues of the form $m\,\epsilon$. The final step is to factorize it by finding an equivalent matrix which is by itself proportional to $\epsilon$, so $\mathbb{A}^{\prime}(x,\epsilon)=\epsilon\,\mathbb{S}(x)$.

A transformation which is constant in $x$ is sufficient for this task. Let $\mathbb{T}(\epsilon)$ be such a transformation, then according to (5), we have:

[TABLE]

or explicitly

[TABLE]

Since $\mathbb{S}_{i}$ in (11) is the same no matter what values $\epsilon$ takes, we can say that:

[TABLE]

from where we obtain a system of linear equations for $\mathbb{T}(\epsilon,\mu)\equiv\mathbb{T}(\epsilon)\,\mathbb{T}^{-1}(\mu)$:

[TABLE]

Solving this system for $\mathbb{T}(\epsilon,\mu)$ we can reconstruct the initial transformation as $\mathbb{T}(\epsilon)=\mathbb{T}(\epsilon,\mu_{0})$, where $\mu_{0}$ can be chosen arbitrary as long as $\mathbb{T}(\epsilon)$ will come out invertible.

In general, the solution for $\mathbb{T}(\epsilon,\mu)$ can have multiple free variables aside from just $\mu$. We choose to set all of them to random small integers, preferably zeros, which keeps the resulting matrix $\mathbb{S}(\epsilon)$ simple.

In Fuchsia the factorization step is implemented as the factorize routine (again, see Section 3.3 for its usage).

2.2.4 Block-triangular form

It often happens that a matrix which defines an ODEs is sparse, i.e. it has many zeros, and can be shuffled into block-triangular form with small blocks. It is possible to exploit this fact to considerably speed up fuchsification and normalization.

Let us consider a block-triangular matrix of this form:

[TABLE]

where $\mathbb{A}^{(ab)}(x,\epsilon)$ is a sub-matrix of size $n_{a}\times n_{b}$.

We start by reducing diagonal blocks of this matrix, $\mathbb{A}^{(ab)}(x,\epsilon)$, to epsilon form. This can be done by treating each diagonal block as an independent matrix, and proceeding as described in the previous sections. Since the characteristic polynomial of a block-triangular matrix is a product of characteristic polynomials of its diagonal blocks, i.e.,

[TABLE]

once we have normalized the diagonal blocks, the whole matrix becomes normalized as well (but not necessarily Fuchsian yet).

Next, we fuchsify off-diagonal blocks given by rectangular matrices $\mathbb{A}^{(ab)}(x,\epsilon)$, $a>b$. To do this, let us look at the parts of (1) related to such a block:

[TABLE]

If $\mathbb{A}^{(b\,a)}(x,\epsilon)$ has a singularity of Poincaré rank $r>0$ at $x=x_{k}$, then to reduce $r$ we can apply this basis transformation:

[TABLE]

where $\mathbb{D}$ is some constant matrix, and $\Delta=x-x_{k}$ if $x_{k}\neq\infty$, or $\Delta=1/x$ if $x_{k}=\infty$. This transformation changes $\mathbb{A}^{(b\,a)}_{k0}(\epsilon)$ into

[TABLE]

If both $\mathbb{A}^{(aa)}_{k0}(\epsilon)$ and $\mathbb{A}^{(bb)}_{k0}(\epsilon)$ have been factorized then it is always possible to solve the right-hand side of this equation for $\mathbb{D}$, thus reducing the Poincaré rank of $\mathbb{A}^{\prime(ba)}$ by one. Moreover, this transformation only affects $\mathbb{A}^{(bi)}$ for $i\leq a$ and $\mathbb{A}^{(ia)}$ for $i>b$, therefore if we will sequentially apply the transformation (17) for each off-diagonal block, starting from the first row to the last, and from the last column to the first, we will obtain a fully Fuchsian (and still normalized) matrix.

Finally, we factorize the whole matrix as described in Section 2.2.3, thus completing the transformation to epsilon form.

In Fuchsia the process of shuffling a matrix to its shortest lower block-diagonal form is performed by the block_triangular_form routine; fuchsification and normalization of diagonal blocks is done by the reduce_diagonal_blocks; fuchsification of the remaining off-diagonal blocks is done by the fuchsify_off_diagonal_blocks.

3 Using Fuchsia

3.1 Installation

To run Fuchsia one needs SageMath [Dev16] version 7.0 or higher to be installed on the computer. This task can be accomplished by following installation instructions available at the website http://www.sagemath.org. 555 Some Linux distributions have SageMath available in their package repositories; we do not recommend using those. A number of Maxima releases contain bugs which Fuchsia is sensitive to, and so far the official SageMath builds have avoided those releases (unlike some Linux distributions).

SageMath is a free and open-source Computer Algebra System licensed under GPL. It is written in Python 2.7 and combines together a number of existing open-source mathematical systems and libraries like Maxima, Singular, and others with the goal of providing the best free CAS. In particular our code heavily relies on the interface to Maxima [Max14].

We also should obtain a file fuchsia.py which can be downloaded from the website https://github.com/gituliar/fuchsia, where many examples of usage and reduced matrices are also collected.

3.2 Usage from the command line

To run Fuchsia use the command666For brevity, we use a shortcut fuchsia which is equivalent to sage -python fuchsia.py.

$ sage -python fuchsia.py <action> <options>where

• •

<action> is one of the algorithms described in the previous section, i.e., fuchsify, normalize, factorize, or auxiliary action transform, which applies a user-defined transformation to the given matrix.

• •

<options> are action-dependent options described in the help message printed with the help of fuchsia -h command.

In the following we provide a complete help information printed by fuchsia -h:

Fuchsia v16.11.14Authors: Oleksandr Gituliar, Vitaly MageryaUsage:    fuchsia [-hv] [--use-maple] [-f <fmt>] [-l <path>] [-P <path>]            <command> <args>...Commands:    reduce [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>        find an epsilon form of the given matrix    fuchsify [-x <name>] [-m <path>] [-t <path>] <matrix>        find a transformation that will transform a given matrix        into Fuchsian form    normalize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>        find a transformation that will transform a given Fuchsian        matrix into normalized form    factorize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>        find a transformation that will make a given normalized        matrix proportional to the infinitesimal parameter    sort [-m <path>] [-t <path>] <matrix>        find a block-triangular form of the given matrix    transform [-x <name>] [-m <path>] <matrix> <transform>        transform a given matrix using a given transformation    changevar [-x <name>] [-m <path>] <matrix> <expr>        transform matrix by susbtituting free variable by a        given expressionOptions:    -h          show this help message    -f <fmt>    matrix file format: mtx or m (default: mtx)    -l <path>   write log to this file    -v          produce a more verbose log    -P <path>   save profile report into this file    -x <name>   use this name for the free variable (default: x)    -e <name>   use this name for the infinitesimal parameter (default: eps)    -m <path>   save the resulting matrix into this file    -t <path>   save the resulting transformation into this file    --use-maple speed up calculations by using Maple when possibleArguments:    <matrix>    read the input matrix from this file    <transform> read the transformation matrix from this file    <expr>      arbitrary expression

Simple and more advanced results obtained with the help of Fuchsia are located in the directory examples. There you will find many examples of original matrices together with generated transformations which lead to the epsilon form of corresponding matrices. Another example of applying Fuchsia to find master integrals for next-to-leading order contributions to splitting functions in QCD were discussed in [GM16].

3.3 Usage from SageMath or Python

You can also use Fuchsia as a library by starting the SageMath prompt and importing the fuchsia module like this:

$ sage    ┌────────────────────────────────────────────────────────────────────┐    │ SageMath Version 7.1, Release Date: 2016-03-20                     │    │ Type "notebook()" for the browser-based notebook interface.        │    │ Type "help()" for help.                                            │    └────────────────────────────────────────────────────────────────────┘    sage: from fuchsia import *

In order to give an example for the API, let us try to reduce a simple matrix. For the list of functions available after import, please, read the next section.

sage: x, eps = var("x eps")    sage: M = matrix([    ....:   [(2-eps)/x, 0, 0],    ....:   [x/(x-1), eps/x, 0],    ....:   [(1+2*eps)/x**3, 0, (1+eps)/x/(x+1)]    ....: ])

First, let us see where the singularities of this matrix are located:

sage: singularities(M, x)    {-1: 0, 0: 2, 1: 0, +Infinity: 1}

So, 4 singularities in total, with the Poincaré rank being 2 at $\texttt{x}=0$, 1 at $\texttt{x}=\infty$ and 0 everywhere else. To get rid of non-zero ranks (thus transforming the system into Fuchsian form) we will need to fuchsify this matrix as:

sage: Mf, Tf = fuchsify(M, x)    sage: Mf    [ -(eps - 2)/x           0                          0]    [   -1/(x - 1)  (eps -1)/x                          0]    [(2eps + 1)/x           0  (eps + 2x + 3)/(x^2 + x)]    sage: singularities(Mf, x)    {-1: 0, 0: 0, 1: 0, +Infinity: 0}

Now, let us take a look at the eigenvalues of Mf residues:

sage: [matrix_residue(Mf, x, x0).eigenvalues()    ....:   for x0 in [-1, 0, 1, Infinity]]    [[-eps - 1, 0, 0],     [-eps + 2, eps - 1, eps + 3],     [0, 0, 0],     [-eps + 1, eps - 2, -2]]

Many of these eigenvalues are not equal to zero in the limit $\texttt{eps}\to 0$, so Mf is not normalized. It is, however, the case that all of the eigenvalues are of the form $n+m*\texttt{eps}$, so there is a chance that we will be able to normalize Mf. Let us try:

sage: Mn, Tn = normalize(Mf, x, eps)    sage: Mn    [-eps/x                                                            ...    [(4eps^3 - 8eps^2 - (4eps^2 - 6eps + 3)x + 5eps)/((4eps^3 - ...    [((2eps + 1)x + 3eps + 1)/(x^2 + x)                             ...    sage: [matrix_residue(Mn, x, x0).eigenvalues()    ....:   for x0 in [-1, 0, 1, Infinity]]    [[-eps, 0, 0],     [-eps, eps, eps],     [0, 0, 0],     [-eps, eps, 0]]

So, the matrix is normalized, but it grew quite a bit larger. This happens. Sometimes it is possible to simplify it a bit:

sage: Ms, Ts = simplify_by_jordanification(Mn, x)    sage: Ms    [               -eps/x      0              0]    [            1/(x - 1)  eps/x              0]    [1/2*(eps + 1)/(x + 1)      0  eps/(x^2 + x)]

That is much better.

Finally, we need to factorize Ms to complete the reduction:

sage: Mr, Tr = factorize(Ms, x, eps)    sage: Mr    [         -eps/x      0              0]    [1/4eps/(x - 1)  eps/x              0]    [5/8eps/(x + 1)      0  eps/(x^2 + x)]

This is the fully transformed matrix. As you can see, it is both proportional to eps and Fuchsian. To make sure we got everything right, we can double-check the full transformation:

sage: T = (TfTnTs*Tr).simplify_rational()    sage: (Mr - transform(M, x, T)).is_zero()    True

Note that we have used the construct (A - B).is_zero() to compare matrices instead of the more obvious bool(A == B). This is a SageMath idiosyncrasy; the more obvious way compares symbolic matrices only structurally.

Of course, you do not need to walk through all these steps yourself every time. Normally, you just need to call this one function to do all the reduction work:

sage: MM, TT = epsilon_form(M, x, eps)    sage: MM    [                          -eps/x      0              0]    [                 1/4eps/(x - 1)  eps/x              0]    [1/4(9epsx + 13*eps)/(x^2 + x)      0  eps/(x^2 + x)]    sage: (MM - transform(M, x, TT)).is_zero()    True

Notice that this matrix is slightly more complex than the one we have obtained step by step above. This also happens. The final form we are computing is not unique, and it will be different depending on the precise sequence of reduction steps you have taken.

Additionally, many of the transformations take a special $seed$ parameter to control the order of operations they perform internally. By supplying different seeds, you will obtain different results as well.

3.3.1 Function reference

$\texttt{epsilon\_form}(\mathbb{M},x,epsilon,seed=0)$** (function)**

Fully reduces a system of equations defined by a matrix $\mathbb{M}$, an independent variable $x$, and an infinitesimal parameter $epsilon$. Returns a pair of values: the transformed matrix $\mathbb{M}^{\prime}$ and the transformation matrix $\mathbb{T}$. Raises $FuchsiaError$, if the system is irreducible.

The reduction is performed by first converting $\mathbb{M}$ to block-triangular form, then reducing the diagonal blocks via fuchsify, normalize and factorize, reducing off-diagonal blocks as described in Section 2.2.4, and finally factorizing $epsilon$ via factorize.

$\texttt{fuchsify}(\mathbb{M},x,seed=0)$** (function)**

Reduces a system defined by a matrix $\mathbb{M}$ and an independent variable $x$ to Fuchsian form. That is, it makes sure that the Poincaré ranks of all singularities of the transformed matrix $\mathbb{M}^{\prime}$ are [math]. Returns a pair of values: the transformed matrix $\mathbb{M}^{\prime}$ and the transformation $\mathbb{T}$. If the system is irreducible, it raises $FuchsiaError$.

$\texttt{normalize}(\mathbb{M},x,epsilon,seed=0)$** (function)**

Transforms a Fuchsian system defined by a matrix $\mathbb{M}$, an independent variable $x$ and, an infinitesimal parameter $epsilon$ to a normalized form. That is, it makes sure that real parts of the eigenvalues of all matrix residues of the transformed matrix $\mathbb{M}^{\prime}$ lie in the range $[-1/2,1/2)$ in the limit $epsilon\to 0$. Returns a pair of values: the transformed matrix $\mathbb{M}^{\prime}$ and the transformation $\mathbb{T}$. If such a transformation can not be found, it raises $FuchsiaError$.

$\texttt{factorize}(\mathbb{M},x,epsilon,seed=0)$** (function)**

Transforms a normalized system defined by a matrix $\mathbb{M}$, an independent variable $x$, and an infinitesimal parameter $epsilon$ so that the transformed matrix $\mathbb{M}^{\prime}$ is proportional to $epsilon$. Returns a pair of values: the transformed matrix $\mathbb{M}^{\prime}$ and the transformation $\mathbb{T}$. If such a transformation can not be found, it raises $FuchsiaError$.

$\texttt{block\_triangular\_form}(\mathbb{M})$** (function)**

Transforms a matrix $\mathbb{M}$ into a lower block-triangular form.

Returns three values: a transformed matrix $\mathbb{M}^{\prime}$, a transformation matrix $\mathbb{T}$, and a list of tuples $(m_{i},n_{i})$, where $n_{i}$ is the size of $i$-th diagonal block, and $m_{i}=\sum_{j=1}^{i-1}n_{j}$. The tuple list represents block structure of $\mathbb{M}^{\prime}$; it is used by the next two functions, where it is passed as the argument $B$.

$\texttt{reduce\_diagonal\_blocks}(\mathbb{M},x,epsilon,B=None,seed=0)$** (function)**

Finds a transformation that reduces diagonal blocks of $\mathbb{M}$ into epsilon form. If $B$ is not provided, it transforms $\mathbb{M}$ into lower triangual form before reduction. Otherwise it assumes $\mathbb{M}$ blocks are described by $B$.

Returns two values: a transformed matrix $\mathbb{M}$ and a transformation matrix $\mathbb{T}$.

$\texttt{fuchsify\_off\_diagonal\_blocks}(\mathbb{M},x,epsilon,r=None)$** (function)**

Given a matrix $M$ with diagonal blocks in epsilon form, it transforms off-diagonal blocks in to Fuchsian form. If $B$ is not provided, it transforms $\mathbb{M}$ into lower triangual form before reduction. Otherwise it assumes that blocks of $\mathbb{M}$ are described by $B$.

Returns two values: a transformed matrix $\mathbb{M}$ and a transformation matrix $\mathbb{T}$.

$\texttt{simplify\_by\_jordanification}(\mathbb{M},x)$** (function)**

Tries to simplify a system defined by a matrix $\mathbb{M}$ and an independent variable $x$ by constant transformations that transform leading expansion coefficients of $\mathbb{M}$ into their Jordan forms. Returns a pair of values: the simplified matrix $\mathbb{M}^{\prime}$ and the transformation $\mathbb{T}$. If none of the attempted transformations reduces the complexity of $\mathbb{M}$ (as measured by matrix_complexity), it returns the original matrix and the identity transformation.

$\texttt{simplify\_by\_factorization}(\mathbb{M},x)$** (function)**

Tries to simplify a system defined by a matrix $\mathbb{M}$ and an independent variable $x$ by a constant transformation that extracts common factors found in $\mathbb{M}$ (if any). Returns a pair of values: the simplified matrix $\mathbb{M}^{\prime}$ and the transformation $\mathbb{T}$.

$\texttt{matrix\_complexity}(\mathbb{M})$** (function)**

This function is used as a measure of matrix complexity by fuchsify and simplification functions. Currently it is defined as the length of textual representation of matrix $\mathbb{M}$.

$\texttt{balance}(\mathbb{P},x_{1},x_{2},x)$** (function)**

Returns a balance transformation between points $x=x_{1}$ and $x=x_{2}$ using the projector matrix $\mathbb{P}$. See eq. (6) for the definition of a balance.

$\texttt{transform}(\mathbb{M},x,\mathbb{T})$** (function)**

Transforms a system defined by a matrix $\mathbb{M}$ and an independent variable $x$ using the transformation matrix $\mathbb{T}$ as specified by eq. (5). Returns the transformed matrix $\mathbb{M}^{\prime}$.

$\texttt{balance\_transform}(\mathbb{M},\mathbb{P},x_{1},x_{2},x)$** (function)**

Same as $\texttt{transform}(\mathbb{M},x,\texttt{balance}(\mathbb{P},x_{1},x_{2},x))$, but implemented more efficiently: since the inverse of $\texttt{balance}(\mathbb{P},x_{1},x_{2},x)$ is $\texttt{balance}(\mathbb{P},x_{2},x_{1},x)$, this function can avoid a time-consuming matrix inversion operation that transform must perform.

$\texttt{singularities}(\mathbb{M},x)$** (function)**

Finds values of $x$ around which the matrix $\mathbb{M}$ has a singularity in $x$. Returns a dictionary with $\{x_{i}:p_{i}\}$ entries, where $p_{i}$ is the Poincaré rank of $\mathbb{M}$ at $x=x_{i}$. The set of singular points can include Infinity, if $\mathbb{M}$ has a singularity at $x\to\infty$.

$\texttt{matrix\_c0}(\mathbb{M},x,x_{0},p)$** (function)**

Returns the 0-th coefficient of the series expansion of a matrix $\mathbb{M}$ around $x=x_{0}$, assuming the Poincaré rank of $\mathbb{M}$ at that point is $p$. If $x_{0}$ is Infinity, it returns minus the coefficient at the highest power of $x$. In other words, it returns $\mathbb{A}_{k0}$ from eq. (2).

$\texttt{matrix\_c1}(\mathbb{M},x,x_{0},p)$** (function)**

Returns the 1-th coefficient of the series expansion of a matrix $\mathbb{M}$ around $x=x_{0}$, assuming the Poincaré rank of $\mathbb{M}$ at that point is $p$. If $x_{0}$ is Infinity, it returns minus the coefficient at the second-to-highest power of $x$. In other words, it returns $\mathbb{A}_{k1}$ from eq. (2).

$\texttt{matrix\_residue}(\mathbb{M},x,x_{0})$** (function)**

Returns a residue of a matrix $\mathbb{M}$ at $x=x_{0}$, assuming that the Poincaré rank of $\mathbb{M}$ at $x=x_{0}$ is [math]. Returns matrix residue at infinity if $x_{0}$ is Infinity.

This is the same as $\texttt{matrix\_c0}(\mathbb{M},x,x_{0},0)$.

$\texttt{export\_matrix\_to\_file}(filename,\mathbb{M},fmt=\texttt{"mtx"})$** (function)**

Writes a matrix $\mathbb{M}$ to a file $filename$ using MatrixMarket array format if $fmt$ is "mtx" (which is the default), or Mathematica format if $fmt$ is "m".

$\texttt{import\_matrix\_from\_file}(filename)$** (function)**

Reads a symbolic matrix from a file $filename$. Both Mathematica and MatrixMarket array formats are supported. The exact format will be autodetected.

$\texttt{setup\_fuchsia}(verbosity=0,use\_maple=\texttt{False})$** (function)**

Modifies some of the Fuchsia inner workings. In particular, it sets $verbosity$ to $2$ to enable verbose logging, $1$ to enable normal logging, and [math] to only log errors.

Set $use\_maple$ to True to enable usage of Maple to speed up calculations when possible; this may be particularly beneficial for big matrices, and for matrices with singularities at complex points.

$FuchsiaError$** (class)**

This is the class of exceptions raised by the Fuchsia routines. It indicates the inability to perform the requested reduction.

4 Summary

In this paper we have presented Fuchsia, a program for reducing differential equations for Feynman master integrals to the epsilon form based on the Lee algorithm [Lee15] which consists of three main computational steps: fuchsification, normalization, and factorization. Fuchsia is open-source nature and depends on free software tools only: the programming language Python and the computer algebra system Maxima/Sage, which makes it available for everyone. Unfortunately, for some heavy-duty computations one needs to switch from Maxima to Maple, which is more powerful for some tasks, hence the access to latter is desirable. Though it helps a lot in some situation, Maple is not a universal solution and is not able to process huge intermediate expressions which arise during the reduction process of some advanced examples. It happens due to inability to work with polynomials which contain algebraic numbers in their coefficients (even in the case of complex coefficients the performance drastically decreases).

Despite of the discussed limitation Fuchsia shows a great performance in many cases. It is possible due to the optimization for block-triangular (or sparse) matrices, which allows to reduce relatively large matrices: the reduction of ${74\times 74}$ matrix with 20 rational and complex singular points and at most $3\times 3$ coupled blocks takes about an hour on a laptop with Intel i5 CPU.

Acknowledgment

We are gratefully thankful for advanced examples of differential equations provided by Roman Lee and Costas Papadopolous. We also appreciate useful suggestions from Sven Moch during writing of this paper.

This work has been supported by the Deutsche Forschungsgemeinschaft in Sonderforschungsbereich 676 Particles, Strings, and the Early Universe and by the Narodowe Centrum Nauki with the Sonata Bis grant DEC-2013/10/E/ST2/00656.