Deriving canonical differential equations for Feynman integrals from a single uniform weight integral

Differential equations are a powerful tool for evaluating Feynman integrals. Their solution is straightforward if a transformation to a canonical form is found. In this paper, we present an algorithm for finding such a transformation. This novel technique is based on a method due to Höschele et al. and relies only on the knowledge of a single integral of uniform transcendental weight. As a corollary, the algorithm can also be used to test the uniform transcendentality of a given integral. We discuss the application to several cutting-edge examples, including non-planar four-loop HQET and non-planar two-loop five-point integrals. A Mathematica implementation of our algorithm is made available together with this paper.


JHEP05(2020)025
1 Introduction Feynman integrals are ubiquitous in perturbative quantum field theory. They are required in order to extract predictions from the theory beyond the leading perturbative order. As a consequence, they are important in many areas. A prominent example is collider physics, where the underlying scattering processes are described by Feynman diagrams, and consequently on-shell momentum space Feynman integrals are needed. Another example are off-shell position-space correlation functions (e.g. in conformal field theory), from which one can determine the scaling dimension of fields, or renormalization group coefficients. Beyond their obvious importance for physics, Feynman integrals also have interesting connections to mathematics. The reason is that Feynman integrals are periods [1,2], and give rise to interesting classes of special functions. Moreover, they can be studied using partial differential equations, and algebraic geometry plays an important role in the analysis of the integrands.
Given these motivations, it is not surprising that this has been an active area of research for decades, and continues to be. It has produced numerous insights and new methods for Feynman integrals, as reviewed in [3][4][5][6]. One important concept that has emerged, initially in the maximally supersymmetric theory N = 4 super Yang-Mills, is that functions of uniform transcendental (UT) weight play a special role. In the case of JHEP05(2020)025 multiple polylogarithms, the weight refers to the (minimal) number of integrations that are needed to obtain the function, starting from a rational function. For example a logarithm has weight one, and a dilogarithm has weight two.
The weight property serves as an important organizing principle, but also significantly simplifies calculations. It was observed that uniform weight integrals satisfy particularly simple, canonical differential equations. One of the present authors conjectured that one may always chose a basis of Feynman integrals such that their differential equations take a canonical form [7]. Indeed, it is known that any Feynman integral satisfies some n-th order (Picard-Fuchs) differential equation. Equivalently, the latter can be rewritten as a coupled n × n first-order system of differential equations. For example, let us assume one kinematic variable x and take the dimension to be D = 4 − 2 , then we have In general for Feynman integrals (but under certain conditions on the basis vector f ), the matrix A(x, ) is rational. However, in general the dependence is very complicated, and hence solving such equations is a difficult task. The statement of [7] is that a basis transformation exists that significantly simplifies the matrix to a form where the solution of the DE can be read off in terms of special functions. Perhaps equally important to this guiding principle is the insight which Feynman integrals evaluate to UT functions [8], and hence provide a suitable basis for the differential equations [7]. The key idea is that the loop integrand contains the necessary information, and that the latter can be extracted systematically by taking (multi-dimensional) residues. If no doubles poles are encountered along the way, and the maximal residues (that localize all integrations, and are called leading singularities) are constant, the corresponding Feynman integral is conjectured to be UT [8]. See refs. [9,10] for further details.
Subsequently, a complementary approach [4,11] was proposed independently by one of the present authors, and by Lee. It consists in systematically constructing basis transformations that simplify A(x, ) in eq. (1.1). The first step focuses on simplifying the x-dependence, trying to make the expected singularity structure of the Feynman integrals manifest. The mathematics behind this has its roots in work by Moser [12], see also [13]. In a second step, one simplifies the dependence. A complete algorithm was first given by Lee [11], and several computer programs implementing this and related ideas exist [14][15][16].
However, both methods have certain limitations. In the leading singularity approach, it is usually simple to find some UT integrals, but obtaining a complete set of basis integrals is more difficult. Also, sometimes it is necessary to analyze leading singularities beyond four dimensions [17]. Similarly, in its current implementation the Moser approach works well only for small matrices, and for few kinematic variables.
Here we present a different method that generalizes previous work by Höschele et al. [18]. The starting point is the knowledge that one integral is UT. The key insight is that this assertion contains a lot of information (in some sense infinite amount of information). We can use the latter to find, algorithmically, a complete UT basis, and hence obtain canonical differential equations.

