Automated Solution of First Order Factorizable Systems of Differential Equations in One Variable

We present an algorithm which allows to solve analytically linear systems of differential equations which factorize to first order. The solution is given in terms of iterated integrals over an alphabet where its structure is implied by the coefficient matrix of the differential equations. These systems appear in a large variety of higher order calculations in perturbative Quantum Field Theories. We apply this method to calculate the master integrals of the three--loop massive form factors for different currents, as an illustration, and present the results for the vector form factors in detail. Here the solution space emerging is given by the cyclotomic harmonic polylogarithms and their associated special constants. No special basis representation of the master integrals is needed. The algorithm can be applied as well to more general cases factorizing at first order, which are based on more general alphabets, iterated integrals and associated constants.


Introduction
The fundamental objects in any gauge theory are the scattering amplitudes or correlation functions, as they allow to compute the scattering cross sections for collider experiments at large facilities like the Large Hadron Collider (LHC) at CERN. Computations of such objects are mostly using the diagrammatic approach. Especially, in the case of perturbative Quantum Chromodynamics (QCD), one calculates these objects by obtaining all Feynman diagrams at each order in the expansion coefficient, the strong coupling constant α s . Decades of dedicated work have made it possible to partially automate this procedure, from generating Feynman diagrams to the momentum-integral structure. The main remaining step consists in the computation of the integrals over loop momenta and to perform the associated Feynman parameter integrals.
Through the reduction of the whole problem by integration-by-parts (IBP) techniques [1][2][3][4][5][6][7][8] one obtains master integrals (MIs). One method to solve these integrals is the method of differential equations [9][10][11][12]. Differentiating with respect to a parameter in the system one obtains coupled systems of ordinary differential equations of master integrals in the uni-variate case, with which we deal with in the following. 1 In the case where these systems factorize at first order, the complete solution can be constructed algorithmically. This has been done before in Ref. [12] mapping to systems of difference equations, which also factorize to first order. The solution has then been performed using difference ring and field technologies [14][15][16][17][18][19][20][21][22][23][24][25][26], implemented in the package Sigma [27,28].
In the present paper, we present an algorithm operating on uni-variate systems of differential equations, which are factorizing at first order, directly. In the case where the factorization of the system leads to higher order sub-systems, elliptic and even more involved structures will appear, cf. e.g. [29][30][31][32][33][34][35][36][37]. Here still iterative solutions can be found. However, the corresponding integrals contain also letters, which are given by non-iterative integrals and therefore these solutions are given by iterative non-iterative integrals [38].
The solution in the first-order factorizing case is given by iterative integrals over a certain alphabet A = {f 1 (x), . . . f m (x)} together with special constants. We will present the algorithm for solving these systems, which does not require a special choice of a basis for the MIs, like the case in [11].
As an illustration, we employ this method of integration for computing the set of MIs which contribute to both the color-planar and complete light quark non-singlet three-loop contributions to the heavy-quark form factors for different currents, namely the vector, axial-vector, scalar and pseudo-scalar currents. The massive form factors for vector and axial-vector currents play an important role in the forward-backward asymmetry of bottom or top quark pair production at electron-positron and hadron colliders. The scalar and pseudo-scalar ones contribute to the decay of a Higgs boson to a pair of heavy quarks. They are also of importance to scrutinize the properties of the top quark [39,40] during the high luminosity phase of the LHC [41] and experimental precision studies at future high energy e + e − colliders [42]. The perturbative QCD contributions to these massive form factors at two loops were first computed in [43][44][45][46]. Later an independent computation was performed in [47] for the vector form factors, additionally including O(ε) terms in the dimensional parameter ε = (4 − D)/2. Recently, the two-loop contributions up to O(ε 2 ) for all the massive form factors were obtained in [48]. At three-loop level, the color-planar contributions to the vector form factors have been computed in [49,50] and the complete light quark contributions in [51]. Using the method described in this paper, we have obtained both the color-planar and complete light quark contributions to the three-loop form factors for the other three currents, namely axial-vector, scalar and pseudo-scalar currents in [52]. In a parallel and independent computation in [53] the same results have been obtained. The asymptotic behaviour of the heavy quark form factors has been studied in [54,55] recently, see also Refs. [47]. The large β 0 limit for massive form factors has been considered in [56] and in [57], where the three-loop scalar and pseudo-scalar form factors were computed in the static limit.
The paper is organized as follows. In Section 2 we describe the algorithm to solve first-order factorizing single-variate differential equation systems and present an illustrative example. In Section 3 we consider the massive three-loop vector form factors in an arbitrary basis and present the corresponding analytic results in Section 4. In Section 5 a numerical representation for the cyclotomic harmonic polylogarithms (HPLs) up to weight w = 6 is given, to allow the numerical evaluation of the massive three-loop form factors. Section 6 contains the conclusions. The complete expressions for the vector form factors, which are very large, are given in ancillary files together with the code CPOLY.f for the cyclotomic harmonic polylogarithms and other material, attached to this paper.

