Algorithmic transformation of multi-loop master integrals to a canonical basis with CANONICA

TL;DR

This paper introduces CANONICA, a Mathematica package that automates transforming multi-loop Feynman integrals into a canonical basis, simplifying their differential equations, and is the first publicly available tool for multi-scale cases.

Contribution

The paper presents the first publicly available implementation of an algorithm for transforming multi-scale Feynman integrals into a canonical basis.

Findings

CANONICA automates the transformation process.

The algorithm's uniqueness of canonical forms is proven.

The package facilitates solving differential equations of Feynman integrals.

Abstract

The integration of differential equations of Feynman integrals can be greatly facilitated by using a canonical basis. This paper presents the Mathematica package CANONICA, which implements a recently developed algorithm to automatize the transformation to a canonical basis. This represents the first publicly available implementation suitable for differential equations depending on multiple scales. In addition to the presentation of the package, this paper extends the description of some aspects of the algorithm, including a proof of the uniqueness of canonical forms up to constant transformations.

Equations

$$\text{d}\vec{f}(ϵ,\{x_j\})=a(ϵ,\{x_j\})\vec{f}(ϵ,\{x_j\}),$$

$$a(ϵ,\{x_j\})=\sum _{i=1}^M{a}_i(ϵ,\{x_j\})\text{d}{x}_i.$$

$$\text{d}a-a\wedge a=0,$$

$$\vec{f}=T(ϵ,\{x_j\}){\vec{f}}^{\prime },$$

$$a^{\prime }=T^{-1}aT-T^{-1}\text{d}T.$$

$$a^{\prime }(ϵ,\{x_j\})=ϵ\text{d}\tilde{A}(ϵ,\{x_j\})=ϵ\sum _{l=1}^N{\tilde{A}}_l\text{d}\log (L_l(\{x_j\})),$$

$$ϵ\text{d}\tilde{A}=T^{-1}aT-T^{-1}\text{d}T,$$

$$\text{Tr}[a]=ϵX(\{x_j\})+Y(ϵ,\{x_j\}),$$

$$\det (T)=F(ϵ)\exp \left({\int }_{\gamma }Y(ϵ,\{x_j\})\right),$$

$$\text{Tr}[\text{d}\tilde{A}]=X(\{x_j\}).$$

$$\text{d}T-aT+ϵT\text{d}\tilde{A}=0.$$

$$\hat{a}=\sum _{k=0}^{k_{\text{max}}}{ϵ}^k{\hat{a}}^{(k)}.$$

$$h(ϵ,\{x_j\})=\sum _{l=l_{\text{min}}}^{l_{\text{max}}}{ϵ}^l{h}^{(l)}(\{x_j\}),l_{\text{min}}\ge 0.$$

$$-\hat{T}\text{d}h+h\text{d}\hat{T}-\hat{a}\hat{T}+ϵh\hat{T}\text{d}\tilde{A}=0.$$

$$\hat{T}=\sum _{n=l_{\text{min}}}^{n_{\text{max}}}{ϵ}^n{\hat{T}}^{(n)},$$

$${\hat{T}}^{(n)}=\sum _{k=1}^{|{\mathcal{R}}_T|}{\tau }_k^{(n)}{r}_k(\{x_j\}),$$

$${\mathcal{R}}_T=\left\{r_1(\{x_j\}),\dots ,r_{|{\mathcal{R}}_T|}(\{x_j\})\right\},$$

$$\tilde{A}=\sum _{l=1}^N{\alpha }_l\log (L_l(\{x_j\})),$$

$$a_{I}=\left(\,\begin{array}[]{|ccc|c|}\cline{1-4}\cr&&&\\ &\epsilon\tilde{c}&&\,\,0\\ &&&\\ \hline\cr&b&&\epsilon\tilde{e}\\ \hline\cr\end{array}\,\right),$$

$$t_{D}=\left(\,\begin{array}[]{|ccc|c|}\cline{1-4}\cr&&&\\ &\mathbb{I}&&\,\,0\\ &&&\\ \hline\cr&D&&\mathbb{I}\\ \hline\cr\end{array}\,\right),$$

$$\text{d}D-ϵ(\tilde{e}D-D\tilde{c})=b-b^{\prime }.$$

$$S\sim \frac{1}{f^n}$$

$$S=R\frac{1}{f^n}.$$

$$\hat{T}\sim \frac{1}{f^n},n\ge 1$$

$$\hat{a}\sim \frac{1}{f^k},k\ge 1.$$

$$-\hat{T}\text{d}h+h(\text{d}\hat{T}+ϵ\hat{T}\text{d}\tilde{A})=\hat{a}\hat{T}.$$

$$\hat{T}\text{d}h\sim \frac{1}{f^k},k\le n.$$

$${\hat{T}}^{(s)}\sim \frac{1}{f^n},$$

$${\hat{T}}^{(s-1)}\sim \frac{1}{f^k},k<n.$$

$$\text{d}{\hat{T}}^{(s)}\sim \frac{1}{f^{n+1}}.$$

Full text

††journal: Computer Physics Communications

HU-EP-17/10

**PROGRAM SUMMARY

** Program Title: *CANONICA

Licensing provisions:* GNU General Public License version 3

Programming language: Wolfram Mathematica, version 10 or higher

Nature of problem: Computation of a rational basis transformation of master integrals leading to a canonical form of the corresponding differential equation.

Solution method: The transformation law is expanded in the dimensional regulator. The resulting differential equations for the expansion coefficients of the transformation are solved with a rational ansatz.

1 Introduction