JHEP05(2020)025
The reason this works can be understood as follows. The statement that an integral is UT in dimensional regularization means that, at a given order k in its Laurent expansion, the coefficient function has uniform transcendental weight k (given some choice for the overall normalization). In other words, we obtain a non-trivial constraint at each order in . On the other hand, a single integral knows about the full set of integrals. Indeed, its derivatives generate (in general) all other function in the integral basis. Our task is then to properly organize the information contained in the derivatives into UT functions as well.
In this paper, we show how to do this algorithmically. As a result, starting from a single UT integral (and completing the basis in an arbitrary way), we derive algorithmically a transformation to a canonical form of differential equations (if such a transformation exists). As a corollary, this provides a test of the UT property of an integral. This can be used to find the canonical form for differential equations for individual integral sectors, which is advantageous since this limits the size of then matrices that are needed. At the same time, we wish to emphasize that our method can also be used with just a single integral in the top sector as input, to transform the equations to canonical form in a single step. The applications we present in this paper suggest that this is feasible in practice even for large systems of differential equations.
Our work considerably simplifies the application of the canonical differential equation method, as it is in general much easier to find one UT integral, as compared to a full UT basis. Several methods and techniques exist for finding the integral to start from. Moreover, conjecturally, loop corrections in N = 4 super Yang-Mills are given by uniform weight functions, and therefore this theory can be used as further inspiration for finding a good 'seed' for generating the UT integrals needed for any quantum field theory. Feynman integrals are related to each other by integration-by-parts relations. They can be reduced to a finite number of master integrals, which satisfy a system of homogeneous linear differential equations [19][20][21][22]. Alternatively, the coupled system of master integrals can be described by an n-th order differential equation, the so-called Picard-Fuchs equation, which defines the linear relation between a certain master integral and its higher-order derivatives. In [7], it was proposed that a canonical basis of master integrals exists, which consists of integrals of uniform transcendental weight (UT). They satisfy a linear system of differential equations in a simple, canonical form: where one needs to find all the linearly-independent UT integrals to bring the first-order differential equations into canonical form. On the other hand, one can derive the Picard-Fuchs equation for any given master integral. The coefficients of the Picard-Fuchs equation have some special characteristics, which allow to test if an integral is UT without referring to the other master integrals in the bases. More importantly, by establishing the relation between the Picard-Fuchs equation and the first-order differential equations, one can map between two bases that share one common integral. In this way we can construct a canonical basis g starting from an arbitrary basis f that contains one UT integral.
Our starting point is one UT integral, which without loss of generality we denote by f 1 . We complete this to a basis of n master integrals, denoted by f . Given the system of differential equations we can reformulate this in terms of the function f 1 and its derivatives only. To this end, we take derivatives of eq. (2.2), As we are interested in the higher-order derivatives of f 1 , we project this equation with the vector v 0 = (1, 0, . . . , 0) and define Comparing to eq. (2.3), we see that the matrix Ψ satisfies Assuming that f 1 , . . . , f (n) 1 are linearly independent, 1 we can then invert the matrix Ψ, and write down an n-th order differential equation for f 1 , JHEP05(2020)025 The coefficients b 1 , . . . , b n are rational functions of x, . By factoring out the denominators, the differential equation ( has weight −k. After assigning transcendental weight −1 to , we group together terms in (2.10) that have equal transcendental weight and thus obtain p max + 1 independent linear equations in the unknowns {f 1 , f (m,k) 1 }. If non-trivial solutions exist, then p max + 1 must be less or equal to the total number of unknowns, which is 1 + n m=1 m. In particular, the equation of weight zero only involves f 1 , whose coefficient must vanish, and thereforeb 0 = O( ). To conclude, we find the following necessary conditions for f 1 being a UT integral: These conditions allow us to test whether an integral is UT or not. Similar conditions were also discussed in [18]. Now let us propose a way of solving the differential system by the method of undetermined exponential function. First, we assume the existence of a canonical basis g, where g 1 = f 1 and all the other members of g are unknown. g satisfies the canonical differential equations (2.1) with unknown matrix B(x), whose general solution is defined through an iterated integral Here B(x) dx is the argument of the path-ordered exponential, which needs to be determined. As we shall see in the following section, the linear relations among g 1 (= f 1 ) and its derivatives allows us to determine the matrix B(x) up to a constant similarity transformation. To start, we define a matrix Φ in complete analogy with the definition of Ψ in eq. (2.6), but with A replaced by B ≡ B(x), The latter satisfies (g 1 , g 1 , . . . , g 1 ) T = Φ(x, ) g. (2.14)

JHEP05(2020)025
Comparing eqs. (2.7) and (2.14), we see that the matrices Φ and Ψ define the transformation matrix T between the two bases f and g, Knowing that f 1 = g 1 , the first line of T or T −1 must equal the unit vector v 0 . Explicitly, Recall that Ψ is explicitly known, while Φ depends on the unknown matrix B. The constraint (2.16) will be important for finding B, and hence the transformation to the UT basis T .

Ansatz for canonical differential equations, and determination of unknowns
By assumption, the basis g is UT, and hence the differential equations are in canonical form (2.1). We write the matrix B appearing in that equation as where the m l are constant matrices. The a l are functions behaving as da l ∼ dy/y near singular points y(x) = 0. Their explicit form depends on the Feynman integrals under consideration. In the case of multiple polylogarithms 2 they take the form This form makes it manifest that Feynman integrals have fuchsian singularities only. The set of (rational or algebraic) functions α l depends on the problem under consideration. For a given system of master integrals, one can read off the set of singularities from the differential equations (2.2). In the case of a rational alphabet {α l } this is all we need. For instance, in most of the applications considered below, we have The precise form or number of singularities is not important for our method. For simplicity of notation and presentation, let us for the following assume the rational form (2.19) with three singularities at finite distance. (These assumptions can be dropped in special cases where we need different forms of ansatz for the alphabet, see discussions in section 2.4).
With this information, we see that eq. (2.16) effectively becomes an equation for the constant ( -and x-independent) matrices m l , and products thereof (projected by the vector v 0 ). We may profit from the fact that the matrices are constant by sampling the equations for different values of x. Moreover, the method can naturally be combined with finite field methods, see e.g. [23,24].

JHEP05(2020)025
In summary, we have the constraint (2.16), together with the ansatz (2.17) or (2.18) for the canonical differential equation matrix B. The latter is parametrized via the set of constant matrices m l , and determines Φ, see eq. (2.13). More precisely, each row of Φ defines an n-dimensional vector, which is a polynomial in , Substituting into (2.16), or equivalently (2.10), we obtain a linear equation with xdependent coefficients at each fixed order in . At O( k ), the equation reads The recursive structure of (2.22) allows us to solve (2.21) iteratively and determine the constant matrices m l . Below we give details of this algorithm. Starting from k = 1 and assuming 3 letters in the ansatz, (2.21) and (2.22) lead to Assuming the α l (x) are rational, (2.23) can be sampled for a set of different generic numbers which gives rise to a system of algebraic equations linear in v 0 m l . Solving for the three unknowns, we find where r 1 is the rank of the linear system (2.23) sampled over a set of numbers and β l j are rational numbers.
. . , N 1 . Then we can repeat the analysis for (2.21) at k = 2, treating v i m l as the unknowns.
To generalize, the algorithm introduces a total number of S k−1 ≡ k−1 j=1 N j free vectors through step 1 to k − 1, (k ≥ 2). Moving on to step k, we first substitute the solutions obtained at the previous step into (2.22) and reduce Φ k m onto a set of vectors Evaluating (2.21) over a set of generic numbers then leads JHEP05(2020)025 to a linear system of equations in Q k of rank r k . Solving for the unknowns, we obtain Again, β l ij are rational constants. Since the algorithm constructs linear equations step-by-step, the relations in (2.25) must be independent from those obtained at previous steps. Iterating the above procedure until we reach a certain step where N k = 0, and S k = S k−1 ≡ |S| − 1, the algorithm terminates. Each row of the Φ−matrix is now completely determined as a linear combination of vectors in Q ≡ { v 0 , v 1 , . . . , v |S|−1 }. Solutions in all previous steps combine into 3|S| independent linear relations Q m l = β l Q , where β l (l = 1, 2, 3) are three |S| × |S| matrices, whose components β l ij are rational numbers. In fact, if non-trivial solutions exist, then we must have |S| = n, the rank of the coupled linear system. If |S| is greater than n, some of the vectors in Q are not linearly independent, which contradicts the claim that the 3|S| linear relations we obtained are independent. If |S| is smaller than n, then the rank of Φ-matrix is smaller than n, which contradicts our assumption that f 1 , . . . , f (n) 1 are linearly independent. Therefore |S| = n, Q is an n × n invertible constant matrix and m l = Q −1 β l Q.
Finally, we need to make some explicit choice for the vectors in Q. For example, we could choose v i to be the (i + 1)-th row of the identity matrix, such that the m l = β l . A different choice of v i generates a constant linear transformation P acting on Q, which preserves its first component. Correspondingly, the matrices m l differ by a constant similarity transformation To summarize, the constraint (2.16) suggests a system of linear equations (2.21). We develop an algorithm to solve them iteratively order-by-order in . Given the assumption that f 1 is UT, the Ψ-matrix has rank n and the ansatz for letters α l is complete, the algorithm will terminate at a certain order in , when it finds a non-trivial solution for the m l that determines Φ(x, ) and hence T (x, ), up to a constant linear transformation. Equations at higher orders in must be trivially satisfied. We provide a Mathematica implementation of this algorithm, see section 4.

Generalization to multi-variable case
Unlike methods based on the Moser algorithm, see e.g. [4,11,[14][15][16], the inclusion of multiple scales does not pose a significant problem within our approach. The main point is that the canonical form of the differential equations for multiple polylogarithms still only depends on a finite number of constant matrices m l that are determined by our procedure. Concretely, the multi-variable generalization of eq. (2.18) is (2.27)

JHEP05(2020)025
Here x = {x 1 , . . . , x m } denotes a set of variables, d = i dx i ∂ x i , and the set {α l (x)} is the alphabet. Although the whole differential depends on multiple variables x, one may always view it as a single-variable problem by treating all but one preferred variable (e.g. x 1 ) as constants. The knowledge about the x 1 -dependent letters provides sufficient information to reconstruct the answer. Let us illustrate here how the method works in the two-variable case. We will present a state-of-the-art multi-variable example in section 3.4. Starting with a system of differential equations for a basis f containing one UT integral we make an ansatz for the alphabet {α l (x, y)}, and hence the differential equation for a UT basis g in canonical form: We now study the Pichard-Fuchs equation in x, treating y as a constant y = y 0 . For convenience, we assume that the first L letters depend on x (and possibly y) and the others depend only on y. Using the short-hand notation B x ≡ ∂ ∂x B, the partial derivative B x only depends on m 1 , . . . , m L . Similar to the one-variable case (see (2.6)), by taking partial derivatives we define the matrix Ψ x through A x and likewise Φ x through B x . They satisfy the following relations, The next step is to solve this system of constraints for m l (l = 1, . . . , L). As described in the previous section, solving (2.31) order-by-order in . We will obtain a set of linear relations (2.32), where β l (l = 1, . . . , L) are n × n rational matrices and Q denotes { v 0 , . . . , v n−1 }, For any given Q, this defines our solution for m l . Let us emphasize that the solution is constant and, in particular, independent of y 0 . This is because, by construction, the coefficients β l obtained by the algorithm are unique (the algorithm runs with a specific If there exists a y-independent solution for m l , then it must satisfy (2.32) at y = y 0 , and we have found it. These constant matrices m l , together with the knowledge about the first L letters in {α l }, allow us to determine the x-and y-dependent matrix B x and hence Φ x . In this way we find the transformation matrix T (x, y, ) ≡ Φ −1 x Ψ x , which brings the partial differential matrix A x into canonical form.

JHEP05(2020)025
Now we will argue that g = T f is a UT basis. Let us assume there exists a UT basis g UT , whose first component is f 1 and which satisfies the differential equation ∂ dx g UT = B x g UT . From the analysis in the previous section (see (2.26)), solutions the m l are related by a similarity transformation which leaves its first row invariant. Therefore, there exists a constant transformation matrix P between B x and our solution for B x , such that (2.33) Following from definition (2.13), the matrix Φ x defined through B x is related to Φ x through a constant linear transformation The same constant matrix transforms the UT basis g UT into g Therefore g itself is a UT basis. The above argument holds as long as the x-derivatives for f 1 couple to all master integrals in the family, such that Φ x and Ψ x are invertible. In the process of searching for UT integrals, the algorithm refers only to the partial derivatives in x, thus becoming very efficient. The algorithm does not require a prior knowledge of the complete set of letters. In practice, the x-independent letters can be determined afterwards by transforming (2.28) into the whole differential in canonical form.

Special cases with degenerate Ψ-matrix and algebraic letters
In previous sections we explained how the algorithm applies to standard one-and twovariable examples. For simplicity of presentation, we assumed 1. One can find a UT integral f 1 whose higher-order derivatives couple to all master integrals in the system.
2. The differential equations contain only rational letters, and hence takes the form of (2.19).
In reality, these assumptions can be dropped when needed. Now we will discuss the subtleties of applying our algorithm when 1) and 2) no longer hold.
Regarding assumption 1), we would like to comment on the situation where derivatives of the first UT integral only couple to a subset of master integrals. For an n × n coupled system, depending on our choice of f 1 , the number of its independent higher-order derivatives could be less then n. The corresponding Ψ-matrix is then degenerate. Typically this could happen when f 1 belongs to a sub-topology of the integral family. Sometimes it occurs even if f 1 is in the top sector. One example would be the scalar integral in three-loop ladder integral family, whose derivatives only couple to 23 out of a total number of 26 JHEP05(2020)025 master integrals, see section 3.2. In these situations the algorithm still works, but we need to search for a second UT integral f 2 , such that the union of independent derivatives of f 1 and f 2 couples to all master integrals. One way to proceed can then be to first bring a linearly independent sub-block of the differential equation into canonical form, and then use the second UT integral to work on the remaining part.
Another very efficient approach can be to use both UT integrals simultaneously: for example, one can find a set of linearly independent derivatives {f 1 , . . . , f (n 1 ) 1 , f 2 , . . . , f (n−n 1 ) 2 }, which contains n 1 and n − n 1 higher order derivatives of f 1 and f 2 , respectively. They form a basis of the coupled system. Starting from a master integral basis f ≡ (f 1 , f 2 , . . . , f n ) T , we then construct a matrix Ψ in the same way as in (2.6): where v 0 , v 1 are the first and second row of the identity matrix. The Ψ-matrix thus defined is invertible and satisfies the relation Given the two UT integrals that are already known, we can search for a UT basis g = (g 1 , g 2 , . . . , g n ) T , where g 1 = f 1 , g 2 = f 2 . As we discussed before, g satisfies (2.1), and therefore (g 1 , . . . , g (n 1 ) 1 , g 2 , . . . , g (n−n 1 ) 2 where Φ is defined in the same way as Ψ in (2.36) with A replaced by B. (2.37) and (2.38) imply the following system of constraints, By solving these constraints, we find the transformation between the two bases f and g.
Next, we will come to assumption 2) about rational letters. In the cases where 2) holds, we should allow the ansatz for the B-matrix to take a more general form compared with the oversimplified version in (2.19). In particular, the differential equations might contain fuchsian singularities at zeros of a certain higher-degree polynomial, which can be algebraic and complex numbers (e.g. sixth-root of unity). More generally, we can assume a factorized form with k singularities,