Description of the method
We consider n master integrals (MIs) I = (I 1 , . . . , I n ) which belong to the same topology and are functions of the dimensional parameter d = (4 − 2ε) and the variable x Here q 2 denotes the virtuality of the current and m is the heavy quark mass. One obtains an n × n system of coupled linear differential equations by taking the derivative for x of each of the MIs followed by the IBP reduction, Here the n×n matrix M consists of entries from the rational function field K(d, x) (or equivalently from K(ε, x)) where K is a field of characteristic 0; in the examples below the entries are even from Q(d, x) (or equivalently from Q(ε, x)). Furthermore, the inhomogeneous part R = (R 1 , . . . , R n ) is composed by simpler master integrals whose evaluations are immediate or can be carried out by other methods, like symbolic summation and integration; see [12] for details and references therein. In simpler situations R turns out to be just the 0-vector. For more involved applications we will assume that each entry R i is expanded into a Laurent series 2 in ε up to a certain order in terms of special functions. More precisely, we assume that the first coefficients R (j) i are given as polynomial expressions with coefficients from K in terms hyperexponential functions and iterative integrals over such functions; for a detailed definition see below. Furthermore, we assume that the unknown integrals I i can be expanded in an ε-expansion This applies to the case that M has no poles in ε. If this is not the case, according index-shifts have to be performed. Here it may happen that part of the equations which have to be solved, are no differential equations but are algebraic. Given such a coupled system, we seek for the first coefficients I (j) i in the form of polynomial expressions in terms of hyperexponential functions and iterative integrals over such functions.
. Such a function may be given in the form for some properly chosen l ∈ K. An iterative integral over hyperexponential functions is an integral of the form where f 1 (x), . . . , f λ (x) are hyperexponential functions and the lower bounds l 0 , . . . , l λ−1 ∈ K are appropriately chosen. The class of hyperexponential functions covers functions of the form q(x) µ where q(x) ∈ K(x) and µ ∈ K. Note that in all our calculations that arose so far, we only dealt with the special case µ ∈ Q. In the following we will use the property that f (x) g(x), 1 f (x) with f = 0 and d dx f (x) are hyperexponential functions provided that f (x) and g(x) are hyperexponential functions. In addition, d dx acting on the iterative integral (4) simply removes the outermost integral. As a consequence, applying the derivative to a polynomial expression in terms of hyperexponential functions and iterative integrals over such functions will lead again to a polynomial expression in terms of such functions. Furthermore, the multiplication of two iterative integrals over hyperexponential functions can be written as a linear combination of iterative integrals over hyperexponential functions due to its shuffle algebra [58]. Consequently also a polynomial expression in terms of hyperexponential functions and iterative integrals over hyperexponential functions can be always written as a linear combination of the form h 1 (x)I 1 (x) + . . . h λ (x)I λ (x) where the I i (x) are iterative integrals over hyperexponential functions and the h i (x) are hyperexponential functions.
A general assumption of our method will be that the degree of uncoupling will be of first order. More precisely, we will apply internally Zürcher's algorithm [59][60][61] implemented in the package OreSys [62] in order to decompose the system into one scalar linear differential equation (sometimes also several such equations) determining all unknown functions. In the case that the scalar equations (evaluated at ε = 0) are first-order factorizable, we proceed.
In general, the dimension n of the system (2) is rather high (e.g., n = 100). In this matter we note that the MIs can be distinguished sector-wise. A sector is defined by a set of maximum number of non-vanishing propagators in a single Feynman graph. Correspondingly, the absence of some propagators defines sub-sectors. The differential equation of a MI hence only contains integrals from the same sector or its sub-sectors. Thus, organizing the integrals in a way such that integrals with a minimum number of propagators are kept at the end of the list, provides an upper-block-triangular form of M, i.e. the diagonal elements of M are square matrices of not only rank one but higher. Each such square matrix represents a completely coupled set of integrals and we call them sub-systems of M. The advantage of arranging the system in this way is that now we can solve the system in a bottom-up approach, i.e. we first solve for the last set of coupled integrals in the list which depend on themselves only (plus inhomogeneous parts from R that are already expanded in terms of special functions), and then solve the second last set of coupled integrals, which depends on themselves and the last integrals and thus going up in the list.
We will now elaborate the different steps of our proposed algorithm. Similar ideas have been utilized already in Refs. [12,[63][64][65] in order to find solutions in terms of iterative sums over hypergeometric products. 1. Let us consider m integralsĨ = (Ĩ 1 , . . . ,Ĩ m ) which constitute a coupled sub-system, where the non-diagonal elements ofM are mostly non-zero and are rational functions from K(d, x), or equivalently from K(ε, x). In particular, we may assume thatM is an invertible matrix; if not, one can derive an alternative system by simple row operations with this property (the new system consists of less unknown integrals and the redundant integrals, that are removed from the system, can be expressed trivially by the integrals that arise in the new system). The inhomogeneityR is formed by contributions from integrals belonging to sub-sectors and the components of R. By construction we succeeded already in calculating the first coefficients of the ε-expansions of these integrals in terms of iterative integrals over hyperexponential functions. Consequently, plugging these results into R yields the ε-expansions are given explicitly in terms of iterative integrals over hyperexponential functions. Now we exploit the fact that for a certain topology and kinematics, the order k of highest pole of an integral inĨ is well-defined, as e.g. the integrals arising in three-loop massive form factors can have at most a pole of 1/ε 3 . Hence, one has the following Laurent expansions In order to determine the first coefficientsĨ i in terms of iterative integrals, we proceed as follows. We plug in (3) with undetermined coefficientsĨ i into (5), perform the series expansion in ε and consider the coefficient of ε k : for k = −l, −l + 1, etc.

2.
At each order in the ε-expansion we have now functions of a single variable x only. Solving order by order, one obtainsĨ as a Laurent series expansion in ε. To accomplish that we start with the coefficient of the leading pole ε −l . The corresponding sub-system is To solve Eq. (8), a natural first step is to reduce this m × m system to a higher order differential equation for a single integral. We will refer to this procedure as 'uncoupling' from now on. By using the package OreSys one obtains Here p l (x) are rational functions in K(x) and for some integer λ. Since the differentiation of iterative integrals over hyperexponential functions yields again iterative integrals over hyperexponential functions, the inhomogeneous part r(x) can be given explicitly in terms of iterative integrals over hyperexponential functions. Besides this scalar differential equation, the package OreSys provides in addition the solutionsĨ in terms of linear combinations ofĨ (−l) 1 (x) and its derivatives: with a k,i ∈ K(x). Like the r(x) in (10) the ρ k (x) can be given in such a form. Consequently, also the ρ k (x) can be expressed explicitly in terms of iterative integrals over hyperexponential functions. In other words, if one succeeds in solving the linear differential equations (9) and obtains the solution forĨ (−l) 1 (x) in terms of iterative integrals over hyperexponential functions, one can plug this closed form into (11) and can extract such an integral representation of the remaining functionsĨ . Concerning the uncoupling we remark that it is not advised using the classical cyclic vector algorithm to achieve the uncoupling, since generally this method provides uncoupled equations with large coefficients. Moreover, it is beneficial that Zürcher's algorithm may find several linear differential equations for several of the unknown functions: they have usually smaller orders than the cyclic vector algorithm (which always finds only one differential equation). The solving tools are now applied to each of the found equations. For simplicity we assume in the following that only one scalar differential equation forĨ 3. From Eq. (7) it is evident that the homogeneous solutions are always the same for any order in ε. The inhomogeneous solutions are different, however, by action of the inhomogeneities. We first consider the homogeneous solutions of Eq. (9). First we check if the differential equation can be factorized into first-order factors of the form withp k being rational functions in K(x) by using algorithms from [66][67][68]; for more details see [69,Chapter 4]. If this is possible, we proceed as follows. Define for 1 ≤ k ≤ m the hyperexponential functions for some appropriate lower bounds l k ∈ K, that are solutions of the kth first-order factor, i.e., Finally, one can read off from the factorization (12) the m solutions where the lower bounds l 0 , . . . , l m−2 are chosen accordingly. These solutions, also called d'Alembertian solutions [70], form iterative integrals over hyperexponential functions h k (x) h k−1 (x) . Since they are linearly independent over K, see [70,Thm. 5], yield the full solution space of the homogeneous recurrence (12). In the calculations presented below the hyperexponential functions h k (x) can be simplified all to rational functions in K. Even more is true after further simplifications: The integrands can be decomposed into the form by partial fractioning with φ l (x) = α l (x) β l (x) e l where the β l (x) are irreducible polynomials in K[x] and α l (x) are polynomials in K[x] with deg(α l (x)) < deg(β l (x)) and e l ∈ N\{0}.
Let us consider a typical example. Usually the functionsp k (x) in Eq. (12) are rational functions, which factor into the letters f l (x) of an alphabet A . If these letters are the ones of the Kummer-Poincaré type [71] the ratios h m /h m−1 in Eqs. (15) are again Kummer-Poincaré letters after partial fractioning. This representation holds as well for cyclotomic harmonic polylogarithms, since these have complex representations by Kummer-Poincaré letters.
By linearity one can now apply the integration sign to each of the summands in (16) and eliminate algebraic relations among the arising integrals utilizing their shuffle relations; these ideas have been elaborated in detail for the sum case [58,72]. In particular, the multiplicity e l can be reduced upon noting Likewise, the cyclotomic letters integrate to structures like In the course of these simplifications also products of these functions and corresponding harmonic polylogarithms arise that can be joined using shuffle relations [58]. In addition to the cyclotomic harmonic polylogarithms also their corresponding values at x = 1 contribute through partial integration. In the case of the harmonic polylogarithms these are the multiple zeta values (MZVs) [73]. In the cyclotomic case they are given by the cyclotomic constants [74][75][76][77]. In [78] relations beyond those known from [74][75][76][77] already have been conjectured using PSLQ [79]. More generally, this tactic can be applied in combination with the Almkvist-Zeilberger algorithm [80] if in addition hyperexponential functions arise that cannot be handled by the simplifications described above. In summary, a homogeneous linear differential equation stored in the variable de in terms of the unknown function f (x) can be solved in terms of d'Alembertian solutions by executing the HarmonicSums command In particular, if a factorization of the form (12) exists for the given differential operator, it will be computed and the full set of solutions (15) will be produced where all the simplifications described above are applied. In all our applications so far, the homogeneous solutions y i (x), i = 1, . . . , m could be expressed in terms of iterative integrals where f b (y) ∈ A are rational functions (or roots of rational functions) taken from a finite alphabet A. For instance, for the massive three-loop form factor discussed below, the alphabet can be chosen by where f b (x) corresponds to the bth entry. Summarizing, in our concrete application below the arising d'Alembertian solutions (15) will be simplified to expressions in terms of the class of harmonic polylogarithms [81] and the cyclotomic harmonic polylogarithms [77]. 4. The solution of the inhomogeneous differential equation (9) can be given explicitly by the following iterative integral [70] g Consequently,Ĩ where the constants C i are implied by (physical) boundary conditions and they are usually determined by separate calculations. Since the inhomogeneous part r(x) can be given in terms of iterative integrals over hyperexponential functions, also g(x) and thusĨ (−l) 1 (x) can be expressed in terms of iterative integral over hyperexponential functions. Furthermore, using our simplification tools from above, these integrals can be simplified further. E.g., within all our calculations we end up at alphabets of the form (21) or variants involving also rooted letters. We want to emphasize an alternative approach to find a particular solution g(x) of (9). If is the Wronskian of the linear differential equation (9) and then for some appropriately chosen l ∈ K yields another particular solution. Note that by (a mild generalization) of Abel's theorem we have that W (x) itself can be written as a hyperexponential function for some constant c ∈ K and an appropriately chosen lower bound l ∈ K; the polynomials come from the linear differential equation (9). Furthermore, the W i (x) are given by polynomial expressions in terms of the homogeneous solutions y i (x). As a conse- in (26) forms a polynomial expression in terms of hyperexponential functions and iterative integrals over such functions. In particular, g(x) yields such a representation. The following extra bonus often makes the formula (26) superior to (4): By reusing the simplified homogeneous solutions y 1 (x), . . . , y m (x) for (26), it is much easier to obtain a simplification of (26) than of (4) in terms of iterative integrals of the form (20) with alphabets like (21). (11) for k = 2, . . . , m. Since the derivation of iterative integrals over hyperexponential functions yields again iterative integrals over hyperexponential functions, all entries in the vectorĨ (−l) 1 (x) can be given within the class of iterative integrals over hyperexponential functions. 6. Finally, we plug this representation ofĨ (−l) (x) in terms of iterative integrals into (7) for k = −l + 1 and obtain a new system of the form (8)

Now we plug this representation ofĨ
). Thus we repeat the game for ε −l+1 and the remaining coefficients in (3) by induction/recursion. We note once more that the formula (23) remains the same, except that in (26) the function r(x) changes. As a consequence one can again reuse the already simplified homogeneous solutions y 1 (x), . . . , y m (x) and just needs to simplify g(x) in (26) with the updated function r(x). Let us illustrate the above algorithm by an example, which concerns the solution of a sub-system in the calculation of the three loop massive form factors in the color planar limit.