The calculation of higher order corrections to the cross-sections measured at the LHC is crucial in order to improve the understanding of both the background reactions as well as the signal processes. The current state of the art are NNLO QCD corrections to $2\rightarrow 2$ processes involving a limited number of mass scales. A major challenge in these computations is the evaluation of the occurring Feynman integrals. While the calculation of Feynman integrals can be attempted with numerous approaches, the method of differential equations [1, 2, 3] has been particularly successful in the recent years [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. This success is due to the observation [4] that the differential equation can be simplified significantly by turning to a so-called canonical basis of master integrals. The differential equation of a canonical basis of master integrals can easily be integrated in terms of iterated integrals such as multiple polylogarithms [39, 40].

It is well known that Feynman integrals exist [41, 42, 43, 44, 45, 46, 47, 48, 49, 32], which do not evaluate to this class of functions. These integrals generally satisfy differential equations of higher order. The solutions of the homogeneous part of these equations have been shown to be constructible by evaluating unitarity cuts [50, 51, 52, 53, 54]111In fact, unitarity cuts have also been used to derive differential equations of Feynman integrals [55].. Some integrals, which exceed the class of multiple polylogarithms, have recently been shown to be iterated integrals of modular forms [56]. However, the concept of a canonical basis has not yet been extended to integrals of this kind.

The class of Feynman integrals, which do admit a canonical basis, is still large and contains many integrals of phenomenological interest. It is therefore desirable to automate the calculation of these integrals as much as possible. The systematic application [57, 58] of integration by parts relations [59, 60] to reduce all scalar integrals to a finite number of master integrals has been automated in a variety of publicly available tools [61, 62, 63, 64, 65, 66, 67, 68, 69]. This leaves the process of choosing a canonical basis as the next step to be automated. A number of different methods [4, 8, 10, 11, 14, 70, 71, 30, 72, 73] have been proposed to construct such a basis or, equivalently, the transformation from a given basis to a canonical basis. Until now, only implementations of the algorithm presented in Ref. [70] are publicly available [74, 75, 76]. However, this algorithm is restricted to ordinary differential equations, which are not sufficient to describe the full functional dependence of Feynman integrals depending on multiple dimensionless scales. A wide class of phenomenologically relevant integrals is thus not covered.

This paper aims to overcome this restriction by introducing an implementation of the algorithm in Ref. [72], which is applicable to multi-scale problems. The accompanying Mathematica package CANONICA allows to calculate a rational transformation to a canonical basis for a given differential equation. In addition, the package provides some supplemental functionality for handling differential equations of Feynman integrals.

The description of the algorithm in Ref. [72] is extended in the present paper by a detailed account of the construction of the set of rational functions used for the ansatz. Moreover, the occurrence of non-linear polynomial equations in the parameters of the ansatz is addressed with a procedure to extract all relevant information by solving only linear equations, while maintaining all of the algorithms generality. The latter relies on the uniqueness of canonical forms up to constant transformations. While it is trivial to show that a constant transformation of a given canonical form leads again to a canonical form, it is not obvious that all possible canonical forms can be obtained in this way, which will be proven in this paper.

The paper is organized as follows. In Section 2 the algorithm in Ref. [72] is briefly reviewed and a description of the procedure for the generation of the ansatz is given. Furthermore, the uniqueness of canonical forms up to constant transformations is proven in this section. Building on this result, the treatment of non-linear polynomial equations in the parameters of the ansatz is discussed. Section 3 introduces the CANONICA package by outlining the installation and the contents of the package, which is followed by a few examples illustrating the usage of the main features of CANONICA. Furthermore, an overview over the hierarchy of the main public functions is given. The conclusions are drawn in Section 4. A brief description of all functions and options provided by the package is contained in A and B. The global variables and protected symbols of the package are listed in C.

2 Algorithm

This section briefly reviews the algorithm introduced in Ref. [72] and presents more details on some aspects of the algorithm. In particular, the procedure used by CANONICA to generate an ansatz is described in detail. Furthermore, canonical forms are proven to be unique up to constant transformations. This result allows to attribute the occurrence of non-linear equations in the parameters of the ansatz precisely to this ambiguity. On this basis, a procedure to calculate the transformation by solving only linear equations is outlined.

2.1 Preliminaries

Let $\vec{f}(\epsilon,\{x_{j}\})$ denote the $m$-dimensional vector of master integrals, which depends on the dimensional regulator $\epsilon$ and a set $\{x_{j}\}$ of dimensionless invariants. By taking the total derivative of the vector of master integrals with respect to the invariants and expressing the result as a linear combination of master integrals, a coupled system of differential equations is obtained:

[TABLE]

with

[TABLE]

Here the $a_{i}(\epsilon,\{x_{j}\})$ denote $m\times m$ matrices of rational functions in the invariants and $\epsilon$. Using the linear independence of the master integrals over the field of rational functions in the invariants and taking the exterior derivative of Eq. (1) implies the following integrability condition:

[TABLE]

which is a valuable consistency check for differential equations in several variables. Transforming the basis of master integrals with an invertible transformation $T$,

[TABLE]

as suggested in Ref. [4], leads to the following transformation law for $a(\epsilon,\{x_{j}\})$:

[TABLE]

It has been observed [4] that with an appropriate change of the basis of master integrals, the differential equation can often be cast in the following form:

[TABLE]

where the $\tilde{A}_{l}$ denote constant $m\times m$ matrices and the functions $L_{l}(\{x_{j}\})$ are called letters. The above form of the differential equation is called canonical- or $\epsilon$-form. In this form, the integration of the differential equation in terms of iterated integrals is reduced to a merely combinatorial task (c.f. e.g., [6, 9, 10, 11, 12]).

2.2 Review of the algorithm

In this section the algorithm presented in Ref. [72] is briefly reviewed. Throughout this section, the existence of a rational transformation of the differential equation into $\epsilon$-form is assumed. The purpose of the algorithm is to compute this transformation for a given differential equation, provided it exists. Any transformation to a canonical basis has to satisfy the following equation:

[TABLE]

for some $\textup{d}\tilde{A}$ of the form in Eq. (6). The resulting differential form $\textup{d}\tilde{A}$ is unknown and thus has to be determined as well. In Ref. [72] it has been proven that the determinant of the transformation is fixed up to a rational function $F(\epsilon)$ by the trace of the differential form $a$

[TABLE]

where $X(\{x_{j}\})$ denotes the sum of dlog-terms with coefficients proportional to $\epsilon$ and $Y(\epsilon,\{x_{j}\})$ denotes the sum of dlog-terms with constant coefficient. Then the determinant is given by

[TABLE]

and the trace of the resulting $\epsilon$-form is determined by

[TABLE]

For invertible transformations $T$, Eq. (7) can equivalently be written as

[TABLE]

The basic idea of the algorithm is to expand this equation in $\epsilon$ and solve for the expansion coefficients of the transformation with a rational ansatz. However, the expansion of $T$ may not be finite and therefore an additional step has to be taken, which is reviewed in the following. As $a(\epsilon,\{x_{j}\})$ is required to be rational in both the invariants and $\epsilon$, a polynomial $h(\epsilon,\{x_{j}\})$ exists such that $\hat{a}=ah$ has a finite Taylor expansion in $\epsilon$

[TABLE]

In order to fix $h(\epsilon,\{x_{j}\})$ up to an irrelevant constant factor, $h(\epsilon,\{x_{j}\})$ is required to satisfy the above condition with the smallest possible number of irreducible factors. The expansion of $h$ is denoted as follows:

[TABLE]

Rewriting Eq. (11) in terms of $\hat{T}=Th$ yields the following equation:

[TABLE]

It has been shown in Ref. [72] that Eq. (14) will have a solution for $\hat{T}$ with finite expansion, if Eq. (11) has a rational solution for $T$. This allows to expand $\hat{T}$ in $\epsilon$

[TABLE]

and solve Eq. (14) order by order in $\epsilon$ for finitely many coefficients $\hat{T}^{(n)}$. The equations at each order are solved by making an ansatz for $\hat{T}^{(n)}$ in terms of rational functions of the invariants $r_{k}(\{x_{j}\})$

[TABLE]

[TABLE]

where the $\tau_{k}^{(n)}$ denote unknown $m\times m$ matrices independent of the invariants and the regulator, which are to be determined by the algorithm. More details on the choice of the set of rational functions $\mathcal{R}_{T}$ are given in Section 2.3. For the unknown $\tilde{A}$ an ansatz of the following form can be used:

[TABLE]

where the $\alpha_{l}$ are considered to be unknown $m\times m$ matrices independent of the invariants and the regulator. The set of polynomials $L_{l}(\{x_{j}\})$ is taken to be the set of irreducible denominator factors of the differential form $a(\epsilon,\{x_{j}\})$ with trivial dependence on the regulator. In Section 2.4, this set is shown to contain all letters of the resulting canonical form.

Inserting the ansatz in the expansion of Eq. (14) and requiring the resulting equations to hold for all non-singular values of the invariants implies polynomial equations in the parameters of the ansatz. For more details on the solution of these equations, see Section 2.7.

It is well known that the differential form $a(\epsilon,\{x_{j}\})$ can be cast in a block-triangular form. This allows to split the computation of the transformation into a recursion over the sectors of the differential equation, which leads to significant performance improvements. With regard to the recursion step, consider a differential equation where all previous sectors have already been transformed into $\epsilon$-form. The first part of the recursion step is to transform the diagonal block of the next sector into $\epsilon$-form with the part of the algorithm described above. After this step, the differential equation is in the following form:

[TABLE]

where $\tilde{c}$ and $\tilde{e}$ are in dlog-form. It has been shown in Ref. [72] that the transformation of a differential equation in this form can be split into two parts. First, the off-diagonal block $b$ is transformed into dlog-form with a transformation of the form

[TABLE]

which is determined by a differential equation for $D$

[TABLE]

Here $b^{\prime}$ is an unknown quantity, which is required to be in dlog-form. The above equation is solved by first multiplying $D$ by appropriate factors to render its expansion finite. The expansion coefficients are then determined by making an ansatz in terms of a set of rational functions $\mathcal{R}_{D}$. For more details, see [72] and Section 2.5, which describes how the set $\mathcal{R}_{D}$ is constructed.

The last part of the recursion step is to employ a procedure proposed in Ref. [70] to compute a rational transformation in the regulator, which transforms the full differential equation into $\epsilon$-form.

2.3 Ansatz for diagonal blocks

In this section, the choice of the set $\mathcal{R}_{T}$ of rational functions in the ansatz is discussed. The basic compromise with the ansatz is to choose it large enough to encompass the solution and as small as possible in order to keep the resulting number of equations small and therefore allow the algorithm to perform well. The goal of this section is to present a procedure to generate a finite set of rational functions for a given differential form $a(\epsilon,\{x_{j}\})$, which can then be used as an ansatz. The first step towards this goal is to determine the set of possible denominator factors of the transformation to a canonical basis. A natural guess is to consider the set of irreducible denominator factors of $\hat{a}$, which is proven in the following to contain all possible factors.

It is useful to first define some notation. Let $f(\{x_{j}\})$ be an irreducible polynomial and $S(\epsilon,\{x_{j}\})$ some matrix-valued rational differential form or function. Then, the notation

[TABLE]

indicates222Throughout this paper, the number 0 is understood to be included in the natural numbers. $n\in\mathbb{N}$ to be the maximal number for which $S$ can be written as

[TABLE]

Here $R$ is required to be nonzero and not to be the product of $f$ and a quantity, which is finite on the set of all zeros of $f$. The set $\mathcal{I}(S)$ of irreducible denominator factors of $S$ is then given by those factors $f$ with $S\sim 1/f^{k}$ and $k\geq 1$.

Claim 1**.**

Each irreducible denominator factor $f(\{x_{j}\})$ of a rational solution $\hat{T}$ of Eq. (14) is an irreducible denominator factor of $\hat{a}$.

The assertion in the claim is equivalent to $\mathcal{I}(\hat{T})\subseteq\mathcal{I}(\hat{a})$, which will be proven by showing that

[TABLE]

implies

[TABLE]

To this end, Eq. (24) is assumed to hold. It is instructive to rearrange the terms in Eq. (14)

[TABLE]

Since $\textup{d}h$ is polynomial, the term $\hat{T}\textup{d}h$ behaves as

[TABLE]

Given Eq. (24), there must be a lowest order $s\in\mathbb{Z}$ in the $\epsilon$-expansion of $\hat{T}$ with

[TABLE]

and consequently

[TABLE]

The derivative raises the power of $f$ by one, which implies

[TABLE]

Note that $\textup{d}\tilde{A}$ is in dlog-form, and therefore

[TABLE]

Taking both Eq. (29) and Eq. (31) into account, it follows

[TABLE]

This implies for the order $s$ of the expansion of $(\textup{d}\hat{T}+\epsilon\hat{T}\textup{d}\tilde{A})$

[TABLE]

which also holds for the full expression $(\textup{d}\hat{T}+\epsilon\hat{T}\textup{d}\tilde{A})$. Due to the minimality requirement, $h$ does not contain any irreducible factors independent of $\epsilon$ and therefore the multiplication of the term in brackets with $h$ in Eq. (26) cannot cancel any power of $f$, because $f$ is independent of $\epsilon$. Since Eq. (27) shows that the term $\hat{T}\textup{d}h$ is of lower order in $f$ than $h(\textup{d}\hat{T}+\epsilon\hat{T}\textup{d}\tilde{A})$, the whole left-hand side of Eq. (26) is of order $1/f^{n+1}$ and consequently the right-hand side as well:

[TABLE]

Since $\hat{T}$ is only of order $1/f^{n}$, it can be concluded

[TABLE]

which proves the claim. Thus, the ansatz can without loss of generality be restricted to the set

[TABLE]

of rational functions with the denominator factors drawn from the set $\mathcal{I}(\hat{a})=\{f_{1},\dots,f_{U}\}$ of irreducible denominator factors of $\hat{a}$.

As was argued in Ref. [72], rational functions may be decomposed in terms of a class of simpler rational functions, called Leinartas functions [77, 78]. Let $\mathcal{L}(\mathcal{Q})$ denote a basis of the $K$-span of $\mathcal{Q}$ in terms of Leinartas functions. While $\mathcal{L}(\mathcal{Q})$ is guaranteed to contain the correct ansatz, it is still an infinite set. Therefore, a constructive procedure is needed to generate a finite subset of $\mathcal{L}(\mathcal{Q})$ for a given $a(\epsilon,\{x_{j}\})$. This procedure should be inexpensive to compute while yielding a correct ansatz for most practical examples. Since the procedure outlined in the following is not proven to generate a correct ansatz, it is important to be able to systematically enlarge the ansatz in a way that is guaranteed to eventually encompass the solution.

The strategy to define a finite subset of $\mathcal{L}(\mathcal{Q})$ is to set restrictions on the powers of the invariants in the numerator as well as on the powers of the denominator factors.

While the powers of those factors occurring in $a(\epsilon,\{x_{j}\})$ may be suspected to be a good indicator for the powers in the transformation, the following simple example demonstrates this to be false. Consider the differential equation

[TABLE]

which contains the factor $x$ with the negative power one. However, for any given integer $\alpha$, the rational transformation to the canonical form reads

[TABLE]

Consequently, the transformation can contain any power of the factor $x$, while the power of the same factor in the differential equation remains fixed. A much better predictor is given by the determinant of the transformation, which in the 1-dimensional example above is identical to the transformation itself and therefore always yields the correct power of the factor $x$. For higher dimensional differential equations, the determinant does not fix the transformation but still carries information on the powers of the irreducible denominator factors of the transformation. Let the determinant of the transformation read

[TABLE]

where $F(\epsilon,\{x_{j}\})$ denotes the product of all irreducible factors with a non-trivial dependence on the regulator. Then, for each factor $f_{i}$ with $\lambda_{i}>0$, there has to be a component $T_{jl}$ of the transformation satisfying

[TABLE]

where $m$ denotes the dimension of the differential equation. Thus, the determinant sets lower bounds on the maximal powers of the denominator factors in the transformation, which have to be taken into account in the construction of the ansatz.

In the following, a finite subset of $\mathcal{Q}$ will be constructed, which then leads to a finite subset of $\mathcal{L}(\mathcal{Q})$ by taking a basis of its $K$-span in terms of Leinartas functions. The powers $\lambda_{i}$ are used to define a set of denominators:

[TABLE]

which has been restricted to at most $M$ denominator factors with $M$ denoting the number of invariants. Any higher number of polynomials in $M$ invariants is algebraically dependent and therefore reducible in terms of Leinartas functions with $M$ or less denominator polynomials. The parameter $\delta_{D}\in\mathbb{N}$ has been introduced to define a way to enlarge the set of $\mathcal{D}(\delta_{D})$ systematically. The default value is going to be $\delta_{D}=0$. The lower bounds in Eq. (40) are satisfied for all allowed values of $\delta_{D}$. For the numerators, consider the set of all possible monomials up to a fixed bound on their total degree

[TABLE]

where the parameter $\delta_{N}$ has been introduced to control the highest total degree of the monomials in $\mathcal{N}(\delta_{N})$. For the default value $\delta_{N}=0$ the highest total degree of the numerator monomials is three. This choice is made based on practical examples and is intended to make the default value $\delta_{N}=0$ work for most cases and at the same time yield a rather small ansatz. Furthermore, it has proven useful to also include the following sets of monomials:

[TABLE]

[TABLE]

in order to capture the correct ansatz in more cases already with the default value $\delta_{N}=0$. Usually, the inclusion of $\mathcal{N}_{\textup{det}}$ and $\mathcal{N}_{a}$ does not significantly enlarge the ansatz, while making the default value work for more examples. Finally, the ansatz $\mathcal{R}_{T}$ is obtained by computing a basis of Leinartas functions of the $K$-span of the set of rational functions drawing their numerators and denominators from the sets defined above:

[TABLE]

The set $\mathcal{R}_{T}(\delta_{D},\delta_{N})$ is finite and contains all elements of $\mathcal{L}(\mathcal{Q})$ necessary to represent the elements of $\mathcal{Q}$ with denominators from $\mathcal{D}(\delta_{D})$ and numerators from $\mathcal{N}(\delta_{N})$. Therefore, by increasing the values of $\delta_{D}$ and $\delta_{N}$, the set $\mathcal{R}_{T}(\delta_{D},\delta_{N})$ can be systematically extended to the whole set of $\mathcal{L}(\mathcal{Q})$, which contains the correct ansatz. While the correct ansatz is necessarily contained in $\mathcal{L}(\mathcal{Q})$, the choice of the finite subset $\mathcal{R}_{T}\subset\mathcal{L}(\mathcal{Q})$ presented here is a heuristic procedure. However, the knowledge of upper bounds on $\delta_{D}$ and $\delta_{N}$ would be enough to turn the algorithm into a computable criterion for the existence of a rational transformation transforming a given differential equation into canonical form.

2.4 Ansatz for the resulting canonical form

The ansatz for the resulting canonical form in Eq. (18) requires the knowledge of a set of polynomials in the invariants that encompasses the set of letters of the resulting canonical form. In this section, these letters will be shown to be a subset of the set $\mathcal{I}(a)$ of irreducible denominator factors of the original differential equation with trivial dependence on the regulator. Consider the transformation law Eq. (7)

[TABLE]

Since the transformation $T$ is rational, the derivative does not alter the set of its denominator factors $\mathcal{I}(\textup{d}T)=\mathcal{I}(T)$. Claim 1 implies $\mathcal{I}(T)\subseteq\mathcal{I}(a)$ and thus

[TABLE]

The denominator factors of $T^{-1}$ can be deduced by writing the inverse as

[TABLE]

The cofactors in the adjugate of $T$ are a sum of products of components of $T$, which implies

[TABLE]

Due to Eq. (9) and Eq. (8), the determinant of $T$ can be written in the form

[TABLE]

which leads to

[TABLE]

Thus, the irreducible denominator factors of the right-hand side of Eq. (46) have been shown to be a subset of $\mathcal{I}(a)$, hence the same holds for those of the left-hand side

[TABLE]

Since $\mathcal{I}(\textup{d}\tilde{A})$ is equal to the set of letters of $\tilde{A}$, this allows to restrict the set of polynomials in the ansatz in Eq. (18) to the set $\mathcal{I}(a)$.

2.5 Ansatz for off-diagonal blocks

The transformation $t_{D}$ in Eq. (20), which transforms the off-diagonal blocks into dlog-form, is determined by Eq. (21). A rational solution of this equation for $D$ is computed by making a rational ansatz. In this section, the set of rational functions $\mathcal{R}_{D}$ to be used for the ansatz is constructed. First, the set of irreducible denominator factors of $D$ will be shown to be a subset of the irreducible denominator factors of $b$. Here and in the following, these factors are assumed not to depend on the regulator unless stated otherwise. The argument proceeds similarly to the one in the proof of claim 1. In this case it will even be possible to derive upper bounds on the powers of the irreducible denominator factors of $D$. In a second step, these global upper bounds will be refined to upper bounds for the individual components of $D$, which reduces the number of rational functions in the ansatz considerably.

The set of possible irreducible denominator factors occurring in a rational solution $D$ of

[TABLE]

can be determined from the denominator factors of $b$. In order to demonstrate this, assume

[TABLE]

where the same notation as in Section 2.3 is used. Then, there exists a lowest order $s\in\mathbb{Z}$ in the expansion of $D$ with

[TABLE]

and therefore

[TABLE]

The derivative raises the order of $f$ by one

[TABLE]

Since $\tilde{e}$ and $\tilde{c}$ are in dlog-form and thus at most of order $1/f$, it follows

[TABLE]

which in turn implies the left-hand side and therefore also the right-hand side of Eq. (53) to be of order $1/f^{n+1}$

[TABLE]

As $b^{\prime}$ is in dlog-form and consequently at most of order $1/f$, it can be concluded

[TABLE]

This result allows to extract upper bounds on the order of the irreducible denominator factors of $D$ from a given $b$. Let $\mathcal{I}(b)=\{f_{1},\dots,f_{U}\}$ denote the set of irreducible denominator factors of $b$ and $\lambda_{i}$ the order of the denominator factor $f_{i}$

[TABLE]

According to the argument above, the upper bounds $\mu_{i}$

[TABLE]

are given by

[TABLE]

Rather than using these bounds to make an ansatz, it is beneficial to reduce the combinatorics of the ansatz by refining the above bounds. The idea is to infer bounds on the powers of the denominator factors of individual components of the solution $D$ rather than for all components at once. Assume

[TABLE]

and let $s\in\mathbb{Z}$ denote the lowest order in the expansion of $D_{ij}$ with

[TABLE]

The derivative raises the power by one

[TABLE]

Consider a component of the order $s$ in the expansion of Eq. (53)

[TABLE]

Since $b^{\prime}$ is in dlog-form, this term cannot cancel the order $1/f^{n+1}$ of the derivative term. Therefore at least one of the following cases must be true:

case 1:

[TABLE]

case 2:

[TABLE]

In case 2, there has to be at least one index $\alpha$ with either

[TABLE]

or

[TABLE]

or there exists at least one index $\beta$ with either

[TABLE]

or

[TABLE]

So far, the assumption $D^{(s)}_{ij}\sim 1/f^{n}$ has been demonstrated to imply either $b^{(s)}_{ij}\sim 1/f^{k}$ for some $k\geq n+1$ (case 1) or that some other component of $D^{(s-1)}$ is of order $1/f^{k}$ with $k\geq n$ (case 2). In case 2, the whole argument can be applied again to the respective components of $D^{(s-1)}$. This can be repeated until either the lowest order in the expansion is reached and therefore case 2 is not possible anymore or at some point only case 1 is possible due to the structure of $\tilde{e}$ and $\tilde{c}$. Thus, all possible chains of this argument necessarily end with case 1. Since $\tilde{e}$, $\tilde{c}$ and $b$ are known prior to the computation of $D$, the chains can be followed backwards in order to derive upper bounds on the powers of the denominator factors of the components of $D$. The idea is to consider all chains at once and start at the last step by reversing case 1 for all components of $D$. Using the powers of the denominator factors of $b$

[TABLE]

case 1 is reversed for all components by setting the upper bounds $\mu_{k,ij}$ of $D_{ij}$ on $f_{k}$

[TABLE]

to

[TABLE]

It can then be deduced from $\tilde{e}$ and $\tilde{c}$ for each component which other components could have implied the current bounds via case 2. For instance, if there exists an $\alpha$ with

[TABLE]

and the current bound for the order in $1/f$ of $D_{\alpha j}$ is $n$, case 2 is reversed by setting the bound on the order of $D_{ij}$ to $n$ as well, unless it is already higher. At each step it is checked for all components $D_{ij}$, whether there is an $\tilde{e}_{i\alpha}$ as in Eq. (70) or Eq. (71) or a $\tilde{c}_{\beta j}$ as in Eq. (72) or Eq. (73). If this is the case, the bounds are updated accordingly. This is repeated until the bounds do not change anymore, and therefore they incorporate all possible cases. Algorithm 1 summarizes the procedure.

Since the values of the bounds $\mu_{k,ij}$ can only increase during each iteration in algorithm 1 and the overall bounds $\mu_{k}$ given in Eq. (63) are upper bounds on the bounds of the components

[TABLE]

it is clear that the algorithm terminates after a finite number of steps. Using the bounds computed with Algorithm 1, the following sets of rational functions can be defined:

[TABLE]

with

[TABLE]

The above argument shows that for high enough $\delta_{N}$, the component $D_{ij}$ is an element of the $K$-span of $\mathcal{R}_{ij}(\delta_{N})$. For the ansatz, a basis of Leinartas functions of the $K$-span of the union of all $\mathcal{R}_{ij}(\delta_{N})$ is taken

[TABLE]

It would be more efficient to make a different ansatz for each component of $D$ using $\mathcal{L}(\mathcal{R}_{ij}(\delta_{N}))$. However, this functionality will only be included in a future version of CANONICA.

2.6 On the uniqueness of canonical bases

The application of a constant invertible transformation $C$ to a differential equation in canonical form obviously preserves the canonical form

[TABLE]

This raises the question whether all canonical forms can be obtained in this way. The following claim shows that indeed every canonical form can be obtained by a constant transformation from any other canonical form. In this sense the canonical form of a given differential equation is unique up to constant transformations.

Claim 2**.**

Let $a(\epsilon,\{x_{j}\})$ be a differential equation of Feynman integrals and $T_{1}(\epsilon,\{x_{j}\})$ and $T_{2}(\epsilon,\{x_{j}\})$ be invertible rational transformations, which transform it into the canonical forms $\epsilon\textup{d}\tilde{A}_{1}(\{x_{j}\})$ and $\epsilon\textup{d}\tilde{A}_{2}(\{x_{j}\})$ respectively. Then there exists a constant invertible transformation $C$ transforming $\epsilon\textup{d}\tilde{A}_{1}(\{x_{j}\})$ into $\epsilon\textup{d}\tilde{A}_{2}(\{x_{j}\})$.

Consider the transformation $T=T^{-1}_{1}T_{2}$, which transforms $\epsilon\textup{d}\tilde{A}_{1}$ to $\epsilon\textup{d}\tilde{A}_{2}$. First, the transformation $T$ has to be shown to be independent of the invariants. The corresponding transformation law reads

[TABLE]

It is instructive to rewrite this equation:

[TABLE]

The summation over the letters is meant to run over the union of the sets of letters of the two canonical forms, since it is a priori not clear that they both have exactly the same set of letters. The letters are assumed to be irreducible polynomials and the union is meant to remove all scalar multiples of letters as well. Since the transformation law is invariant under the multiplication of $T$ with any rational function $g(\epsilon)$, the $\epsilon$-expansion of $T$ can be assumed to start at the order $\epsilon^{0}$. Then the first order in the expansion of the above equation reads

[TABLE]

and therefore $T^{(0)}$ has to be constant. At any order $n>0$ the expansion of the above equation is given by

[TABLE]

Assuming $T^{(n-1)}$ to be constant, this equation can easily be integrated

[TABLE]

Since $T_{1}$ and $T_{2}$ are assumed to be rational in $\epsilon$ and the invariants, the same holds for $T$ and therefore the coefficients of its $\epsilon$-expansion have to be rational as well. This implies

[TABLE]

and consequently $T^{(n)}$ has to be constant. By induction, these arguments imply that all coefficients of the $\epsilon$-expansion of $T$ are constant and therefore $T=T(\epsilon)$. As $T$ is independent of the invariants, the transformation law Eq. (83) has the form

[TABLE]

It can be concluded that $T(\epsilon)$ transforms $\textup{d}\tilde{A}_{1}$ to $\textup{d}\tilde{A}_{2}$ for all non-singular values of $\epsilon$, because the left-hand side does not depend on $\epsilon$. Upon choosing such a value $\epsilon_{0}$, a constant invertible transformation $C=T(\epsilon_{0})$ is obtained, which concludes the proof of the claim. The same argument also holds for the more general case of an algebraic dependence of $T_{1}$ and $T_{2}$ on $\epsilon$ and the invariants. Altogether, canonical forms have been shown to be unique modulo $GL(m,K)$ transformations.

This result explains the origin of the non-linear parameter equations, which are treated in Section 2.7. Moreover, the above result can be utilized for the comparison of two different canonical forms of the same problem, provided they are expressed in the same set of invariants. In this situation, claim 2 asserts the existence of a constant transformation relating the two canonical forms. This can be tested by checking whether the following system of linear equations

[TABLE]

has a non-singular solution for the components of $C$.

The uniqueness of the canonical form also manifests itself in Eq. (21), which governs the transformation of the off-diagonal blocks into dlog-form. In the following, the rational solution of this equation is proven to be unique up to the addition of terms depending solely on the regulator. In practice, this result allows to exclude terms with trivial dependence on the invariants from the ansatz without losing generality.

For a given $b$, let $D$ and $b^{\prime}$ satisfy Eq. (21)

[TABLE]

with $b^{\prime}$ understood to be in dlog-form. Adding a term $C(\epsilon)$ to the solution $\tilde{D}=D+C(\epsilon)$ solves the same equation

[TABLE]

with

[TABLE]

which is also in dlog-form, since $\tilde{e}$ and $\tilde{c}$ are in dlog-form. This argument establishes the freedom to add terms independent of the invariants to a solution of Eq. (21). The following argument proves this to be the only possible relation between two solutions of Eq. (92). Let $D_{1}$ and $D_{2}$ satisfy Eq. (92) for a given $b$

[TABLE]

Then the difference $\bar{D}=D_{1}-D_{2}$ satisfies

[TABLE]

Let $\hat{\bar{D}}=\bar{D}\epsilon^{\tau}$ be defined such that the expansion of $\hat{\bar{D}}$ starts at the constant order. The equation for $\hat{\bar{D}}$ reads

[TABLE]

with $B=\epsilon^{\tau}(b_{2}^{\prime}-b_{1}^{\prime})$, which is in dlog-form. The first order in the expansion of Eq. (98) reads

[TABLE]

which integrates to

[TABLE]

As $D_{1}$ and $D_{2}$ are assumed to be rational, $\hat{\bar{D}}$ has to be rational as well and therefore

[TABLE]

which implies that $\hat{\bar{D}}^{(0)}$ is constant. Consider the expansion of Eq. (98) at some order $n>0$

[TABLE]

The right-hand side is in dlog-form for constant $\hat{\bar{D}}^{(n-1)}$ and therefore $\hat{\bar{D}}^{(n)}$ can only be rational if it is constant as well. This proves by induction that $\hat{\bar{D}}$ is independent of the invariants. Consequently, the difference of the solutions $\bar{D}=\hat{\bar{D}}\epsilon^{-\tau}$ has to be independent of the invariants as well. Altogether, the argument establishes the uniqueness of a rational solution for $D$ of Eq. (21) up the addition of terms that are independent of the invariants. This fact can be used in practice to fix this freedom without losing generality.

2.7 Treatment of non-linear parameter equations

In the course of applying the algorithm, the ansatz in Eq. (16) and Eq. (18) is inserted in the expansion of Eq. (14). By requiring the resulting equations to hold for all allowed values of the invariants, a system of equations in the unknown parameters is obtained at each order of the expansion. Due to the term $\epsilon T\textup{d}\tilde{A}$ in Eq. (11), these equations can be non-linear. Instead of directly solving these non-linear equations, it will be shown in the following how they can be reduced to linear equations by imposing appropriate constraints.

In Section 2.6, it has been shown that the resulting canonical form is uniquely fixed up to an invertible constant transformation. Exactly this ambiguity leads to the non-linear equations, because if $\textup{d}\tilde{A}$ was fixed, the term $\epsilon T\textup{d}\tilde{A}$ would not generate non-linear equations. Therefore, the non-linear equations can be turned into linear equations by fixing the degrees of freedom in the ansatz corresponding to a subsequent invertible constant transformation. In order to fix these degrees of freedom directly, they would have to be disentangled from those which are determined by the equations in the parameters. Since this would require a parameterization of the solution set of the non-linear equations, which is essentially equivalent to solving them, a more indirect approach of fixing the freedom is taken in the following.

To this end, suppose for the moment that the parameters of the ansatz can be separated in those which are fixed by the parameter equations $\{\tau\}$ and those which correspond to the remaining freedom $\{\tau^{\prime}\}$. Let $T(\epsilon,\{x_{j}\},\{\tau\},\{\tau^{\prime}\})$ be a solution of Eq. (11), provided the parameters $\{\tau\}$ solve the parameter equations. According to the proof of claim 2, this transformation can be thought of as the product of some fixed transformation $T_{1}(\epsilon,\{x_{j}\},\{\tau\})$, which transforms the original differential equation to some canonical form $\epsilon\textup{d}\tilde{A}_{1}$, and a transformation $C(\epsilon,\{\tau^{\prime}\})$ parameterizing the transformation of $\textup{d}\tilde{A}_{1}$ to any other possible canonical form $\epsilon\textup{d}\tilde{A}_{2}(\tau^{\prime})$

[TABLE]

[TABLE]

It should be noted that $\tilde{A}_{2}$ does in general only depend on a subset of the parameters $\{\tau^{\prime}\}$, because some parameters can correspond to a non-trivial $\epsilon$-dependence of $C(\epsilon,\{\tau^{\prime}\})$. As mentioned above, the goal is to fix the resulting differential equation $\textup{d}\tilde{A}_{2}$ by fixing the corresponding parameters of $\{\tau^{\prime}\}$. This can be achieved by demanding $C(\epsilon,\{\tau^{\prime}\})$ to equal some fixed constant invertible transformation at some non-singular value $\epsilon=\epsilon_{0}$. Since the left-hand side of Eq. (104) does not depend on $\epsilon$, this completely fixes $\textup{d}\tilde{A}_{2}$ irrespective of the particular value of $\epsilon_{0}$. However, fixing $C(\epsilon,\{\tau^{\prime}\})$ directly would require the computation of the factorization in Eq. (103), which is only possible if the separation of the parameters into the sets $\{\tau\}$ and $\{\tau^{\prime}\}$ is known. Instead, $C(\epsilon,\{\tau^{\prime}\})$ can be fixed indirectly by demanding

[TABLE]

to hold at some non-singular point $\{x\}=\{x_{j0}\}$, $\epsilon=\epsilon_{0}$. This is equivalent to fixing $C(\epsilon,\{\tau^{\prime}\})$ as follows:

[TABLE]

The constraints given by Eq. (105) can be imposed without being able to separate the parameters into $\{\tau\}$ and $\{\tau^{\prime}\}$. Moreover, these constraints are linear in both the $\{\tau\}$ and the $\{\tau^{\prime}\}$, since the ansatz in Eq. (16) is linear in all parameters. Therefore, the additional constraints in Eq. (105) can be used to completely fix the resulting canonical form, which turns the non-linear parameter equations into linear equations.

Recall that the parameter equations are generated order by order in the expansion of Eq. (14) and at each order it is tested whether the series terminates at the current order. The constraints in Eq. (105) can only be imposed if the full $T(\epsilon,\{x_{j}\},\{\tau\},\{\tau^{\prime}\})$ is known. Thus, the computation must have reached the order at which the series terminates. However, non-linear equations can already occur at lower orders in the expansion, i.e. before Eq. (105) can be imposed to turn them into linear equations. The strategy described in the following overcomes this point by essentially just solving the linear equations at each order and keeping the non-linear equations until the constraints can be imposed.

At each order in the expansion of the transformation, the linear equations are solved first and then their solution is inserted into the non-linear ones, which possibly turns some of them into linear equations. These newly generated linear equations can again be solved. This procedure is iterated until no further linear equations are generated. The remaining non-linear equations are kept unsolved. It is then tested whether the series terminates at the current order by generating the additional parameter equations corresponding to this assumption (cf. Ref. [72]). The linear equations of these additional equations are then iteratively solved as described above, while the previously obtained still unsolved non-linear equations are taken into account as well. If it turns out during this iteration that the system has no solution, the algorithm proceeds with the next order in the expansion of the transformation. If this is not the case and some non-linear equations remain at the end of the iteration, Eq. (105) is imposed. If the series does terminate at the current order, the additional constraints will turn the remaining non-linear equations into linear ones, which then determine the transformation. If either non-linear equations remain or the linear ones have no solution, it can be concluded that the series does not terminate at the current order and the algorithm proceeds with the next order in the expansion.

Altogether, this procedure allows to compute a transformation to a canonical form by only solving linear equations at each order without sacrificing the generality of the algorithm.

3 The CANONICA package

This section introduces the CANONICA package, which implements the algorithm proposed in Ref. [72]. After describing the installation and the contents of the package, the main functionality is illustrated with short usage examples. A detailed description of all functions and options of CANONICA can be found in the interactive manual ./manual.nb.

3.1 Installation

CANONICA is a Mathematica package and requires an installation of version 10 or higher of Mathematica. The CANONICA repository can be copied to the local directory with

git clone https://github.com/christophmeyer/CANONICA.git

Alternatively, an archive file can be downloaded at

https://github.com/christophmeyer/CANONICA/archive/v1.0.tar.gz

which may be extracted with

tar -xvzf CANONICA-1.0.tar.gz

There is no further installation necessary, in particular, there are no dependencies other than Mathematica. In a Mathematica session, the package can be loaded by

Get["CANONICA.m"];

provided the file CANONICA.m is placed either in the current working directory or in one of the search paths. If this is not the case, Get either has to be called with the full path of the file CANONICA.m, or its location has to be added to the list of Mathematica’s search paths, which is stored in the global variable $Path, by e.g.

AppendTo[$Path,"/path/to/CANONICA/src/"]

Changes to $Path can be made permanent by adding them to the initialization file init.m.

3.2 Files of the package

The root directory of the CANONICA package contains the following files and directories.

./src/CANONICA.m

Contains all of the source code of the program, in particular, all function definitions as well as short usage messages for the public functions and options.

./manual.nb

An interactive manual in the Mathematica notebook format explaining the usage of all functions and options with short examples.

./examples

Several examples are provided in this directory. The directory of each example contains a .m file with the corresponding differential equation and a .nb notebook file illustrating the application of CANONICA to this example. The calculation of the full transformation can also be run in terminal mode with the script RunExample.m. The script is started by calling

math -run "<<RunExample.m"

or

math -script RunExample.m

Some basic information about the examples, such as the master integrals and the definition of the kinematic invariants, is provided in the ./examples/examples.pdf file.

./LICENSE

A copy of the third version of the GNU General Public License.

./README

A README file providing basic information on the package.

3.3 Usage examples

In this section, the main features and their usage are illustrated with short examples. A similar, but more extensive account of this functionality can be found in the manual notebook file, which comes with the package.

The most common input required by CANONICA is a differential equation of the form in Eq. (1), which is determined by the differential form $a(\epsilon,\{x_{j}\})$. Consider the following example:

[TABLE]

depending on the invariants $x$ and $y$. The differential form is represented in CANONICA as a list of the matrix-valued coefficients of the differentials of the invariants. The dimensional regulator $\epsilon$ has to be denoted by the protected symbol eps. For the above example, the input reads

a = { {{-(2+eps)/x, 0}, {0, -(1+eps)/x}} , {{0, 0}, {((-1+eps)x)/((-1+y)y), (1-eps(1+y))/((-1+y)y)}} };

The order of the coefficient matrices has to be specified by a list of the corresponding invariants.

invariants = {x, y};

The algorithm to compute a transformation of a differential equation to canonical form as outlined in Section 2.2 and in Ref. [72] is implemented in the function TransformDiagonalBlock. As the name suggests, this function is intended to be used for calculating the transformation of diagonal blocks of differential equations. For the example above, the function is called as follows:

res=TransformDiagonalBlock[a, invariants]

which returns the output

  {
   {{(1-2eps)/x^2, (1-2eps)/x^2}, {(1-eps)/x, (1-eps)/(xy)}}
   ,
   {
    {{-(eps/x), 0}, {0, -(eps/x)}}
    ,
    {{-eps/(-1+y),-eps/(y-y^2)}, {eps/(-1+y),eps/(y-y^2)}}
   }
  }

The output of TransformDiagonalBlock is a list with two entries. The first contains the transformation and the second contains the resulting differential equation in $\epsilon$-form. The resulting differential equation is of course redundant, since it can be computed by applying the transformation to the original differential equation. However, the resulting differential equation is generated anyway in the course of the computation of the transformation and applying the transformation can be a costly operation in itself for larger differential equations.

The application of a transformation to a differential equation, according to the transformation law Eq. (5), is implemented in the function TransformDE, which for some transformation

trafo = res[[1]]

is called as

TransformDE[a, invariants, trafo]

and returns the resulting differential equation. In order to apply Eq. (5), the inverse of the transformation needs to be computed, which can consume significant computation time for larger matrices, when done with the build-in Mathematica command. However, the transformations usually exhibit a block-triangular structure, which is exploited by TransformDE leading to a considerably better performance.

The function TransformDiagonalBlock is in principle applicable to differential equations of any size. However, the performance can be improved significantly by splitting the computation according to the block-triangular structure of the differential equation and performing the computation in a recursion over the sectors of the differential equation. The main function in CANONICA for this purpose is RecursivelyTransformSectors. In addition to the two arguments related to the differential equation itself, this function expects an argument that defines the boundaries of the diagonal blocks. The differential equation in the example above actually splits into two blocks of dimension one, and in this case the boundaries read

boundaries = {{1, 1}, {2, 2}};

Each entry of the boundaries list corresponds to one diagonal block, which is specified by the position of its lowest and highest integral. Instead of using TransformDiagonalBlock to transform a into $\epsilon$-form all at once, the following command

RecursivelyTransformSectors[a, invariants, boundaries, {1, 2}]

computes the transformation in a recursion over the sectors, as described in Section 2.2. The last argument determines the sectors at which the computation starts and ends. The output is of the same format as described above for TransformDiagonalBlock. If some lower sectors have already been transformed into $\epsilon$-form and the computation should therefore not start at the first sector, the differential equation of the lower-sectors in $\epsilon$-form and the transformation leading to it need to be provided as two additional arguments.

CANONICA also has functionality to extract the boundaries of the diagonal blocks from the differential equation. The function SectorBoundariesFromDE extracts the most fine grained boundaries compatible with the differential equation. For instance, in the example above

SectorBoundariesFromDE[a]

returns

{{1, 1}, {2, 2}}

The boundaries obtained in this way may be too fine for the algorithm to find the solution, since the solution space could be constrained too much by splitting the transformation into smaller blocks. It is safer to choose the boundaries according to the sector-ids of the integrals, which in general yields coarser grained boundaries. For a given list of integrals specified by their propagator powers

masterIntegrals={Int["T1", {0, 1, 0, 1, 0, 1, 0, 0, 0}], Int["T1", {0, 1, 0, 1, 1, 1, 0, 0, 0}], Int["T1", {0, 0, 1, 1, 1, 1, 0, 0, 0}], Int["T1", {0, 1, 1, 1, 1, 1, 0, 0, 0}], Int["T1", {1, 1, 0, 0, 0, 0, 1, 0, 0}], Int["T1", {1, 1, -1, 0, 0, 0, 1, 0, 0}]};

the boundaries for the corresponding differential equation can be computed with

SectorBoundariesFromID[masterIntegrals]

provided the integrals are ordered with respect to their sector-id. SectorBoundariesFromID then returns the sector boundaries derived from the sector-ids of the integrals in the above format

{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 6}}