JHEP05(2020)025
where the x i are roots of a degree-k polynomial P k (x) with real and rational coefficients. For the purpose of analyzing equations (2.21) by the finite-field method, we need to avoid writing down an ansatz for the Φ-matrix that explicitly contains algebraic numbers. Thus, it is advantageous to make the ansatz for B in the following form, Now that (2.41) guarantees that only rational functions and numbers appear in (2.21), we expect to find solutions for the constant matrices m l whose elements are real and rational numbers. We can again use finite fields to search for such solutions, as explained in section 2.2. So far we have limited our discussion to differential equations with rational letters, where one can always find a rational transformation matrix T (x, ) that brings (2.2) to canonical form. There are also situations where a rational transformation does not exist and hence square roots are required. For rational alphabets, the individual alphabet letters can be read off from the original differential equation. The case of algebraic dependence on the kinematics is more subtle, as the latter are not immediately apparent. However, studying the singular behavior of the differential equation (in particular, the critical exponents in each limit) provides this additional information. For example, if the differential equation (2.2) is in fuchsian form, we can expand it around any singular point The exponential of the coefficient matrix B i determines the asymptotic behaviour of the solution to (2.2) as x → x i , At = 0, if B i contains non-integer eigenvalues, they cannot be transformed away via the so-called balance transformation [11]. Put differently, the square roots in (x − x i ) B i cannot be removed by a rational linear transformation on the basis, and therefore √ x − x i must be included in our ansatz for the alphabet. Repeating the analysis for all fuchsian singular points should allow us to find all the square-root letters in the alphabet. In many cases, it is possible to find a set of variables that rationalizes all letters, (see e.g. [25]). After such a change of variables the problems simplify into those with only rational letters. In this way we apply the algorithm to work out a four-variable example with square-root letters, as we will demonstrate in 3.4.
Let us emphasize that the method of our algorithm should work for algebraic letters without rationalization. In this case the Pichard-Fuchs equation will still be rational, but the Φ-matrix, and hence the coefficients of system of linear equations for the m l , will be algebraic. Nevertheless, the constant m l matrices will not depend on the square roots. One may still apply (a variation of) our algorithm to search for rational solutions and provide a quick test of the conjectured algebraic alphabet letters. JHEP05(2020)025