Example.
We consider the following system of differential equations for the integrals {J 1 , J 2 , J 3 } ∈ I: where c ij 's are rational functions in d, resp. ε, and x as given by The functions R i (ε, x) contain the inhomogeneous contributions from sub-sectors. They can be expanded into a Laurent series expansion in ε up to the required order and read Here we use the convention H a (x) ≡ H a and ζ l = ∞ k=1 1/k l , l ∈ N, l ≥ 2 denote the values of Riemann's ζ-function. The harmonic polylogarithms [81] are defined by and the letters f c are The HPLs are dual, by the Mellin transform, to the harmonic sums [82,83]. The solutions J i are calculated in terms of the following expansion in ε First, one obtains the determining equation for the leading pole O(1/ε 3 ): Using the incomplete Zürcher algorithm, one of four algorithms implemented in OreSys, we obtain for Eq. (35) two uncoupled differential equations, one of order two for J 3 (x) and another of first order for J 2 (x). J 1 (x) can directly be obtained from the solution for J 3 (x). The second order differential equation for J 3 (x) with the inhomogeneous part needs to be solved first. The differential operator in Eq. (36) can be written in factorized form In this case the rational functionsp i (x) are linear combinations of the letters spanning the HPLs. One now uses the method of the variation of the constants to obtain the solutions for the differential equation. The homogeneous solutions y 1 (x), y 2 (x) are given by and we obtain the solution where the constants C i can be determined from the physical boundary conditions and are known from a separate calculation. With the Wronskian W given by the integrals in (42) are easily evaluated For the remaining constants we find and thus The solution for integral J −3 1 (x) can directly be obtained from this result With these results at hand we can obtain the first order differential equation for where the inhomogeneous part is given by The homogenous solution is given by and thus we can obtain the full solution by evaluating The integral is easily evaluated and after fixing the constant of integration we obtain the final result To summarize, the results obtained so far, the solutions J are given by the numbers: Now, once the sub-system is solved for the highest pole, we consider the next order in ε. By construction, the homogeneous structure of the sub-system remains the same for any order in ε and hence the uncoupling procedure. The only change that takes place for different orders in ε is in the inhomogeneous parts which also constitute the contributions from the already-known previous orders in ε-expansion. Thus, for higher orders in ε, the only step is to iterate Eqs. (42) and (51). In this way we obtain the following solutions for the functions J 1 , J 2 , J 3 up to O(ε 0 ). Here we use the basis for the harmonic polylogarithms defined in [84].
In total the solution of Eq. (27)