While the main function of CANONICA is RecursivelyTransformSectors, it is in some cases useful to be able to perform only certain steps of the algorithm. For this reason, there is a hierarchy of functions available in CANONICA allowing to break the calculation of the transformation into smaller steps. The hierarchy of these lower-level functions is illustrated in Fig. 1. For more information on specific functions, see the manual notebook included in the package.

3.4 Tests and Limitations

CANONICA has been successfully tested on a variety of non-trivial single- and multi-scale problems, some of which are included as examples in the package. All tests have been performed with the Mathematica versions 10 and 11 on a Linux operating system. The limitations of CANONICA are mostly limitations of the algorithm itself. In particular, the algorithm is limited to differential equations for which a rational transformation to a canonical form exists. However, it is well known that rational differential equations may require non-rational transformations to attain a canonical form. The following example illustrates this behavior:

[TABLE]

where the transformation to a canonical form is given by

[TABLE]

In this situation it is often possible to render the transformation rational with a change of coordinates. For instance, in the above example, the differential form transforms under the change of variables

[TABLE]

to

[TABLE]

which has the following rational transformation to a canonical form

[TABLE]

While a change of coordinates can remove non-rational letters in more complicated examples as well [6, 9, 11, 12], this has neither been proven to be always possible, nor is a general method to construct such coordinate changes known. In fact, the existence of such a procedure appears to be unlikely, given that the number of independent roots can largely outgrow the number of variables in a problem [32].