Examples and applications
In this section we present examples and applications of our algorithm. We use relatively standard notation for Feynman integrals in the context of differential equations. The integral families are defined as where L is the number of loops, γ E is the Euler-Mascheroni constant, the denominator factors D i are defined in the respective subsections and we are using the mostly minus metric (+ − · · · −). In subsection 3.1 we apply the algorithm to two three-loop four-point integral families that first have been computed in ref. [26]. In subsection 3.2 we present a new result for a four-loop four-point integral and in subsection 3.3 we bring a non-planar four-loop sector that appears in the computation of the angle-dependent cusp-anomalous dimension into canonical form. Finally, as a multi-variable example, we apply the algorithm to the top sector of a non-planar two-loop five-point family which was computed in [17,27]. A summary of the performance of our implementation on a standard desktop computer with twelve logical cores can be found in table 1.
The IBP reductions necessary for computing the initial differential equations were done either by the integral reduction program FIRE6 [28] or by using the tools available in the FiniteFlow framework [24]. Both codes depend on LiteRed [29] to generate the IBP identities.

Full differential equation for planar three-loop integrals
As a first example, we apply the algorithm to the two three-loop four-point integral families shown in figure 1. The definition of the factors in eq. (3.1) is for integral families 1(a), and 1(b), respectively. The top sector is defined by 1,1,1,1,1,1,1,1,1,0,0,0,0,0 in both cases. The integrals were computed previously in ref. [26] with the differential equations method. In this case it is relatively straightforward to find a complete UT basis as in [26], or even a complete dlog integrand basis [30]. Nevertheless, we find it instructive to benchmark our new method using these sets of integrals. We will see that, for each integral family, a single UT integral from the top sector is sufficient to derive the full canonical differential equation. The corresponding matrices are of size 26×26 and 41×41, respectively. A suitable initial integral is easily found using [30], or by taking inspiration from the perturbative expansion of N = 4 super Yang-Mills [31].
Concretely, we took the following integrals as our starting point, for integral families 1(a), and 1(b), respectively. We completed them to basis by taking linearly independent integrals suggested by the integral reduction programs. In other words, no optimization was done on the other integrals. The integrals depend on one dimensionless variable x = t/s. See [26] for more details. The differential equation matrix A(x) has the singular points x = 0, −1, ∞, and consequently we take the alphabet in eq. (2.18) to be α = {x, 1 + x}. With this as input, our algorithm effortlessly found the transformation matrix T , see table 1.

