The Three-Loop Splitting Functions $P_{qg}^{(2)}$ and $P_{gg}^{(2, N_F)}$

We calculate the unpolarized twist-2 three-loop splitting functions $P_{qg}^{(2)}(x)$ and $P_{gg}^{(2,\rm N_F)}(x)$ and the associated anomalous dimensions using massive three-loop operator matrix elements. While we calculate $P_{gg}^{(2,\rm N_F)}(x)$ directly, $P_{qg}^{(2)}(x)$ is computed from 1200 even moments, without any structural prejudice, using a hierarchy of recurrences obtained for the corresponding operator matrix element. The largest recurrence to be solved is of order 12 and degree 191. We confirm results in the foregoing literature.


Introduction
The three-loop anomalous dimensions and splitting functions P (2) ij (x) govern the scale-evolution of the parton distribution functions in Quantum Chromodynamics (QCD) at next-to-next-toleading order (NNLO) and are of importance for precision predictions at ep-and hadron colliders, such as HERA, the Tevatron and the LHC. They are instrumental for the precise prediction of the Higgs-, top-quark, and jet-production cross sections, concerning both, the parton distribution functions [1,2] and the strong coupling constant α s (M Z ) [3], thus allowing the exploration of the hitherto heaviest systems in the Standard Model in greatest possible detail. High precision predictions are also of importance for planned ep-facilities as EIC and LHeC [4], operating at high luminosity. A first computation of the three-loop splitting functions was performed in Refs. [5,6], after the calculation of fixed moments for these quantities in Refs. [7][8][9]. In these references the anomalous dimensions were calculated using massless processes.
In the present paper we compute the splitting function P (2) qg and the N F -dependent part of P (2) gg independently using massive three-loop operator matrix elements (OMEs). They are necessary for the computation of the heavy flavor contributions to deep-inelastic scattering in the region of virtualities Q 2 much larger than the heavy quark mass squared m 2 . While we obtain the complete expression for P (2) qg , only the N F -part of P (2) gg emerges, since the respective process is proportional to the color factor T F . 1 We calculate P (2,N F ) gg directly by applying the techniques described in Ref. [10]. In the case of P (2) qg , we apply the recently proposed method of large Mellin moments [11]. In previous papers we have obtained the corresponding contributions ∝ T F for the splitting functions P (2,N F ) gq , P (2,NS,N F ) qq , P (2,NS,TR,N F ) qq [12,13] and the complete splitting function P (2),PS qq [14]. A series of Mellin moments for these quantities has been obtained before in Ref. [9]. Our results agree with those of the previous calculation in Refs. [5,6]. Recently, parts of the three-loop anomalous dimensions ∝ T F have also been obtained for a series of anomalous dimensions in two-mass calculations [15].
The paper is organized as follows. In Section 2 we describe the calculation of the anomalous dimensions γ (2) qg and γ (2),N F gg in Mellin-N space. Section 3 contains the conclusions. In the Appendices we present the splitting functions in x-space and give some technical details in calculating master integrals. In all representations we project onto the algebraic basis [16], both for the harmonic sums [17,18] and the harmonic polylogarithms (HPLs) [19].
Furthermore, we use the convention The pole terms O(1/ε 3 ) and O(1/ε 2 ) contain the complete O(a s ) and O(a 2 s ) anomalous dimensions, and the single pole term allows to obtain the terms ∝ T F N F of the O(a 3 s ) anomalous dimensions γ (2) qg and γ (2) gg , which in the case of γ (2) qg is the complete anomalous dimension. The Feynman diagrams contributing to the OMEs A gg,Q (2.2) were generated by the code QGRAF [25] 2 . The color configurations were calculated using the package Color [26] and the Feynman integrals were reduced to master integrals using the package Reduze 2 [27,28]. 3 . There are different techniques available to calculate the master integrals, which have been summarized in Refs. [11] and [10].

