Magnus and Dyson Series for Master Integrals

We elaborate on the method of differential equations for evaluating Feynman integrals. We focus on systems of equations for master integrals having a linear dependence on the dimensional parameter. For these systems we identify the criteria to bring them in a canonical form, recently identified by Henn, where the dependence of the dimensional parameter is disentangled from the kinematics. The determination of the transformation and the computation of the solution are obtained by using Magnus and Dyson series expansion. We apply the method to planar and non-planar two-loop QED vertex diagrams for massive fermions, and to non-planar two-loop integrals contributing to 2 ->2 scattering of massless particles. The extension to systems which are polynomial in the dimensional parameter is discussed as well.


Introduction
The method of differential equations (DE's), developed by Kotikov, Remiddi and Gehrmann [1][2][3] and reviewed in Ref. [4,5], is one of the most effective techniques for computing dimensionally regulated multi-loop integrals and has led to significant achievements in the context of multi-loop corrections. Within the continuous dimensional regularization scheme, Feynman integrals can be related by using integration-by-parts identities (IBP-id's) [6,7], Lorentz invariance identities [3], Gram identities [8], and quasi-Shouten identities [9]. These relations can be exploited in order to identify a set of independent integrals, dubbed master integrals (MI's), that can be used as a basis of functions for the virtual contributions to scattering amplitudes.
The MI's are functions of the kinematic invariants constructed from the external momenta, of the masses of the external particles and of the particles running in the loops, as well as of the number of spacetime dimensions. Remarkably, the existence of the aforementioned relations forces the MI's to obey linear systems of first-order differential equations in the kinematic invariants, which can be used for the determination of their expression. When possible, these systems are solved exactly for generic values of the space-time dimension D. Alternatively, they can be Laurent-expanded around suitable values of the dimensional parameter up to the required order, obtaining a system of chained differential equations for the coefficients of the expansions. In the most general case, the latter are finally integrated by using the method of Euler's variation of constants.
The nested structure of the Laurent expansion of the linear system leads to an iterative structure for the solution that, order-by-order in ǫ = (4 − D)/2, is written in terms of repeated integrals, starting from the kernels dictated by the homogeneous solution. The transcendentality of the solution is associated to the number of repeated integrations and increases by one unit as the order of the ǫ-expansion increases. The solution of the system, namely the MI's, is finally determined by imposing the boundary conditions at special values of the kinematic variables, properly chosen either in correspondence of configurations that reduce the MI's to simpler integrals or in correspondence of pseudo-thresholds. In this latter case, the boundary conditions are obtained by imposing the regularity of the MI's around unphysical singularities, ruling out divergent behavior of the general solution of the systems.
For any given scattering process the set of MI's is not unique, and, in practice, their choice is rather arbitrary. Usually MI's are identified after applying the Laporta reduction algorithm [10]. Afterward, convenient manipulations of the basis of MI's may be performed.
Proper choices of MI's can simplify the form of the systems of differential equations and, hence, of their solution, although general criteria for determining such optimal sets are not available. An important step in this direction has been recently taken in Ref. [11], where Henn proposes to solve the systems of DE's for MI's with algebraic methods. The key observation is that a good choice of MI's allows one to cast the system of DE's in a canonical form, where the dependence on ǫ, is factorized from the kinematic. The integration of a system in canonical form trivializes and the analytic properties of its general solution are manifestly inherited from the matrix associated to the system, which is the kernel of the representation of the solutions in terms of repeated integrations.
This novel idea has been applied in a number of cases by Henn, Smirnov, and Smirnov [13][14][15], showing the effectiveness of this approach. As pointed out in [11], finding an algorithmic procedure which, starting from a generic set of MI's, leads to a set MI's fulfilling a canonical system of DE's is a formidable task. In practice, the quest for the suitable basis of MI's is determined by qualitative properties required for the solution, such as finiteness in the ǫ → 0 limit, and homogeneous transcendentality, which turn into quantitative tools like the unit leading singularity criterion and the dlog representation in terms of Feynman parameters [12].
In this article, we suggest a convenient form for the initial system of MI's, and we propose an algorithm to find the transformation matrix yielding a canonical system. In particular, we choose a set of MI's obeying to a system of DE's which has a linear ǫ-Let us assume that H can be split in two terms as where H 0 is a solvable Hamiltonian and ǫ ≪ 1 is a small perturbation parameter. We may move to the interaction picture by performing a transformation via a unitary operator B.
In this representation any operator A transforms according to In the interaction picture one imposes that only H 1 (H 0 ) enters the time evolution of the states (of the operators), thus B is obtained by imposing In the interaction picture the Schrödinger equation can be cast in a canonical form, where the ǫ-dependence is factorized. If the Hamiltonian H 0 at different times commute, the solution of Eq. (2.5) is The important remark in this derivation is that, as a consequence of the linear ǫdependence of the original Hamiltonian Eq. (2.2), the states fulfill an equation in a canonical form by means of a transformation matrix B that obeys the differential equation (2.5). This simple quantum mechanical example contains the two main guiding principles for building canonical systems of differential equations for Feynman integrals: • choose a set of Master Integrals obeying a system of differential equations linear in ǫ; • find the transformation matrix by solving a differential equation governed by the constant term.
In this example H 0 (t) and B(t) commute. In the case of Feynman integrals, no assumption can be made on the properties of the matrix associated to the systems of DE's built out of IBP-id's. Therefore, in the following, we need to consider the generic case of noncommutative operators.

Magnus series expansion
Consider a generic linear matrix differential equation [18] in the scalar case, the solution can be written as In the general non-commutative case, one can use the Magnus theorem [16] to write the solution as, where Ω(x) is written as a series expansion, called Magnus expansion, The proof of the Magnus theorem is presented in the Appendix A, together with the actual expression of the terms Ω n . The first three terms of the expansion (3.4) read as follows: , , We remark that if A and its integral commute, the series (3.4) is truncated at the first order, Ω = Ω 1 , and we recover the solution (3.2). As a notational aside, in the following we will use the symbol Ω[A](x) to denote the Magnus expansion obtained using A as kernel.

Magnus and Dyson series expansion
Magnus series is related to the Dyson series [18], and their connection can be obtained starting from the Dyson expansion of the solution of the system (3.1), (3.6) in terms of the time-ordered integrals Y n . Comparing Eq. (3.3) and (3.6) we have and the following relations The matrices Q (j) n are defined as In the following, we will use both Magnus and Dyson series. The former allows us to easily demonstrate how a system of DE's, whose matrix is linear in ǫ, can be cast in the canonical form. The latter can be more conveniently used for the explicit representation of the solution.

Differential equations for Master Integrals
We consider a linear system of first order differential equations where f is a vector of MI's, while x is a variable depending on kinematic invariants and masses. We suppose that A depends linearly on ǫ, and we change the basis of MI's via the Magnus series obtained by using A 0 as kernel, which, analogously to the quantum-mechanical case, Eq. (2.5), implies that the new basis g of MI's fulfills a system of differential equations in the canonical factorized form, The matrixÂ 1 is related to A 1 by a similarity map, and does not depend on ǫ. The solution of Eq. (4.5) can be found by using the Magnus theorem with ǫÂ 1 as kernel where the vector g 0 corresponds to the boundary values of the MI's. Therefore, the solution of the original system Eq. (4.1) finally reads, Let us remark that the previously described two-step procedure is equivalent to solving, first, the homogeneous system whose solution reads, and, then, to find the solution of the full system by Euler constants' variation. In fact, by promoting g to be function of x, and by requiring f to be solution of Eq. (4.1), one finds that g(ǫ, x) obeys the differential equation (4.5).
The matrix B 0 , implementing the transformation from the linear to the canonical form, is simply given as the product of two matrix exponentials. Indeed one can split A 0 into a diagonal term, D 0 , and a matrix with vanishing diagonal entries N 0 , The transformation B is then obtained by the composition of two transformations whereN 0 is given byN In the last step of Eq. (4.13) we have used the commutativity of the diagonal matrix D 0 with its own integral. The leftmost expansion performs a transformation that "rotates" away D 0 , while the second expansion gets rid of the O(ǫ 0 ) contribution coming fromN 0 , i.e. coming from the image of N 0 under the first transformation.
In the examples hereby discussed it is possible, by trials and errors, to find a set of MI's obeying a system of DE's linear in ǫ. Moreover in these cases one finds that Ω[N 0 ] contains just the first term of the series, except for the non-planar box, where also the second order is non vanishing.

One-Loop Bhabha scattering
The calculation of the one-loop Bhabha scattering within the DE's method was discussed in [28,29], and more recently in Ref. [14]. A selection of the Feynman diagrams contributing to this process is depicted in Fig. 1. In this section, we compute a set of MI's with a slightly different definition from the ones in [14], which will be also adopted for the the one-loop × one-loop subtopologies of the QED vertices in the next section. The diagrams depend on the invariants s = (p 1 +p 2 ) 2 , t = (p 1 +p 3 ) 2 , u = (p 2 +p 3 ) 2 and on the fermion mass m. Momentum conservation and the on-shellness of the external legs render these variables not independent as they are related by the condition s + t + u = 4m 2 . The integrals can be expressed in terms of the Landau auxiliary variables x and y, defined as follows We identify the following basis f of scalar integrals, in terms of the integrals T in Fig. 2. The basis f fulfills the following systems of differential equations (σ = x, y) Both systems are linear in ǫ and in both cases the O(ǫ 0 ) term, D σ,0 , is diagonal. The systems can be brought in the canonical form by performing the transformation The new basis g, fulfills the canonical systems ∂ x g(ǫ, x, y) = ǫÂ x,1 (x, y) g(ǫ, x, y) , ∂ y g(ǫ, x, y) = ǫÂ y,1 (x, y) g(ǫ, x, y) , The two systems of DE's in Eq.(5.6) can be combined in a full differential form, along the lines of Ref. [14], dg(ǫ, x, y) = ǫ dÂ 1 (x, y) g(ǫ, x, y) , where the matrixÂ 1 fulfills the relations, and the integrability condition The matrixÂ 1 is logarithmic in the variables x and y, The position of the non-zero entries of the sparse matrices M i agrees with the result obtained in Ref. [14]. The actual value of the non-zero entries, however, are different, owing to the different normalization of the elements of the basis of MI's. The solution of the system (5.8) can be computed along the lines of Ref. [14]. In particular, the solution is computed in the Euclidean region 0 < x, y < 1 by using the analytic structures of the g i and then extended in the physical region by analytic continuation [29].

Two-Loop QED Vertices
A basis of MI's for the electron form factor at two loops in QED [20] was computed in Ref. [19], for arbitrary kinematics and finite electron mass. The diagrams contributing to such corrections are depicted in Fig. 3 and depend on s = (p 1 + p 2 ) 2 and p 2 the auxiliary variable x, defined through The canonical form can be obtained performing the transformation described in Section 4, The new basis g is given by where The new basis of MI's obeys a system of DE's in the canonical form,    Figure 5: Non-planar two-loop diagram with massless internal propagators, and massless external particles. The internal momenta shown in the diagram are oriented according to the arrows. All the external momenta are incoming.
The solution of the system can be expressed as Dyson series, as well as Magnus series, in terms of one-dimensional Harmonic Polylogarithms (HPL's) [30]. The requirements that the MI's are real-valued in the Euclidean region and regular in x = 1 (s = 0), or simply the matching against the known integrals at x = 1, fix all but three boundary conditions, corresponding to the constant MI's g 1 , g 4 and g 7 (that do not depend on x). The integrals g 1 and g 4 can be easily computed by direct integration, while g 7 can be determined from the results of Ref. [31]. Our results were checked analytically, using the code HPL [32,33], against the results available in the literature [19]. The expressions of the transcendentally homogenous MI's g are shown in Appendix B, and collected in the ancillary file <vertex2L.m>.

Two-Loop non-planar Box
The evaluation of the two-loop non-planar box diagram in Fig. 5, contributing to the 2 → 2 scattering among massless particles, has already been considered in the literature [21,22]. Recently, for its planar partner, a set of MI's with homogeneous transcendentality was presented in Ref. [11]. Our method can be easily applied to it, but instead of showing the case of the ladder-box diagram, in this section, we compute the additional MI's required for determining the non-planar contribution, having expressions with manifest homogeneous transcendentality as well.
The integrals, in this case, are functions of the invariants s = (p 1 + p 2 ) 2 , t = (p 1 + p 3 ) 2 , and u = (p 2 + p 3 ) 2 , with p 2 i = 0, and s + t + u = 0. We adopt the following initial choice of MI's, where the integrals T i correspond to the diagrams in Fig. 6. We notice that the integrals f 1 , . . . , f 9 are common to the two-loop planar box diagram [11]. The set f of MI's obeys a system of differential equations the variable x, defined as, which is linear in ǫ. According to the procedure in Section 4, we can build the matrix B 0 (x) ruling the change of basis f (ǫ, x) = B 0 (x)g(ǫ, x), so that the new MI's, 3) obey the canonical system, with T g (s, t) Figure 6: MI's for the two-loop diagram in Fig. 5. All the external momenta depicted are incoming. In the last integral the loop momenta have to be fixed according to Fig. 5 and a term (k 2 + p 1 ) 2 enters the numerator of its integrand. A dot indicates a squared propagator.
The solution of the system can be expressed as Dyson series, as well as Magnus series, in terms of one-dimensional HPL's [30]. All MI's have been computed in the scattering kinematics, i.e. s > 0, t < 0, u < 0 with |s| > |t|, which gives 0 < x < 1. As long as the planar sub topologies are concerned, one can fix the boundary conditions using the regularity properties of the integrals in some special kinematical points. On the other hand, the analyticity structure of the crossed box is more complicated, since it involves at the same time cuts in all three Mandelstam variables s, t, u. Nevertheless, in this particular case, the boundaries can be fixed by direct comparison with the results presented in [21,22]. The expressions of the transcendentally homogeneous MI's g are shown in Appendix C, and collected in the ancillary file <xbox2L.m>.

Polynomial ǫ dependence
The cases discussed above admitted an initial choice of MI's f obeying a system of differential equations linear in ǫ. We cannot be sure that this feature is general, and holds for any scattering process in dimensional regularization. Nevertheless, the use of Magnus series enables us to generalize our algorithm to the case of systems of DE's whose matrix is a polynomial in ǫ. In fact, let us consider a system of equations where A is of degree κ in ǫ, By iterating the algorithm described in Section 4, the solution of the differential equation (8.1) can be expressed in terms of a chain of products of Magnus exponentials, where the kernelÂ k is defined aŝ It is worth to observe that, within our construction, the solution f is given by repeated transformations. Starting from we iteratively write f k as, The generalization of the canonical system Eq. (4.5) is obtained at the last step of the iteration, when k = κ − 1, It is important to remark that the complete factorization of ǫ is achieved only if κ = 1, i.e. if the system is linear in ǫ, because, althoughÂ 1 is independent of ǫ,Â k acquires a dependence on ǫ for k > 1, cfr. Eq. (8.3).
The algorithm here described has a wide range of applicability and can be used to compute generic sets of MI's, provided that the matrix associated to the system of DE's can be Taylor expanded around ǫ = 0. In this case, the MI's are obtained perturbatively by truncating the ǫ expansion of the matrices associated to the systems of DE's.

Conclusions
In this article we elaborated on the method of differential equations for Feynman integrals within the D-dimensional regularization scheme.
The freedom in the choice of the MI's allowed us to analyze the paradigmatic case of systems of differential equations whose matrix is linear in the dimensional parameter, ǫ = (4−D)/2. We showed that these systems admit a canonical form, where the dependence on ǫ is factorized from the kinematic variables, as recently suggested by Henn.
We used Magnus series to obtain the matrix implementing the transformation from the linear to the canonical form. The solution of the canonical system is obtained by using either Dyson series or Magnus series. Both series require multiple integrations which allow one to naturally express the MI's in terms of polylogarithms and of their generalization.
We demonstrated that the one-loop Bhabha scattering, the two-loop electron form factors in QED and the two-loop 2 → 2 massless scattering exhibit a basis of MI's leading to linear systems of DE's. We then obtained the corresponding canonical bases, in terms of uniform transcendentality functions.
Finally, we have shown that our procedure can be extended to the more general case of systems of DE's that are polynomial in ǫ.
The range of applicability of the algorithm is rather wide and can be used to compute generic sets of MI's, provided that the matrix associated to the system of DE's can be Taylor expanded around ǫ = 0.

A. Magnus Theorem
We closely follow the discussion of Ref. [35]. Given an operator, Ω, we define the derivative of Ω k w.r.t. Ω by its action on a generic operator H: This definition guarantees that, when Ω = Ω(x) and H = ∂ x Ω, The definition (A.
The last equation can be obtained by induction using the relation The exponential of a matrix Ω is defined via a series expansion: The derivative and the inverse of the exponential of a matrix can be straightforwardly obtained by using the previous results: Lemma A.1 (Derivative of the exponential) The derivative of the matrix exponential can be derived from its action on a generic operator H and reads as follows Lemma A.2 (Inverse of the exponential) If the eigenvalues of ad Ω are different from 2ℓπi with ℓ ∈ {±1, ±2, . . .}, then d exp Ω is invertible, and where β k are the Bernoulli numbers, whose generating function is We have now all the ingredients to prove the following [16] Theorem A.1 (Magnus) The solution of a generic linear matrix differential equation can be written as where Ω(x) can be computed by solving the differential equation, Proof Let us consider the derivative of (A.11). Using the definition (A.6) and the property (A.2) we have The relation (A.12) is thus proven by applying the operator d exp −1 Ω to both sides of (A.13).
The differential equation for Ω explicitly reads, [Ω, [Ω, A(x)]] + . . . , (A. 14) and the solution can be written as a series, called Magnus expansion, The coefficients β j are the Bernoulli numbers while the integrands S

B. Master Integrals for the two-loop QED vertices
In this Appendix we collect the 17 MI's of the two-loop QED vertices introduced in Eq. (6.5). In Section 6, we have obtained them starting from the integrals T i depicted in Fig. 4, which are normalized according to the integration measure (Minkowskian metric is understood) The MI's exhibit uniform transcendentality. In the following we present the expression of the coefficients of their expansion around ǫ = 0 up to O(ǫ 4 ). The coefficients g (a) i are defined as follows:

C. Master Integrals for the two-loop non-planar Box
In this Appendix we present the expression of the 12 MI's of the two-loop non-planar Box in Eq. (7.3). They are obtained according to the procedure described in Section 7 starting from the integrals T i in Fig. 6. The latter are normalized according to the integration measure (Minkowskian metric is understood) In the following we collect the coefficients g The MI's have uniform transcendentality, as can be explicitly checked by using the expressions of the coefficients g (a) i . ,