New result for a four-loop four-point integral
Let us now present an application to previously unknown four-loop integrals. The definition of the integral family is JHEP05(2020)025 The sector shown in figure 2 is G 1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0 and, together with the subsectors, there are 19 MI. This case is particularly interesting for the following reason. After taking certain residues, the scalar integrand exhibits a double pole. This can be understood from a power counting argument and comes from the fact that the integral has relatively few propagators. As a consequence, the scalar integral, or integrals with the same propagator structure, and numerators, do not have a dlog form in four dimensions. There are possible remedies to this, including writing down a dlog basis for the 'top' sector shown in figure 2, and then solving for the full family of integrals.
Instead, here we wish to use the new tool for testing UT integrals to find directly UT integrals in the sector of figure 2. We proceed as follows. First, we look directly at the scalar integral G 1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0 . Although this turns out to be a suitable integral (it is UT up to an overall normalization in ), only 18 of its derivatives are linearly independent. One way to proceed can be to first bring an 18 × 18 block of the matrix into canonical form as described in section 2.4. After this, one can apply the ideas of [4,11] to the last row and obtain a canonical form for the full differential equation.
Instead, we choose to make use of the ability of our algorithm to test for further suitable candidates. Heuristic rules for finding UT integrals from ref. [26] suggest that doubling a propagator is promising for this type of integral with off-shell triangle subintegrals. Looking at figure 2, we see that there are two inequivalent ways of doubling one propagator. Checking both with our algorithm, we find that each of them is a UT integral (up to an overall normalization in ) but that only one of them has 19 linearly independent derivatives: 1,1,1,1,1,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0 (3.7) Note that the factor of (1 + x) can easily be found through testing possible factors or by integrating out the 0 part ofb 0 , as described in section 2.2. Using this integral as our starting point, our algorithm takes less than a minute to find the canonical basis.  Solving the DE (2.1) in terms of iterated integrals as in (2.12) is then straightforward. The alphabet, α = {x, 1+x}, suggests that the solution can be written in terms of harmonic polylogarithms (HPLs) [32] with indices {0, −1} only.