Application to the heavy quark form factors
We consider the decay of a virtual massive vector boson of momentum q into a pair of heavy quarks of mass m, momenta q 1 and q 2 and color c and d, through the vertex Γ µ V,cd . We follow the notation used in Ref. [48]. Here q 2 = (q 1 + q 2 ) 2 is the center of mass energy squared. The general form of the amplitude is given bȳ whereū c (q 1 ) and v d (q 2 ) are the bi-spinors of the quark and the anti-quark, respectively, σ µν = The form factors are obtained from the amplitudes by multiplying appropriate projectors as provided in [48] and performing the trace over the color and spinor indices. n l and n h are the numbers of light and heavy quarks. For convenience, we use the Landau variable [86] x = q 2 − 4m 2 − q 2

Details of the calculation
The calculation of the three-loop massive vector form factors proceeds in a similar way as has been outlined in Refs. [48,52]. As usual the packages QGRAF [87], Color [88], Q2e/Exp [89,90] and FORM [91,92] have been used to generate the Feynman diagrams, calculate their color and Dirac-structure, and to determine moments for comparisons. The reduction to master integrals has been performed using Crusher [8]. Finally, we have obtained 109 MIs, out of which 96 appear in the color-planar case, as indicated in [52].

Figure 1: The color-planar topologies
In the color-planar limit, the families of integrals can be represented by eight topologies, shown in Figure 1, whereas for the complete light quark contributions, three more topologies are required, cf. Figure 2. Note that, only sub-topologies with a maximum of eight propagators contribute in the latter scenario.