In addition to this limitation of the algorithm, the calculations are in practice limited by the size of the available memory. This imposes limits on the size of the systems of linear parameter equations. The main factors determining the sizes of these systems of linear equations are the size of the differential equation itself and the size of the ansatz. Thus, run time and memory consumption of CANONICA are highly problem dependent. For instance, the most complicated example that is provided with the package is a two-loop double box topology depending on three dimensionless scales. It has a run time of about 20 minutes and a memory consumption of less than 8 GB.

4 Conclusion

The description of the algorithm in Ref. [72] has been extended in the present paper by providing details on the choice of an ansatz for both the diagonal blocks and the off-diagonal blocks. In both cases it has been shown that the irreducible denominator factors can be chosen from the set of irreducible denominator factors occurring in the differential equation. The ansatz for the off-diagonal blocks has been further restricted by proving upper bounds on the powers of the denominator factors in the solution.

Furthermore, canonical forms have been proven to be unique up to constant transformations, which allowed to attribute the occurrence of non-linear equations in the parameters to the freedom in the choice of this constant transformation. By fixing this freedom in a specific way, it has been argued that only linear equations need to be solved in the course of applying the algorithm.

The main focus of this publication has been the presentation of the Mathematica package CANONICA, which implements the aforementioned algorithm and allows to compute rational transformations of differential equations into canonical form. This represents the first publicly available implementation of an algorithm applicable to problems depending on multiple scales. CANONICA has been successfully tested on a number of state of the art multi-scale problems, including previously unknown integral topologies. CANONICA may thus provide a valuable contribution to the ongoing efforts to automatize multi-loop calculations.