Canonical form for non-planar four-loop sector with 17 master integrals
Here we discuss one of the most complicated applications of our algorithm. In the previous cases, the maximal size of individual sectors (i.e., the number of coupled integrals) was three in section 3.1, and twelve in section 3.2. In contrast, here we will apply the algorithm to a case of a sector with 17 coupled master integrals. It is shown in figure 3. The definition JHEP05(2020)025 of the integral family is where we consider the sector G 1,0,1,1,0,1,1,1,1,0,0,1,0,0,0,0,0,0 . Inspecting the differential equation for the cut integral, we identify the alphabet of the sector to be {x, The algorithm needs one UT integral to start with. We used the following procedure to find this integral: 1. We use heuristic rules to find likely UT candidates (see [26,36] for more information).
2. We use our algorithm to test the UT property for each integral individually. If only the appropriate normalization factor is missing, we find it in the same ways as mentioned in the previous subsection.

Four-variable example: non-planar double pentagon integrals
Here we apply the algorithm to the cutting-edge example of a coupled system of two-loop five-point integrals, whose differential equation depends on four dimensionless variables [17,27]. Consider the non-planar double-pentagon integral family figure 4, where the inverse propagators are,

JHEP05(2020)025
We focus on the top sector G 1,1,1,1,1,1,1,1,0,0,0 which contains 9 master integrals. We start with an ansatz for the alphabet containing 31 letters {W i }, suggested by [37,38]. The five-point external kinematics can be parametrized via (a variation of) momentum-twistor variables b = {b 1 , . . . , b 5 }, which rationalize all letters of the alphabet, see e.g. [39]. 12) where b 1 sets the overall kinematic scale. The differential equation depends on four dimensionless variables b 2 , . . . , b 5 . The algorithm takes derivatives with respect to one preferred variable. It is convenient for us to choose this variable to be b 2 . Given the ansatz {W i }, the canonical partial differential matrix B 2 (b, ) ≡ ∂ ∂b 2 B depends on 22 independent letters, (3.14) To proceed, we need to select a suitable UT integral as an input to the algorithm. In order to investigate the parity dependence of the differential equations, we tested the algorithm starting with both a parity-even and parity-odd integral taken from the canonical basis given in [17]. Fixing the value of b 3 , b 4 , b 5 to be certain constants, we execute the algorithm to find a solution for m 1 , . . . , m 22 , and hence for B 2 (b, ). In order to determine the transformation matrix T (b, ), as defined in (2.15), normally one would compute the Ψ and Φ matrix with full analytic depedence on all variables. This is the most time-consuming step of the algorithm since taking higher derivatives generates rational coefficients that depend on high-degree polynomials. In this case, we find it more convenient to first compute the matrices B i ≡ ∂ ∂b i B, i = 1, . . . , 5, which allow us to determine the whole differential in (2.27) analytically.
We proceed in the following way: • Set b 4 and b 5 to constants. Compute Ψ 2 , Φ 2 and hence the T matrix with analytic dependence on b 3 . Transform the original partial differential equation in b 3 into canonical form. Thus we obtain B 3 where all other variables are set to constants.
• By matching onto the ansatz for the relevant alphabet letters, we extract the corresponding matrix-residues in B 3 . In this way we reconstruct B 3 with analytic dependence on all variables.
• Construct an intermediate integral basis h to relate the canonical basis g to the initial basis f . For example, we set h = ( f 1 ) which contains only first-and second-order partial derivatives.
• Find the linear transformation between h and g through the B i matrices. Likewise, work out the linear relations between h and the initial basis f through the original differential equations. The way we construct h guarantees that these relations can be easily obtained analytically.
• Compute the transformation matrix T (b, ) between g and f through their relations to h.
Starting from either a parity-even or parity-odd UT integral, the algorithm finds the solution to transform the differential equations on the maximal cut into canonical form, which depend on 17 letters W i , i ∈ {1, . . . , 5, 11, 16, . . . , 20, 26, . . . , 31}. The total running time for solving the double-pentagon example is of the order of minutes. As a starting point, the algorithm needs to know one UT integral. This information can be obtained from D-dimensional Baikov representation [17]. Alternatively, one could start with a dlog integral in six dimensions [27]. Compared with the methods in the literature, our algorithm requiries minimum input from the integrand analysis [30]. This feature makes the algorithm particularly suitable for dealing with multi-variable problems.
The performance of our algorithm on all examples of this section is summarized in table 1.