Ultraviolet renormalization and infrared structure
The UV renormalization of the form factors has been performed in a mixed scheme. The heavy quark mass and wave function have been renormalized in the on-shell (OS) renormalization scheme, while the strong coupling constant is renormalized in the MS scheme, where we set the universal factor S ε = exp(−ε(γ E −ln(4π)) for each loop order to one at the end of the calculation. The required renormalization constants are available and are denoted by Z m,OS [101][102][103][104][105], Z 2,OS [101][102][103]106] and Z as [107][108][109][110][111][112][113] for the heavy quark mass, wave function and strong coupling constant, respectively. The renormalization of the heavy-quark wave function and the strong coupling constant are multiplicative, while the renormalization of massive fermion lines has been taken care of by properly considering the counter terms.
Considering the high energy limit, the universal behavior of infrared (IR) singularities of the massive form factors was first investigated in [114]. Later in [115], a general argument was provided to factorize the IR singularities as a multiplicative renormalization constant. Its structure is constrained by the renormalization group equation (RGE), as follows, where F fin I is finite as ε → 0. The RGE for Z(µ) reads d d ln µ ln Z(ε, x, m, µ) = −Γ(x, m, µ).
Here Γ is the corresponding cusp anomalous dimension, which is by now available up to threeloop order [116,117]. Note that Z is independent of the current. Both Z and Γ can be expanded in a perturbative series in α s as follows and the solution for Eq. (71) up to O(α 3 s ) is