Acknowledgments

The author would like to thank Peter Uwer for useful discussions and comments on the manuscript. This research was supported by the German Research Foundation (DFG) via the Research Training Group 1504.

Appendix A List of functions provided by CANONICA

CalculateDlogForm:

CalculateDlogForm[a, invariants, alphabet] returns a list of matrices of the same dimensions as a, where each matrix is the matrix-residue of one of the letters. The ordering is the same as the one in alphabet. Returns False if a cannot be cast in a dlog-form with the given alphabet.

CalculateNexta:

CalculateNexta[aFull, invariants, sectorBoundaries, trafoPrevious, aPrevious] applies trafoPrevious to aFull and returns the differential equation of the next sector. aPrevious is used to recycle the transformation of lower sectors.

CalculateNextSubsectorD:

CalculateNextSubsectorD[a, invariants, sectorBoundaries, previousD] computes the $D_{k}$ of the next sector, prepends it to previousD and returns the result. The ansatz to be used can be specified with the optional argument userProvidedAnsatz. If no ansatz is provided, an ansatz is generated automatically. The size of the automatically generated ansatz can be controlled with the option DDeltaNumeratorDegree.

CheckDlogForm:

CheckDlogForm[a, invariants, alphabet] tests whether the differential equation a is in dlog-form for the given alphabet. Returns either True or False.