Public implementation
We provide a Mathematica package INITIAL (an INitial InTegral ALgorithm) which utilizes FiniteFlow [24] to perform operations of our algorithm over finite fields. The package is publicly available at https://github.com/UT-team/INITIAL It relies on the FiniteFlow library [24] and its dependencies. The examples mentioned in the previous section can also be downloaded from the same repository.

Conclusion and outlook
The automated calculation of Feynman integrals in quantum field theory is of considerable interest to the scientific community. In this paper, we developed further an idea due to [18], which in turn relies on the method of canonical differential equations [7].

JHEP05(2020)025
In other approaches, one needs a full set of uniform weight integrals to obtain the canonical differential equations. In the new approach, only one such integral is needed: the remaining ones are obtained from the former algorithmically, if such a basis exists. As a corollary, the new approach provides an algorithm to test whether the candidate integral has uniform weight. A necessary condition is given in (2.11). Moreover, the existence of a canonical transformation verifies the uniform weight property of the candidate integral.
We expressed the necessary equations in matrix form. The equations are easily handled, and we explained how to solve them systematically. We found that this implementation can be readily used for state-of-the-art problems. We used it to find the canonical form of complicated systems of differential equations, e.g. with 17 coupled master integrals in one integral sector (on the cut). We explained how the algorithm deals with multivariable differential equations in an efficient way. For demonstration, we applied it to a cutting-edge two-loop five-point example. Given one UT integral, the algorithm finds the canonical transformation for the non-planar double-pentagon integrals on the maximal cut, which depends on four dimensionless variables.
Our work opens up several interesting directions for further developments.
Canonical forms for elliptic polylogarithms: the idea of the canonical form of differential equations has also been explicitly applied for elliptic polylogarithms. Specifically, there are two approaches that seem promising in this regard. Firstly, it has been argued [4] (see also [40]) that a pre-canonical form should exist for all Feynman integrals, of the following type, d g(x) = [dA 0 (x) + dA 1 (x)] g(x) , where the matrices A 0 and A 1 only involve logarithms (i.e., the fuchsian property of the differential equations is manifest for all singular points). In the polylogarithmic case, 'integrating out' the A 0 part can be done using algebraic functions only, but in the elliptic case (and beyond) this leads to special functions. In the latter case, it has been shown explicitly that allowing non-algebraic transformation matrices one finds the canonical form of [7] also in the elliptic case, where A(x) contains functions beyond logarithms [41,42].
Application to finite integrals: in most applications, for example when computing finite cross sections, or scaling dimensions of operators, one is ultimately interested in a finite, four-dimensional quantity. In some situations, it is possible to explicitly separate divergent parts of the calculation and express the remainder in terms of manifestly finite integrals. It is therefore interesting to apply the method for finite integrals [43]. The latter occur e.g. when making the infrared properties of scattering amplitudes manifest. Indeed, it is very natural within the dlog integrand approach to classify integrals according to their infrared properties, with the infrared finite integrals typically being the most interesting JHEP05(2020)025 ones [7,8]. Moreover, finite integrals occur frequently when dealing with problems with several mass scales. The main simplification when dealing with finite integrals is that the weight expansion truncates: there is a maximal weight occurring in the calculation. The way this happens in practice is that the differential equation matrix becomes nilpotent. Moreover, the number of master integrals reduces in this case [44] (some integrals decouple), and it has also been observed that the required function alphabet may simplify. Due to these simplifications, in ref. [43] a full massive three-loop calculation was possible. We find it promising to combine this method with our new approach.
Uniform weight as guiding principle for recurrence relations: finally, let us mention that Feynman integrals satisfy dimensional recurrence relations [45]. For single-scale integrals, where differential equations can only be applied indirectly [46], they are one promising method to evaluate Feynman integrals. See e.g. [47,48] for recent applications. We find it likely that the uniform weight information of a single UT integral provides important input for that method as well, and provides a useful organizing principle for the calculation.
We wish to close with a discussion on the role played by N = 4 super Yang-Mills in our understanding of Feynman integrals. In many calculations of physical quantities in that theory it was observed that the answer typically is given by UT functions. Whether or not this is true in general is an open question. One may wonder whether the integrals encountered in N = 4 super Yang-Mills are simpler or more complicated, with respect to QCD. Of course, when talking about full amplitudes or finite quantities in N = 4 super Yang-Mills, the latter often have additional symmetries or hidden properties, and hence are definitely special. However, what about individual Feynman integrals in dimensional regularization? The answer is that generic QCD integrals are just as nice, or simple, as the ones in N = 4 super Yang-Mills, at least if one organizes them in a suitable way, as we have learned in the course of the last ten years. With the present work, it becomes clear that a stronger statement is possible: not only are generic integrals as simple as the ones in N = 4 super Yang-Mills, but in fact they can all be obtained from the knowledge of the former!