The Three-Loop Vector Form Factors
We apply our algorithm to the single scale and first order factorizable system of differential equations which are relevant for the color-planar and the complete light quark contributions to the heavy quark form factors. One obtains the solutions for all contributing integrals in Laurent series expansion up to the required order in ε . Finally, using the results for the integrals, we obtain the color-planar and complete light quark (n l ) non-singlet contributions of the three-loop massive form factors for the vector current. Due to the substantial length of the expressions, we provide them as supplemental material along with this publication. In the following we only present expansions of the form factors in different kinematic limits and give numerical results for the whole kinematic region. Here the following abbreviation is used c 1 = 12ζ 2 ln 2 (2) + ln 4 (2) + 24Li 4 as mentioned in [48] and Li k (x) denotes the polylogarithm [118,119]. In V,1 ) is given by Here , T F = 1/2 and N C denotes the number of colors for SU (N C ) with N C = 3 in case of QCD. The magnetic component of the vector form factor (F V,2 ) is the anomalous magnetic moment of a heavy quark in this limit allowing for a cross-check of our computation with Ref. [122]. The form factor F (3) V,2 up to order y 2 reads

High energy region x → 0
Here we present the expansion of the form factors F V,1 and F V,2 in the asymptotic limit i.e. for x → 0 + up to O(x 2 ). The abbreviation L has been used to indicate ln(x).
In this limit, the electric component of the vector form factor i.e. F V,1 satisfies the Sudakov evolution equation. This behavior has been studied in detail in [54,55,114] accounting for the components known. The complete three-loop result has been given in [54] and partial four-loop results in [54,55].