CheckEpsForm:

CheckEpsForm[a, invariants, alphabet] tests whether the differential equation a is in $\epsilon$-form with the given alphabet. Returns either True or False.

CheckIntegrability:

CheckIntegrability[a, invariants] tests whether a satisfies the integrability condition $\textup{d}a-a\wedge a=0$ and returns either True or False.

CheckSectorBoundaries:

CheckSectorBoundaries[a, sectorBoundaries] tests whether the sectorBoundaries are compatible with a and returns either True or False.

ExtractDiagonalBlock:

ExtractDiagonalBlock[a, boundaries] returns the diagonal block of the differential equation a specified by the boundaries argument. boundaries is expected to be of the format {nLowest, nHighest}, where nLowest and nHighest are positive integers indicating the lowest and highest integrals of the diagonal block, respectively.

ExtractIrreducibles:

ExtractIrreducibles[a] returns the irreducible denominator factors of a that do not depend on the regulator. The option AllowEpsDependence->True allows the irreducible factors to depend on both the invariants and the regulator.

FindAnsatzSubsectorD:

FindAnsatzSubsectorD[a, invariants, sectorBoundaries, previousD] takes a differential equation a, which is required to be in $\epsilon$-form except for the off-diagonal block of the highest sector. Needs to be provided with all previous $D_{k}$ in the argument previousD and computes the ansatz $\mathcal{R}_{D}$ for the computation of the next $D_{k}$. Takes the option DDeltaNumeratorDegree to enlarge the ansatz. For more details, see Section 2.5.