The Anomalous Dimension γ (2)
qg In calculating γ (2) qg we have used the method of large moments [11]. First the local operator insertions (∆.p) N are formally resummed into propagator like terms by The system of differential equations in the variable x obeyed by the master integrals can be obtained by integration by parts identities through Reduze 2. The master integrals can then be expanded into a formal Laurent-Taylor series, which can be inserted in the differential equations, leading to associated difference equations, with known initial conditions. One can then exploit these difference equations to generate a large set of moments for each color and ζ-value coefficient, expressed as a sequence of rational numbers. In the present approach it is required that the master integrals have to be known to a higher order in the dimensional parameter ε. Usually standard techniques allow to derive the corresponding expansion. In one particular case we had to spend some effort to obtain the initial conditions for some master integrals, which is described in Appendix B.
With the corresponding series of rational numbers at hand, we use the algorithm of guessing [31] to derive their associated minimal recurrence. Using the package Sigma [32,33] this recurrence is then solved in terms of nested sums and the solutions are afterwards transformed to special function formats using the package HarmonicSums [34][35][36][37][38].
The recurrences derived for the anomalous dimension are first order factorizable and have therefore sum-and product solutions according to difference field theory [39][40][41][42][43][44][45][46]. Moreover, the solutions turn out to be given even in terms of harmonic sums only. To obtain a solution conform with a later analytic continuation one should consider even moments only, which would require a twice as large set of moments as input. In many problems this number is so high that it might be difficult to reach. Alternatively, one can start by taking all moments and derive a first expression at general values of N , from which only the even moments are then projected. This is algebraically possible by considering (−1) N as a variable y and performing division with remainder w.r.t. the polynomial 1 − y (or 1 + y). In this way we obtain the N th moment of the single pole term. The use of also odd moments is possible in other approaches and recommended if a lower number of moments allows to obtain the difference equation to be determined. It yet corresponds to a different problem in general. However, given a recurrence for odd and even values, one can now generate easily a very high number of moments using this recurrence, can extract afterwards the needed even moments and can repeat the guessing process with the now sufficiently large number of even moments.
In the calculation we computed only the irreducible contributions and added the other terms known already as a function of N , cf. Refs. [9,[20][21][22][23][24]. Here the master integrals are not 2 See Ref. [9] for the implementation of the local operators. 3 The package Reduze 2 uses the packages FERMAT [29] and Ginac [30]. solved in terms of sums, but only their moments are determined and those are inserted into the expressions of A (3) Qg , for which, again, only moments are calculated at first. We have then projected the irreducible contributions, using the known terms, on to γ (2) qg using Eq. (2.1). The order and degree of the recurrences associated to the different color and ζ-factors, as well the times to guess the corresponding recurrences using Guess [31] and to solve them using the packages Sigma [32,33] and HarmonicSums [34][35][36][37][38] are summarized in Table 1. We have generated 2000 even moments for the unrenormalized expressionÂ (3) qg up to its 1/ε term using the package SolveCoupledSystem in which the method of large moments is implemented [11]. The computation of the moments for the master integrals took 9.30 days per kernel, where we parallelized to 15 kernels. The creation of all pole terms inÂ (3) qg took 10.67 days/kernel where we parallelized to 25 kernels. The moments of the OME formed the input of the subsequent calculation. We used the even moments to guess the corresponding recurrences using Guess [31]. The implementation of the guesser in Sage [47,48] needs a much shorter time to find a recurrence than an earlier one written in Mathematica [31]. The total guessing time of ∼ 89 sec was smaller than the corresponding one using the Mathematica-version by a factor 75. The Sage-version of the guesser does currently not report the number of moments used. For the most demanding example, the C 2 A T F N F -term, 1000 moments were not sufficient, but 1200 moments were. In our earlier study [49], slightly different recurrences were obtained, since there we determined them from even and odd moments. The present recurrences are somewhat shorter.
The computation time for the last two steps to obtain γ (2) qg amounted to 44.6 min, which included the time to create the different recurrences through guessing and their solution to obtain the closed expressions as a function of the Mellin variable N .
This computational time and memory requirement are much shorter than that one using the traditional method, where first all master integrals are calculated in explicit form and are then inserted to obtain the final result. The costs to determine the different master integrals as functions of N would be much higher, if possible at all, because elliptic structures may arise at higher orders in ε.
In any case, the time needed to reduce the OME to master integrals is the longest part of the calculation. Using Reduze 2 [27,28] the reduction time amounted to several months. Similar or larger run times are expected also for other reduction programs.  The O(a s ) and O(a 2 s ) anomalous dimensions are read from the O(1/ε 3 ) and O(1/ε 2 ) terms and γ (2) qg is determined from the O(1/ε) term. The expressions are reduced to the polynomial basis in the harmonic sums [17,18] One obtains with the polynomials The anomalous dimensions γ (0,1) qg agree with results given in Refs. [50][51][52][53][54][55][56][57], or with their Mellin transform.
For the three-loop anomalous dimension γ (2) qg , we obtain (2.14) The anomalous dimension γ (2) qg , Eq. (2.11), agrees with the corresponding expression given in Ref. [6]. The terms ∝ N 2 F were confirmed before in [58]. gg,Q up to the 1/ε terms in the traditional way, i.e. the contributing master integrals are first calculated in explicit form. A subset of master integrals can be obtained using representations through (generalized) hypergeometric functions [59,60] or using Mellin-Barnes integrals [61,62]. The resulting definite sums have been simplified afterwards with symbolic summation using the packages Sigma [32,33], EvaluateMul-tiSums [63], and HarmonicSums [34][35][36][37][38]. A further large set is calculated using the method of differential equations [10,64]. Here we use the algorithm described in Ref. [10], which is available in SolveCoupledSystem using the packages Sigma, HarmonicSums, EvaluateMultiSums, SumProduction [63] and OreSys [65]. This algorithm allows to solve the differential equations analytically in the one-parameter case without the need to use any specific basis of master integrals. In the present calculation of anomalous dimensions, the corresponding difference equations are factorizing in first order and can be solved analytically using the package Sigma. In a single case the computation time turned out to be large. Therefore, we used the Almkvist-Zeilberger theorem [66] and our packages MultiIntegrate [35] to derive a linear recurrence. Solving the recurrence with Sigma and HarmonicSums provided then an immediate solution. The master integrals are obtained as expressions in the variable x, see Eq. (2.5). The solution in Mellin N -space is then given by the N th expansion coefficient of these expressions, which is found by a routine of HarmonicSums.
For the contribution ∝ N F to γ (2) gg at three-loop order, γ (2),N F gg , we obtain: where the polynomials are The anomalous dimension γ (2),N F gg , Eq. (2.55), agrees with the corresponding expression given in Ref. [6]. We have confirmed the leading N F -terms ∝ N 2 F in [67,68], which were also given in [69].

Conclusions
We have computed the 3-loop anomalous dimension γ (2) qg (N ), the contributions ∝ N F for γ (2) gg (N ) and the associated splitting functions in a massive calculation, which is fully independent of the earlier computation in Ref. [6]. We agree with the previous results. In the case of γ (2) qg (N ), we have performed a fully automatic calculation generating the Feynman diagrams, reducing them to master integrals and calculating them using differential and difference equations. We used formal Taylor series representations in terms of the subsidiary parameter x resumming the local operator insertions. The new method of high Mellin moments [11] has been used. Here the moments are calculated recursively using the difference equation systems associated to the differential equations. We did not compute the master integrals as explicit functions of N , but kept them as vectors of moments. The method of guessing has then been used to obtain one difference equation for each color/ζ-value factor. All difference equations obtained factorize at first order and one finally obtains a representation in harmonic sums.
The contributions ∝ N F to γ (2) gg (N ) have been calculated using more traditional methods. After the reduction to the master integrals, we have calculated these as functions of the Mellin variable N . The corresponding techniques are discussed in detail in Ref. [10] and include (generalized) hypergeometric and associated higher transcendental function representations [59,60], the method of hyperlogarithms [70][71][72], differential equation techniques [10,64], and integration using the Almkvist-Zeilberger theorem [35,66]. The final solution has then been obtained by inserting the master integrals and extracting γ (2),N F gg (N ) from the single pole term ofÂ (3) gg . In this process, quite a series of sums cancel, including all non-harmonic sums.
The universality of the QCD anomalous dimension allows to compute them within various setups. In the present calculation, they were obtained from the pole structure of massive OMEs. The calculation of these OMEs is part of an ongoing project with the final goal to compute the massive Wilson coefficients for deep-inelastic scattering in the region Q 2 ≥ m 2 .
The anomalous dimensions and splitting functions presented in this paper, together with our earlier results in Refs. [12-14, 58, 67, 68], are given in Mathematica.m format in the attachment, together with the expressions for I B5a 101 for D = 4 + ε to D = 30 + ε up to O(ε 6 ) presented in appendix B. Their conversion to maple or Form-inputs [73] is straightforward.

B Calculation of the initial values
A sufficient number of initial values has to be provided to solve the recurrence relations which determine the master integrals. These are the moments of the master integrals for low values of the Mellin variable N . Up to O(ε 0 ) they can be calculated using the MATAD [74] setup from [9,75]. Beyond this order, other methods have to be used. If a multi-sum representation has already been derived to the necessary order in ε, the required initial values can be derived from it. However, since this involves some degree of manual intervention and since an appropriate representation cannot always be found, a fully algorithmic method is desirable. In the following, we describe an algorithmic method which works for all master integrals arising in the ongoing project to calculate the massive OMEs to three-loop order. The method has been briefly described in [67] and applied in several calculations of this project [10,14,67]. It is based on the α-parameterization for Feynman integrals, see, e.g., [76][77][78], and the idea [79,80] of rewriting integrals with scalar products in the numerator in terms of scalar integrals in shifted dimensions with raised propagator powers. This allows us to map fixed moments of the master integrals with operator insertions to operator-less tadpole integrals. These can be reduced using integration-by-parts (IBP) identities [27,28,81] and lead to linear combinations of only three tadpole master integrals. Once these are calculated in shifted dimensions to sufficient order in ε, we can assemble the moments of the master integrals with operator insertions. The steps described here can be implemented efficiently in Mathematica.

B.1 Structure of the integrals
The scalar integrals with operator insertions for which we need initial values are two-point integrals with an on-shell, massless external momentum p 2 = 0. The Feynman rules of the operators involve a light-like vector ∆ and they introduce symbolic powers of scalar products ∆.k, where k is a linear combination of loop momenta and p. Due to the structure of the Feynman rules, we can deal with these scalar products elegantly by introducing an ordinary generating function for the Mellin moments of the integrals. In the simplest case of a scalar product (∆.k) N , we introduce a tracing variable x, multiply by x N and sum over N , cf. Eq. (2.5). Similar terms appear also for more complicated operator insertions. The operator insertions always lead to products of terms like (1−x ∆.k) −1 , which resemble propagators that are linear in the momentum k.
All scalar integrals of this project can be expressed in terms of 24 integral families with nine quadratic propagators (which depend quadratically on a linear combination of loop and external momenta) and three linear propagators that encode the structure of the operator insertion. We use the notation where P i denotes the propagator denominators of integral family f . The first nine indices refer to the quadratic propagators and the last three to the linear propagators. The 24 twelvepropagator integral families are based on six nine-propagator integral families (denoted B1, B3, B5, C1, C2 and C3), which are supplemented by different sets of linear propagators (indicated by an additional letter in the family identifier). For example, the integral family B5a has the propagators We first discuss the properties of the quadratic propagators. The corresponding ninepropagator families differ by the placement of the mass m and by the routing of the external momentum p. If no external linear propagator is present, they can only depend on the mass m and the external momentum p. Thus, they reduce to massive tadpole integrals, since p 2 = 0, and we can remove p from the propagators altogether. Then three quadratic propagators become linearly dependent of the remaining six (P 2 , P 4 and P 9 in the example above) and the six nine-propagator families reduce to only two six-propagator tadpole integral families.
Let us briefly recapitulate the basic facts about the α-parameterization and its use for expressing integrals with irreducible numerators in terms of scalar integrals in shifted dimensions. For more details, see for example [76][77][78][79][80]. The α-parameterization for the above integrals without linear propagators reads after applying a Wick rotation where we have explicitly pulled out a factor (−1) from each propagator P i,E due to the Wick rotation. Below we suppress the subscript E and work with Euclidean momenta unless explicitly stated otherwise. Each propagator is at most quadratic in the loop momenta. Hence, we can write any linear combination of propagators in the form The matrix A is real and symmetric, and the vectors b i will be linear combinations of external momenta. Performing the loop integrals leads to The elements of the matrix A are the coefficients of the scalar products of the loop momenta, while the elements of the vector b are the coefficients of the terms linear in the loop momenta -i.e. scalar products with the external momenta (or other vector quantities). The scalar term c finally contains the masses of the propagators. For our purposes, we have to consider the following form of the propagators: The δ coefficients are either 1, 0 or −1 depending on how these terms enter the propagator under consideration. Each propagator enters the exponential multiplied by its corresponding α parameter if it is present in the integral. We can model this by replacing α i with β i = α i Θ(ν i − 1 2 ), which vanishes if the propagator is not present. Thus, we obtain Here, we have already used that p 2 = 0. For the integral families of this project, the explicit form of the matrix A and the vector b depend only on whether the integral belongs to one of the families B1, B3 or B5 or to C1, C2 or C3. Using the notation introduced above, the quantities read The scalar term c is different for each of the six integral families. However, we skip it here since its specific form does not play a role in the following.

B.2 Treatment of the linear propagators
If the integrals contain scalar products involving the loop momenta in the numerator, we can include them in the α-parameterization by introducing an exponential generating function for them, i.e. We introduce one auxiliary parameter r for each scalar product. Thus, auxiliary parameters associated to scalar products of two loop momenta, k i .k j , end up as additional terms in the matrix elements A ij (e.g., r 1 k 1 .k 2 in the exponential leads to the replacement A 12 → A 12 + r 1 ). The auxiliary parameters which accompany scalar products with one loop momentum, p.k i , are added to the parameter b i (e.g., the scalar product r 2 p. To derive a parameter representation for fixed moments of the master integrals with operator insertions, we first derive the α-parameterization for the corresponding operator-less integral. Afterwards, we add in the linear propagators. Each of them has the structure (1−x ∆.k) −1 , where k is a linear combination of loop and external momenta. As discussed above, it is introduced as an ordinary generating function for the moments of the master integrals given in Eq. (2.5). If there are multiple linear propagators, all of them have the same tracing parameter x, such that the expansion involves finite sums arising from the Cauchy products. For example, two linear propagators yield which provides the generating function of the Feynman rule for an operator insertion on a vertex. Analogously, for three linear propagators, we obtain In case a linear propagator is raised to a higher power, we use Thus, given a number of linear propagators, each of them raised to a given power n i , i.e. (1 − x y i ) −n i , we can extract the coefficient of x N of the series expansion in x as follows: We introduce a summation index j i and a factor for each linear propagator, multiply them together and sum over j i from 0 to j i+1 . Of course, there is no sum over j 0 and the outermost index (which is never summed over) must be renamed to N .
To include the operators into our description, we think of the sums we just introduced as a fully expanded polynomial in the scalar products of the linear propagators. Each term has a number of scalar products each raised to a definite power and for each term, we use the exponential generating function representation, described in Eq. (B.10). To summarize: For each linear propagator we get a binomial coefficient which encodes the combinatorics of multiple, identical linear propagators, and we replace each power of a scalar product by a corresponding power of the derivative operator w.r.t. the auxiliary parameter r i (which belongs to the ith linear propagator).
Given all these ingredients, we can assemble a general formula for the master integrals with operator insertions for the case of three linear propagators Here, A is completely determined by the operator-less part of the integral (A B and A C from Eqs. (B.7), (B.8)). b has one term, proportional to p, which is contributed by the operator-less part (b B and b C from Eq. (B.9)) and one part, proportional to r∆, which stems from the linear propagators. We describe this part in the next paragraph below. The scalar term c contains the mass term, but it remains unchanged by the operators, unless the linear propagator involves the external momentum p.
The coefficients b of the linear term are modified by the introduction of the exponential generating functions for the operators: To each b i we have to add a term r j ∆ where the linear propagators contain the loop momentum k i . Let us clarify this with an example: For family B1a we have the linear propagators We introduce three parameters, r 1 , r 2 and r 3 , for the three propagators and arrive at the exponential generating function and thus we obtain for b Similar terms arise for the other integral families.
If the linear propagators contain the external momentum p, we have to also modify the scalar term c. In that case we have to add a term proportional to ∆.p. As an example, we use family B3b, which has the linear propagators and leads to the exponential generating function and hence The central idea of the dimensional shifts is now that derivatives with respect to r i will act on the exponential function and introduce factors of (A −1 ) ij = (det(A)) −1Ã ij , whereÃ ij is the (i, j) cofactor of the matrix A. Since the dimension D enters the parametric representation Eq. (B.15) only in the exponent of the determinant (except for constant prefactors), we can reinterpret the integrals arising from taking the derivatives as integrals in shifted dimensions, according to the exponent of the determinant. The cofactorsÃ ij are polynomials in the α parameters and can be absorbed into the powers of (α k ) ν k −1 which encode the powers to which the quadratic propagators are raised. More concretely, the powers of the propagators enter the α representation in the form which means that an additional factor of α k j in the integrand corresponds to raising the power of propagator P j by k (i.e., ν j = ν j + k) and multiplying the expression by, cf. e.g. [60], to compensate for the changed (−1) ν i and Γ(ν i ). Here (ν i ) k denotes the Pochhammer symbol. Thus, each derivative shifts the dimension and raises the powers of some propagators. After applying all derivatives, we obtain a linear combination of scalar integrals in shifted dimensions with dotted propagators and no operator insertions. These can then be mapped to tadpole integrals.
Let us now have a closer look at how the derivatives act on the parameter representation. For the master integrals we are concerned with, the auxiliary parameters r i only enter b and c. The vector b has two terms-a term proportional to p and one proportional to ∆. Since b only appears in the bilinear form b T A −1 b and since p 2 = 0 and ∆ 2 = 0, only the mixed terms proportional to ∆.p survive. Therefore, all terms are linear in the r i . Thus, each derivative acting on this part of the exponential brings down exactly one matrix element of A −1 and one power of ∆.p. This means that each derivative shifts the dimension of the integral from D to D + 2. For the scalar term c we also have a linear dependence on the r i . Here, however, we do not get any shift of dimension from a derivative acting on this part of the exponential, since it only brings down a power of 2∆.p. This leads to linear combinations of integrals with different dimensional shifts. In all cases, each derivative introduces a factor ∆.p, so that the N th moment of any master integral will always be proportional to (∆.p) N .
Finally, we have to undo the Wick rotation. The only scalar product left in the final expression is (∆.p) N which introduces an additional, global factor of (−1) N into (B.15).

B.3 Reduction of the tadpole integrals
The procedure described above maps a scalar integral with an operator insertion in D = 4 + ε dimensions to a (possibly large) linear combination of operator-less tadpole integrals in shifted dimensions. To cope with these, we use IBP relations derived with Reduze 2 to map them onto linear combinations of master integrals. Their coefficients are rational functions in the dimension D and the mass squared m 2 . The reductions are done for general values of D so that they can be used for integrals in any dimension. Particular care must be taken if the initial values are mapped to a linear combination of tadpole integrals in different dimensions.
As already observed in [82], all tadpole integrals can be reduced to just three master integrals. They correspond to and their diagrams are shown in Figure 1.

B.4 Required tadpole integrals
The first two tadpole integrals are straightforwardly obtained and can be expressed completely in terms of Γ-functions, for example using their Feynman parameterization. The first one is just the one-loop tadpole cubed, which reads These expressions can be expanded in ε in any dimension D = 2n + ε. When expanding in ε, we keep the normalization N , the mass scale m 3D−2ν and an additional factor of exp( 3 2 εγ E ) unexpanded to simplify the expressions. The latter is the same regardless of the integer dimension around which we expand. The third integral is the three-loop banana graph with fully massive lines and reads (B.30) In [82] this integral was found to correspond to an expression involving Γ-functions and a generalized hypergeometric function 3 F 2 for general values of D. Here, we are interested in the expansion around ε = D − 2n for different values of n ∈ N up to 2n = 30. To calculate it, we follow the method used in [83][84][85], which we briefly describe below. We calculate this integral as a special case of the same integral with two different masses, by deriving a differential equation in the ratio of the masses and then specializing to the case of equal masses. The two-mass diagram is depicted in Figure 2a. The inhomogeneous part r(x) of the differential equation arises from the expansion of T 1 (x) and T 2 (x) in ε as well as lower orders of the solution to J 1 (x). We then use the variation of constants to calculate a particular solution to the inhomogeneous differential equation, cf. [87], is the Wronskian of the homogeneous solutions. For D = 4 it is W (x) = x 2 (1 − x)(1 + x). The structure of these integrals tells us that they only lead to iterated integrals over the letters x, 1 − x and 1 + x, i.e. the harmonic polylogarithms [19]. Performing these integrals can be done efficiently using the HIntegrate command from the package HarmonicSums. Finally, we fix the constants of integration using the expansion of J 1 (0) = I B1a 112 along with the relation Eq. (B.32). Since only HPLs occur, the relation maps them to linear combinations of HPLs again. If both sides are algebraically reduced, we can solve for the constants by comparing coefficients of the HPLs. Of course, this requires the algebraic relations for HPLs of argument x to be available. In intermediate steps of the calculation at hand, we needed these relations up to weight w = 9, which can be calculated, e.g., using HarmonicSums, cf. also [17,88]. The expansion of the solution for J 1 (x) is given completely in terms of HPLs with polynomial coefficients in x. For D = 4 + ε we obtain where we suppressed the higher order terms in ε due to their length.
To obtain the solution for I B5a 101 = J 1 (1), we express the HPLs at x = 1 in terms of multiple zeta values and use their relations [88] to reduce the occurring constants to an algebraically independent basis. The full result for D = 4 + ε reads s 6 = σ −5,−1 , s 7a = σ −5,1,1 , s 7b = σ 5,−1,−1 , s 8a = σ 5,3 , s 8b = σ −7,−1 , s 8c = σ −5,−1,−1,−1 , s 8d = σ −5,−1,1,1 , where the σ a are multiple zeta values, defined by σ a = lim N →∞ S a (N ). Except for the zeta values ζ k and B 4 [82] the remaining multiple zeta values cancel in the final results for the operator matrix elements. Similar results are obtained for this integral in all even dimensions D = 4 + ε to D = 30 + ε up to and including the O(ε 6 ) terms. This allows the calculation of moments up to N = 13 for the master integrals with operator insertions. Since this exceeds the number of required initial values for the recurrences, we can use the remaining moments as a cross check for their solutions.