Checks
By maintaining the gauge parameter ξ to first order, a partial check on gauge invariance has been obtained. After appropriately considering α s -decoupling, the UV renormalized results satisfy the universal IR structure, confirming again the correctness of all pole terms, see [54]. Finally, we have compared our results with those of Ref. [50,51], in the region x ∈ [0, 1] which have been obtained using a different method, and agree by adjusting the respective conventions. We also agree now with the results in [50,51] for the expansions given in [125].

Numerical Results
The color planar parts to the three-loop vector form factors F

Numerical implementation for Harmonic and Cylotomic Harmonic Polylogarithms
The color-planar part of the three-loop massive form factors depends on 206 cyclotomic harmonic polylogarithms (HPLs) [77] up to weight w=6 and correspondingly the harmonic polylogarithms also up to weight w = 6. In intermediary results for both cases HPLs of w = 8 appear. A FORTRAN-implementation of the harmonic polylogarithms to w = 8 has been given in [84]. 4 The space of the cyclotomic HPLs already up to w = 6 is very large and therefore we will rather represent the contributing individual functions numerically, and do not refer to an associated basis representation. The main argument of the cyclotomic HPLs, x, is located in the interval [−1, 1] in the present physical application. This is going to be the range we are considering in the following. In Ref. [77] the range x ∈ [0, 1] was considered. Here the cyclotomic HPLs are real-valued. In the extension to x ∈ [−1, 0[ some of the cyclotomic HPLs will become complex, as we will show below. The cyclotomic HPLs are given as iterated integrals over the letters (84) and those present in the usual HPLs, cf. (33). In the following {6, 0} and {6, 1} (resp. {3, 0} and {3, 1}) encode the corresponding cyclotomic letters. One obtains e.g.
The following cyclotomic HPLs contribute:  On the other hand, the former one is limited to double precision, while the latter one can be extended to arbitrary precision. We have tested the numerical implementation of the real-valued cyclotomic HPLs in CPOLY.f comparing to the corresponding results obtained by corresponding numerical results provided by Ginac. The code is compiled by gfortran CPOLY.f. The representation has an accuracy of ∼ 2 · 10 −15 (97) and better. The reading of the data needed in CPOLY.f requires 5.8 · 10 −2 sec. The calculation of all 206 cyclotomic HPLs at a given value of x is performed in 2.3 · 10 −3 sec or faster. In using the code CPOLY.f only the subroutines UCPOLYIN and UCPOLY are user routines to provide further input and to perform the calculation, respectively. Any use of the code CPOLY.f requires to quote the present paper. 5

Conclusion
We presented an algorithm to solve single-variate systems of differential equations, factorizing at first order and depending on the dimensional parameter ε, analytically. Here no choice of a special basis representation is required. The Laurent expansion in the parameter ε leads to a one-variable problem. We considered differential equations with rational coefficients in x and ε. The algorithm solves these systems in terms of iterative integrals over finite alphabets and rational terms to any order in the dimensional parameter ε. This method can be applied to a wide range of problems in Quantum Field Theory, after one knows whether the corresponding systems factorize to first order, which is checked by the present algorithm.
In the example of the massive three-loop form factors the emerging letters are those forming the HPLs and the cyclotomic HPLs at cyclotomy c = 3,4 and 6. The homogeneous solutions are the same for any order in ε. The corresponding inhomogeneities then determine the respective inhomogeneous solutions using the variation of constants. The iterative-integral structure is preserved by the latter operation, as can be shown by integration-by-parts. Besides the harmonic and cyclotomic harmonic polylogarithms up to weight w = 6 also associated special constants appear. In the cyclotomic case not all their relations have been proven yet by analytic methods. However, a series of relations has been conjectured by using PSLQ [78]. Assuming that these relations would hold, the results at three-loop order presented in this paper can finally be expressed by very few multiple zeta values only ln(2), ζ 2 , ζ 3 , Li 4 1 2 , ζ 5 (98) and no special cyclotomic constants contribute. However, cyclotomic constants remain in the expansion around x = −1.
Our result for the vector form factors agree with those given in Ref. [49][50][51]. We provide the FORTRAN-code CPOLY.f which allows to calculate the cyclotomic harmonic polylogarithms contributing to all massive three-loop form factors in the color-planar limit.