FindAnsatzT:

FindAnsatzT[a, invariants] takes a differential equation a in the invariants and computes an ansatz $\mathcal{R}_{T}$ as described in Section 2.3. The ansatz can be enlarged with the options TDeltaNumeratorDegree and TDeltaDenominatorDegree.

FindConstantNormalization:

FindConstantNormalization[invariants, trafoPrevious, aPrevious] calculates a constant diagonal transformation to minimize the number of prime factors present in the matrix-residues. The transformation is composed with trafoPrevious and returned together with the resulting differential equation.

FindEpsDependentNormalization:

FindEpsDependentNormalization[a, invariants] calculates a diagonal transformation depending only on the dimensional regulator in order to attempt to minimize the number of orders that need to be calculated in a subsequent determination of the transformation to a canonical form. Returns the transformation together with the resulting differential equation.

RecursivelyTransformSectors:

RecursivelyTransformSectors[aFull, invariants, sectorBoundaries, {nSecStart, nSecStop}] calculates a rational transformation of aFull to a canonical form in a recursion over the sectors of the differential equation, which have to be specified by sectorBoundaries. The arguments nSecStart and nSecStop set the first and the last sector to be computed, respectively. If nSecStart is greater than one, the result of the calculation for the sectors lower than nSecStart needs to be provided in the additional arguments trafoPrevious and aPrevious. RecursivelyTransformSectors returns the transformation of aFull to a canonical form for the sectors up to nSecStop and the resulting differential equation. The ansaetze for the individual blocks are generated automatically. The sizes of the ansaetze for the diagonal blocks can be controlled with the options TDeltaNumeratorDegree and TDeltaDenominatorDegree. Similarly, the sizes of the ansaetze for the off-diagonal blocks are controlled by the option DDeltaNumeratorDegree.

SectorBoundariesFromDE:

SectorBoundariesFromDE[a] returns the most fine grained sector boundaries compatible with a.

SectorBoundariesFromID:

SectorBoundariesFromID[masterIntegrals] takes a list of masterIntegrals, which need to be ordered by their sector-id and returns the sector boundaries computed from the sector-ids.

TransformDE:

TransformDE[a, invariants, t] applies the transformation t to the differential equation a. Returns $a^{\prime}=t^{-1}at-t^{-1}\textup{d}t$. The option SimplifyResult->False deactivates the simplification of the result.

TransformDiagonalBlock:

TransformDiagonalBlock[a, invariants] calculates a rational transformation to transform a into canonical form and returns the transformation together with the resulting differential equation. With the optional argument userProvidedAnsatz, the user can specify the ansatz to be used. If no ansatz is provided, an ansatz is generated automatically. The size of the automatically generated ansatz can be controlled with the options TDeltaNumeratorDegree and TDeltaDenominatorDegree.

TransformDlogToEpsForm:

TransformDlogToEpsForm[invariants, sectorBoundaries, trafoPrevious, aPrevious] computes a transformation depending only on the regulator in order to transform aPrevious from dlog-form into canonical form (cf. Ref. [70]). The transformation is composed with trafoPrevious and returned together with the resulting differential equation. Per default, the transformation is demanded to be in a block-triangular form induced by sectorBoundaries. This condition can be dropped with the option EnforceBlockTriangular->False.

TransformNextDiagonalBlock:

TransformNextDiagonalBlock[aFull, invariants, sectorBoundaries, trafoPrevious, aPrevious] calls TransformDiagonalBlock to compute the transformation of the next diagonal block into canonical form and composes it with trafoPrevious. Returns the composed transformation together with the resulting differential equation. With the optional argument userProvidedAnsatz, the user can specify the ansatz to be used. If no ansatz is provided, an ansatz is generated automatically. The size of the automatically generated ansatz can be controlled with the options TDeltaNumeratorDegree and TDeltaDenominatorDegree.

TransformNextSector:

TransformNextSector[aFull, invariants, sectorBoundaries, trafoPrevious, aPrevious] transforms the next sector into canonical form, composes the calculated transformation with trafoPrevious and returns it together with the resulting differential equation. With the optional argument userProvidedAnsatz, the user can specify the ansatz to be used for the diagonal block. If no ansatz is provided, an ansatz is generated automatically. The size of the automatically generated ansatz for the diagonal block can be controlled with the options TDeltaNumeratorDegree and TDeltaDenominatorDegree. Similarly, the sizes of the ansaetze for the off-diagonal blocks are controlled by the option DDeltaNumeratorDegree.

TransformOffDiagonalBlock:

TransformOffDiagonalBlock[invariants, sectorBoundaries, trafoPrevious, aPrevious] assumes aPrevious to be in canonical form except for the highest sector of which only the diagonal block is assumed to be in canonical form. Computes a transformation to transform the off-diagonal block of the highest sector into dlog-form. This transformation is composed with trafoPrevious and returned together with the resulting differential equation. Proceeds in a recursion over sectors, which can be resumed at an intermediate step by providing all previous $D_{k}$ (cf. Ref. [72]) in the optional argument userProvidedD. The sizes of the automatically generated ansaetze for the off-diagonal blocks are controlled by the option DDeltaNumeratorDegree.

Appendix B List of options

AllowEpsDependence:

AllowEpsDependence is an option of ExtractIrreducibles controlling whether irreducible factors depending on both the invariants and the regulator are returned as well. The default value is False.

DDeltaNumeratorDegree:

DDeltaNumeratorDegree is an option controlling the numerator powers of the rational functions in the ansatz used for the computation of $D$ for the transformation of off-diagonal blocks. The default value is 0. For more details, see Section 2.5. DDeltaNumeratorDegree is an option of the following functions: CalculateNextSubsectorD, FindAnsatzSubsectorD, RecursivelyTransformSectors, TransformNextSector, TransformOffDiagonalBlock.

EnforceBlockTriangular:

EnforceBlockTriangular is an option of TransformDlogToEpsForm controlling whether the resulting transformation is demanded to be in the block-triangular form induced by the sectorBoundaries argument. The default value is True.

FinalConstantNormalization:

FinalConstantNormalization is an option of RecursivelyTransformSectors controlling whether FindConstantNormalization is invoked after all sectors have been transformed into canonical form in order to simplify the resulting canonical form. The default value is False.

PreRescale:

PreRescale is an option of TransformDiagonalBlock controlling whether FindEpsDependentNormalization is called prior to the main computation in order to attempt to minimize the number of orders that need to be calculated in a subsequent determination of the transformation to a canonical form. The default value is True.

SimplifyResult:

SimplifyResult is an option of TransformDE controlling whether the resulting differential equation is simplified. The default value is True.

TDeltaDenominatorDegree:

TDeltaDenominatorDegree is an option controlling the denominator powers of the rational functions in the ansatz used for the computation of the transformation of diagonal blocks. The default value is 0. For more details, see Section 2.3. TDeltaDenominatorDegree is an option of the following functions: FindAnsatzT, RecursivelyTransformSectors, TransformNextDiagonalBlock, TransformNextSector.

TDeltaNumeratorDegree:

TDeltaNumeratorDegree is an option controlling the numerator powers of the rational functions in the ansatz used for the computation of the transformation of diagonal blocks. The default value is 0. For more details, see Section 2.3. TDeltaNumeratorDegree is an option of the following functions: FindAnsatzT, RecursivelyTransformSectors, TransformNextDiagonalBlock, TransformNextSector.

VerbosityLevel:

VerbosityLevel is an option controlling the verbosity of several main functions. Takes integer values from 0 to 12 with a value of 12 resulting in the most detailed output about the current state of the computation and a value of 0 suppressing all output but warnings about inconsistent inputs. The default value is 10. The following functions accept the VerbosityLevel option: CalculateNextSubsectorD, FindConstantNormalization, RecursivelyTransformSectors, TransformDiagonalBlock, TransformDlogToEpsForm, TransformNextDiagonalBlock, TransformNextSector, TransformOffDiagonalBlock.

Appendix C List of global variables and protected symbols

$ComputeParallel:

$ComputeParallel is a global variable that needs to be set to True to enable parallel computations. The number of kernels to be used is controlled by$NParallelKernels.

$NParallelKernels:

$NParallelKernels is a global variable setting the number of parallel kernels to be used.$NParallelKernels has no effect if $ComputeParallel is not set to True. If$ComputeParallel is True and $NParallelKernels is not assigned a value, then all available kernels are used for the computation.

eps:

eps is a protected symbol representing the dimensional regulator.