Renormalization of non-singlet quark operator matrix elements for off-forward hard scattering

We calculate non-singlet quark operator matrix elements of deep-inelastic scattering in the chiral limit including operators with total derivatives. This extends previous calculations with zero-momentum transfer through the operator vertex which provides the well-known anomalous dimensions for the evolution of parton distributions, as well as calculations in off-forward kinematics utilizing conformal symmetry. Non-vanishing momentum-flow through the operator vertex leads to mixing with total derivative operators under renormalization. In the limit of a large number of quark flavors $n_f$ and for low moments in full QCD, we determine the anomalous dimension matrix to fifth order in the perturbative expansion in the strong coupling $\alpha_s$ in the $\overline{\mbox{MS}}$-scheme. We exploit consistency relations for the anomalous dimension matrix which follow from the renormalization structure of the operators, combined with a direct calculation of the relevant diagrams up to fourth order.


Introduction
Within the gauge theory of the strong interaction, quantum chromodynamics (QCD), important nonperturbative information about the hadron structure is obtained from matrix elements of local operators between states with the same or different momenta. Depending on the momentum transfer, such operator matrix elements (OMEs) are used to describe the parton distribution functions (PDFs) in forward kinematics, their off-forward counter-parts, the generalized parton distributions (GPDs) or distribution amplitudes (DAs) for vacuum-to-hadron transitions, such as the pion form factor. Extensive reviews on the description of hadron structure using GPDs and DAs can be found in e.g. [1] and [2]. PDFs or GPDs can be extracted from data for hard inclusive and semi-inclusive reactions with identified particles in the final state, collected for example in deep-inelastic scattering (DIS) in the past at the HERA collider [3,4] or in the future through the research program at the planned Electron Ion Collider (EIC) [5,6]. Various DAs can be obtained from data taken in high-intensity, medium energy experiments, for instance at Belle II [7]. All such nonperturbative quantities are also accessible from first principles in QCD simulations of the hadron structure on the lattice, see e.g. [8,9] and for recent progress [10][11][12][13][14][15].
The scale-dependence of such distributions is governed by the renormalization group equations for the corresponding operators and can be computed order by order in the strong coupling α s in perturbative QCD. We focus here on the non-singlet quark OMEs relevant in DIS in the chiral limit and include mixing with operators involving total derivatives. In a given basis of local operators this implies renormalization with a triangular mixing matrix where the diagonal entries are the forward anomalous dimensions, well-known from inclusive DIS, but the off-diagonal elements of the renormalization matrix require a separate calculation.
The anomalous dimensions of non-singlet quark operators as functions of the Mellin moment N, which coincides with the Lorentz spin of the operator, are completely known up to three loops [16][17][18], and, likewise the corresponding splitting functions in x-space. At the four-loop level, low fixed moments [19][20][21] and complete results in the limit of a large number n f of quark flavors [22,23] have been obtained. More recently, fixed moments up to N = 16 have been computed and used to determine complete analytic results for a general SU(n c ) gauge theory in the planar limit, i.e. in the limit of a large n c [24]. Partial information consisting of the large n f limit [22] and of low moments [25] is available even at five loops.
In off-forward kinematics, the three-loop evolution kernel for flavor-non-singlet operators is known as well [26]. The computation has exploited conformal symmetry [27] of the QCD Lagrangian which was first utilized in pioneering work for the two-loop radiative corrections [28,29]. In d = 4 − 2ε dimensions and adopting the modified minimal subtraction (MS) scheme conformal symmetry in QCD is exact for large n f at the critical coupling. In the physical four-dimensional theory the renormalization group equations then inherit a conformal symmetry so that the generators of the conformal transformations and the evolution kernel commute [30]. Consistency relations following from the conformal algebra allow to restore the l-loop triangular mixing matrix for the off-forward anomalous dimensions from the l-loop forward anomalous dimensions and an (l − 1)-loop result for a so-called conformal anomaly [26,31,32]. In moment space, the offforward anomalous dimension depends on the Lorentz spin N as well as on the number k of total derivatives acting on the operator. In momentum space, the corresponding conjugate variables are the momentum fraction x carried by the struck parton and an additional kinematical variable such as the skewedness parameter of the process.
Fixed moments of non-singlet quark OMEs at three-loop order have also been computed in the MS renormalization scheme as well as in alternative ones, such as the regularization invariant (RI) scheme, which are suitable for a direct application to available lattice results [33,34]. The choice for a basis of the renormalized local operators in these works is different from the one in the conformal approach [28], so that a comparison between the fixed moment results of [33,34] and those of [26,28,29] is not immediate, requiring additional computational steps.
In the present article, we review the renormalization of non-singlet quark operators with particular emphasis on possible choices for their bases. In this way, we connect to different results that have appeared as of yet unrelated in the literature. We perform explicit computations of the relevant OMEs up to four loops for non-zero momentum transfer through the operator vertex and derive a number of consistency relations for the respective anomalous dimensions which govern the mixing associated with total derivative operators. This checks and extends previous calculations for the leading-n f terms in the evolution of flavor non-singlet operators in off-forward kinematics up to five loops.
The article is organized as follows: In Sec. 2 we set the stage, review the different basis choices for spin-N local non-singlet quark operators used in the literature and discuss their renormalization together with particular properties of the anomalous dimensions. Details of the computation are given in Sec. 3 and results for the moments of the complete mixing matrix up to five loops for the leading-n f terms are presented in Sec. 4. In Sec. 5, we list some results beyond the leading-n f limit, including the second order mixing matrix in the planar limit and fixed numerical moments in full QCD up to fifth order. The corresponding moments for a general gauge group SU(n c ) are deferred to App. A. Conclusions and an outlook are given in Sec. 6.

Theoretical framework
We consider the renormalization of spin-N local non-singlet quark operators where ψ represents the quark field, D µ = ∂ µ − ig s A µ the standard QCD covariant derivative, and λ α the generators of the flavor group SU(n f ). Since we focus on the leading-twist contributions, we symmetrize the Lorentz indices and take the traceless part, indicated by S. This projects the twist-two contribution, see e.g. [35].
The anomalous dimensions of interest are determined by considering spin-averaged OMEs obtained by inserting the respective operators in two-point functions, with quarks and anti-quarks of momenta p 1 and p 2 as external fields and all momenta are incoming, Choosing p 3 = 0 leads to zero momentum-flow through the operator and the OMEs are renormalized with the standard forward anomalous dimensions, see e.g. [24]. For the general case p 3 = 0 there is a momentum transfer through the operator vertex, which implies mixing between the operators in Eq. (2.1) and additional total derivative operators.
Next we will introduce three sets of bases commonly used in the literature. One is built from an expansion in Gegenbauer polynomials, which is used in the conformal approach, and the other two are based on counting powers of (total) derivatives.

The Gegenbauer basis
We start our discussion by introducing the renormalized non-local light-ray operators [O]. These act as generating functions for local operators, see [26], as Here, n is an arbitrary light-like vector. For simplicity, the x-dependence will be omitted in the following, writing O(z 1 , z 2 ) ≡ O(0;z 1 , z 2 ). The renormalization group equation for these light-ray operators can be written as with µ the renormalization scale, a s = α s /(4π) the strong coupling and β(a s ) the QCD betafunction. The evolution operator H (a s ) is an integral operator and acts on the light-cone coordinates of the fields [36] H with z α 12 ≡ z 1 (1 − α) + z 2 α and the evolution kernel h(α, β). The moments of the evolution kernel correspond to the anomalous dimensions of the local operators in Eq. (2.3) 6) where N represents the total number of covariant derivatives appearing in the operator, i.e. N = m + k in Eq. (2.3). The forward anomalous dimensions will be represented by γ N,N , and there is a shift N → N + 1 compared to the literature which is related to the operator definition.
The light-ray operators in Eq. (2.3) admit an expansion in a basis of local operators in terms of Gegenbauer polynomials, see e.g. [29,37], where k ≥ N is the total number of derivatives and we use the superscript G to denote the operators in the Gegenbauer basis. For any given N, the operator of lowest dimension, O G N,N , is a conformal operator. The Gegenbauer polynomials can be written in terms of the hypergeometric function 2 F 1 as [38] where the series representation in terms of the Gamma function will be more convenient for our purposes. The Gegenbauer polynomial with the differential operators then gives obey the evolution equation where the mixing of the operators is manifest. The anomalous dimension matrix 1 , denoted byγ G N , is triangular, i.e. its elements γ G N, j = 0 if j > N. It can be computed in QCD perturbation theory and its diagonal elements correspond to the standard forward anomalous dimensions γ N,N [24]. Furthermore, the superscript G can be dropped for γ N,N since they do not depend on the particular choice for the basis of additional total derivatives. The matrixγ G N is currently known up to the three-loop level [26] and we will discuss these results in more detail in Sec. 4.

The total derivative basis
Another basis for the quark operators, which directly generalizes Eq.
see e.g. [39,40] and [35,41]. Here the superscript D indicates that the operators are written in the total derivative basis and the indices p, q and r count the powers of the respective derivatives. In 1 We represent the elements of the mixing matrix for spin-(N + 1) operators as γ N,k and the matrix itself asγ N+1 .
practical calculations, all Lorentz indices are contracted with a light-like vector ∆, such that the indices p, q and r in Eq. (2.11) identify a given operator uniquely. Because of the chiral limit, the partial derivatives act as For the bare operators this defines a recursion, which is solved by Another consequence of the chiral limit is that left and right derivative operators renormalize with the same renormalization constants The anomalous dimensions γ D N,k governing the scale dependence of these operators are derived from the Z-factors by The mixing matrix is also triangular in the total derivative basis (γ D N,k = 0 if k > N) and, as was the case for the Gegenbauer basis, the diagonal elements are the standard forward anomalous dimensions γ N,N [24], where we have also dropped the superscript D due to the basis independence.
It is possible to relate the two operator bases using the light-ray operators in Eq. (2.3) as generating functions. Starting from the bare operators we have The evolution equation for the operators O G N,k in Eq. (2.10) relates the anomalous dimension matrices in the two bases Using Eq. (2.18) for renormalized operators the left-hand side can be expanded in the operators of the total derivative basis, 20) and upon comparing the coefficients of the operators [O D N−l,0,l ], we find To get the expression on the left-hand side, we have evaluated the sum

The Geyer basis
The local operators introduced by B. Geyer in [35,41] are expressed in a basis different from the one in Eq. (2.14), but the quoted anomalous dimensions can be related to γ D N,k in Eq. (2.16). The basis for the local operators used in those references is with N, k odd. The superscript B indicates that the operators are written in the Geyer basis. The contraction with an arbitrary light-like vector is understood, i.e.
We now want to relate the operators, and the corresponding anomalous dimensions, in the B-basis to those in the derivative basis. For this we use The corresponding relation for the renormalized operators reads Like in the Gegenbauer basis, we now focus on the diagonal operators O B N,N , the evolution equation for which is with N odd. Like in the other operator bases, the anomalous dimension matrixγ B N is triangular (γ B N,k = 0 if k > N) and in analogy to Eq. (2.21) γ B N,k can be related to the mixing matrices in the Gegenbauer and the total derivative bases. Details will be given below in Sec. 4.1.

Constraints on the anomalous dimensions
The elements of the mixing matrices for the operators in the total derivative basis are not all independent, but subject to particular constraints, which define useful relationships between them in the chiral limit. Starting from Eq. (2.12) we can derive the following relation by acting N times with a partial derivative on the bare operator O D Upon using the renormalization equations (2.14), this leads to a connection between the renormalized operators has to vanish for each value of k, this becomes a relation between the renormalization constants ∀k : which is equivalent to a relation between sums of elements of the mixing matrixγ D N , ∀k : This is a general statement, valid to all orders in a s . It can be considered as a consequence of parity conservation. If one takes a basis with ∂ z 1 − ∂ z 2 as in Eq. (2.7), i.e., "left minus right (l-r)" derivatives, each such derivative changes parity so that the operators with an even number of "l-r" derivatives do not mix with operators with an odd number of "l-r" derivatives. 2 By putting k = N − 1 we can relate the next-to-diagonal elements of the mixing matrix to the forward anomalous dimensions γ N,N , cf. Eq. (2.6) 3 , where, as above, the superscript D can be omitted for γ N,N .
The case k = 0 in Eq. (2.33) gives which relates the sum of the elements in the N-th row of the mixing matrix to the conjugate C γ D N,0 of the last column defined as see e.g. [42].
With the help of Eq. (2.35) for the last column of the mixing matrix another relation between the anomalous dimensions in the total derivative basis and the Gegenbauer one can be obtained. Substituting Eq. (2.35) into Eq. (2.21) one finds Considering now arbitrary k in Eq. (2.33) results in This is a useful relation for several reasons. Its derivation only relies on the chiral limit, which imposes constraints on the renormalization structure of the operators. It can, therefore, be used as an order-independent consistency check, which any expression for γ D N,k has to obey. Alternatively, Eq. (2.38) provides a path for the construction of the full mixing matrixγ D N from the knowledge of the forward anomalous dimensions γ N,N and the last column γ D N,0 at any order of perturbation theory.
Thus, with partial information being available even to five-loop order, one can construct an ansatz for the off-diagonal elements and use Eq. (2.38) to test its self-consistency. In this case, γ D N,0 serves as boundary condition and Eq. (2.34) as a cross-check. This leads to the following 4-step algorithm for the construction of the mixing matrix: 1. Starting from the forward anomalous dimensions γ N,N and the bare OMEs in Eq. (2.2), one determines the all-N expressions for the next-to-diagonal γ D N,N−1 and the last column entries γ D N,0 of the mixing matrix. In the next section, it will be detailed how to relate the last column entries to the calculation of Feynman diagrams.
2. Next, one calculates the sum Based on the structure of the result, one can then make an ansatz for the off-diagonal elements.

One calculates the double sum
with the chosen ansatz and collects everything into Eq. (2.38). This leads to a system of equations in the unknown coefficients of the ansatz. It should be mentioned that Eq. (2.38) by itself is not enough to determine a unique solution for γ D N,k . In particular, some structures might be self-conjugate, like e.g.
Uniqueness of the solution is saved however by the fact that there is a boundary condition, namely the expression for γ D N,k has to agree with the previously found expression for γ D N,0 from step 1.
4. If the system of equations can be solved consistently, one has determined the final expression for the off-diagonal terms of the mixing matrixγ D N . If not, the structure of the remaining terms in Eq. (2.38) can be used to adapt the ansatz, leading one back to step 3. This approach will be applied and described in further detail below.

Calculation in Mellin N-space
We will describe the computation of Feynman diagrams for the matrix elements of operators in Eq. (2.1) along with the necessary details for the renormalization. The constraints for the anomalous dimensions in Eq. (2.33) and the constructive algorithm to determine the mixing matrix lead to particular types of symbolic sums and techniques for their evaluation will be reviewed as well.

Calculating Feynman diagrams
We calculate the spin-averaged OMEs defined in Eq. (2.2). To ensure tracelessness and symmetry of the Lorentz indices, the OMEs are contracted with a tensor of light-like ∆, and ∆ 2 = 0. OMEs in off-forward kinematics require momentum-flow through the operator vertex, p 3 = 0, see Fig. 1. For simplicity but without loss of generality, one can choose p 2 = 0, so that the calculated OMEs are of the form With this momentum configuration the OMEs in Eq. (2.2), which are a priori three-point functions reduce effectively to two-point ones, i.e. one is left with massless propagator-type Feynman diagrams. These can be efficiently calculated using computer algebra methods as will be detailed next.
The computations follow a well-established workflow. The Feynman diagrams are generated using QGRAF [43], the output of which is then directed to a FORM [44,45] program to determine the topologies and to compute the color factors of the diagrams, the latter part based on the algorithms presented in [46]. The necessary Feynman rules are given for example in [24]. For computational efficiency, diagrams of the same topology and color factor are grouped together in so-called meta-diagrams, cf. [47] for further details on the use of FORM in Feynman diagram calculations. The handling of the meta-diagrams is done with the database program MINOS [48]. The actual diagram calculations are then performed using the FORCER program [49], which can efficiently deal with massless propagator-type diagrams in d = 4 − 2ε dimensional regularization [50,51] up to four loops. In this way we can obtain fixed moments of the OMEs in Eq. (3.2).
Finally, the calculations are done using a general covariant gauge. Although the operators considered are gauge invariant, the OMEs will in general depend on the gauge parameter. Using a general covariant gauge, the independence of the anomalous dimensions of the gauge parameter then provides a check on our calculations.
For the renormalization we use the MS-scheme [52,53]. In this scheme, the evolution of the strong coupling is governed by with β(a s ) the standard QCD beta-function and β 0 = (11/3)C A − (2/3)n f . Writing out Eq. (2.14), the operators are then renormalized as where the quark wave function renormalization factor Z ψ takes care of the self-energy corrections for the off-shell external quarks and we have abbreviated, cf. Eq. (3.1), (3.5) In this notation, O 1 simply represents the vector current Since this is a conserved quantity, its anomalous dimension vanishes and Z 0,0 = 1. For the full set of operators up to spin (N + 1), the renormalization takes the form of a matrix equation . . .
Note that the bottom-right n × n submatrix (with n < N) represents the full mixing matrix for the renormalization of the spin-n operators. For example, the 2 × 2 submatrix of Eq. (3.7) is exactly the matrix appearing in the renormalization of the set of spin-2 operators {O 2 , ∂O 1 }.
The anomalous dimensions γ D N,k emerge from the Z-factors according to Eq. (2.16) and can be expanded in a power series in a s The explicit form of the Z-factors in a perturbative expansion in dimensional regularization in the MS-scheme can be obtained from Eq. (2.16) in terms of the anomalous dimensions and the coefficients of the QCD beta-function. For illustration we quote them up to order a 3 s , but the expansions can readily be generalized to higher orders. We present separately the diagonal factors Z N,N , which are just those that renormalize the forward operators and the off-diagonal ones Z N,k , where k = N is understood.
The mixing under renormalization is manifest in the appearance of sums over anomalous dimensions, starting at order a 2 s in the expression for Z N,k .
The computation of the bare OMEs for the operators O N+1 in Eq.
The computation of the bare OMEs uses fixed moments and Eq. (2.35) can be rewritten as a relation for the last column γ D N,0 , which illustrates the bootstrap in N. For a given moment N, the element γ D N,0 of the last column is expressed in terms of the elements γ D j,0 with j < N and the sum of the elements in the N-th row of the mixing matrix. Using in addition Eq. (2.36) it is straightforward to rewrite this as a recursion

Calculating sums
At the l-loop level, the forward anomalous dimensions in general consist of harmonic sums of maximum weight w = 2l − 1 and denominators in N + α (with α ∈ N) up to the same maximum power. More specifically, the maximum weight of a specific term in the anomalous dimension depends on both its color-structure and on the number of powers of n f multiplying it. Terms multiplied by the color factors C A = n c and C F = (n 2 c − 1)/(2n c ) can have weight up to 2l − 1, while each additional factor of n f will decrease this maximum weight by one, down to weight l for the leading-n f terms. A similar reasoning holds for terms multiplied by values of the Riemann-zeta function ζ n defined as (3.15) The maximum weight of a term containing ζ n is in general 2l − 1 − n, and l − n for the leading-n f terms specifically. Harmonic sums at argument N are recursively defined by [42,54] S ±m (N) = The constructive use of the conjugation in Eq. (2.38) requires the evaluation of single and double sums over such structures. These sums are non-standard in the sense that they are outside the class of sums that can be solved by the algorithms encoded in the SUMMER program [42] in FORM, which has been a standard in the calculation of the forward anomalous dimensions. However, the type of sums appearing can be dealt with using the principles of symbolic summation, see e.g. [55,56] for extensive overviews. In particular, the MATHEMATICA package SIGMA [57] is very helpful. For this reason, we briefly review aspects of (creative) telescoping and the key features of SIGMA.
Suppose one wants to find a closed form for some summation (3.17) Often it is possible to solve this problem using the telescoping algorithm. The task is then to find a function g(N) such that the summand can be written as Here ∆ represents the finite-difference operator. Whenever this is possible, the summation problem can be solved as . was constructed by Zeilberger and is called creative telescoping [58]. Here, the task is to find functions c 0 (n), . . . , c d (n) and g(n, k) such that g(n, k + 1) − g(n, k) = c 0 (n) f (n, k) + · · · + c d (n) f (n + d, k) . Summing both sides of the equation then gives, applying telescoping to the left-hand side Hence one gets an inhomogeneous recurrence for the original sum of the form q(n) = c 0 (n)S(n) + · · · + c d (n)S(n + d) . (3.23) The creative telescoping algorithm is the main method used in SIGMA to solve summation problems. After the recurrence is generated, SIGMA first solves this equation in terms of the solution of the homogeneous version and a particular solution. The final closed form expression for the initial sum is then given as a linear combination of the solutions of the recurrence that has the same initial values as the sum.
4 Results up to five loops in the leading-n f limit Here we explore the large-n f limit of the mixing matrices up to five loops using the consistency relations for the anomalous dimension matrices discussed above along with explicit results for the forward anomalous dimensions [22,23] and direct computations of the relevant OMEs in Eq. (3.2) up to four loops. Up to the three-loop level, the results for γ G N,k in the Gegenbauer basis are known [26]. In the total derivative basis, only fixed low-N results for γ D N,k are available [33,34].

One-loop anomalous dimensions
At the one-loop level, the mixing matrixγ G N in the Gegenbauer basis is diagonal [37,59], i.e. 5 γ G,(0) This implies that the operators with spin N = 5 just renormalize multiplicatively at the one-loop level with 6γ G,(0) In contrast, in the total derivative basis the off-diagonal elements are non-zero already at this order. Application of the algorithm described in Sec. 2.4 leads to γ D,(0) which is consistent with the result in the Gegenbauer basis according to Eq. (2.37). We also find agreement with the fixed moments presented in [33]. For spin-5 operators in the total derivative basis, the mixing matrix is then At this order, another check can be made, namely by comparing with the anomalous dimensions in the Geyer basis [35,41]. The off-diagonal piece of the mixing matrix at one-loop is [35,41] . (4.5) Using Eq. (2.28) for k = N we find the following relation between the anomalous dimensions in the different operator bases where the + -sign holds for even values of N and the − -sign for odd N. We have checked that our result Eq. (4.3) obeys these relations. For spin-5 operators in the Geyer basis, the one-loop mixing matrix is thenγ

Two-loop anomalous dimensions
Beyond the one-loop level, also the mixing matrix in the Gegenbauer basis gets non-zero offdiagonal contributions. These can be calculated in general as [26,28]  The discrete step-function in the last term is defined as ϑ N,k ≡ 1 if N − k > 0 and even 0 else. (4.12) The N − k even condition originates from the fact that, in the Gegenbauer basis, only CP-even operators are considered. The conformal anomaly,ŵ(a s ), can be written as a power series in the strong couplingŵ (a s ) = a sŵ (0) + a 2

sŵ
(1) + . . . . For the determination of the mixing matrix at order a l s , the conformal anomaly is only needed up to order a l−1 s [31]. At the two-loop level, the general relation Eq. (4.8) leads to [26] γ G,(1) (4.14) For the leading-n f contributions, only the term proportional to β 0 is relevant. The anomalous dimension can also be written in terms of harmonic sums and denominators, giving for the leadingn f piece For the spin-5 operators, the leading-n f piece of the mixing matrix is then In the total derivative basis we find

Three-loop anomalous dimensions
The Gegenbauer anomalous dimensions, using again Eq. (4.8), arê An analytic expression for the two-loop conformal anomalyŵ (1) depending on N and k in the space of Mellin moments is currently not available [32]. Also, in the literature, e.g. [26], only implicit relations of the form Eq. (4.19) are given. The expansions in terms of harmonic sums here and at higher orders in the following sections are new.
Again, the leading-n f piece again just comes from the term proportional to β 0 . In terms of harmonic sums we find and explicitly for the spin-5 operatorŝ γ G,(2) In the total derivative basis our method gives

Four-loop anomalous dimensions
We start by presenting the results in the total derivative basis. There are two contributions to the anomalous dimensions, namely terms with and without a factor of ζ 3 . Since ζ 3 already has weight-3, the structure multiplying it will just be weight-1. Hence in complexity the contribution of such a term, using our algorithm, is equivalent to a one-loop calculation. We find and for the spin-5 operatorŝ The expression for the ζ 3 -independent terms is Next we turn to the results in the Gegenbauer basis. These cannot be found in the literature, but can be calculated in the same way as the lower-order results. Since the one-loop beta-function and leading-n f three-loop forward anomalous dimensions do not depend on ζ 3 , there are no ζ 3 -terms in the off-diagonal part of the Gegenbauer mixing matrix 28) which again parallels the one-loop case. Hence for N = 5 we just havê

Discussion of the results in the total derivative basis
From the structure of the anomalous dimensions γ D N,k presented here, we can make general predictions for the ζ-independent terms of γ D,(l) N,k .
1. The prefactor of 1 (N+2) l can be written as where we take 2. The prefactor of the difference S l−1 (N) − S l−1 (k) can be written as 3. The prefactor of powers of [S 1 (N) − S 1 (k)] l−1 can be written as 5. We can also say something about powers of denominators. If at the l-loop level we have a structure 1 (N + α) δ (4.37) then the l + 1 expression will have the term −β 0 (N + α) δ+1 , (4.38) where α = 1, 2 and δ = l.
The above considerations imply that the l-loop expression for the leading-n f off-diagonal piece of the mixing matrix only has a small number unknowns, namely the coefficients of the weight w = 1 terms and [S α (N) − S α (k)]. Also at higher loops we could expect terms of the form [S α (N) − S α (k)] β with α, β > 1. The other higher-weight terms can be reconstructed from the lower-loop expressions. Hence, by considering just a few low (N, k)-pairs, it should be possible to use the conjugation relation, Eq. (2.38), to completely fix the off-diagonal part of the mixing matrix. When combined with the all-order expression for the leading-n f terms of the forward anomalous dimensions, which was calculated in [22] using the critical exponents of the Wilson operators, this implies that the ζ-independent terms of the leading-n f piece of the full mixing matrix can in principle be reconstructed to all orders.

Five-loop anomalous dimensions in the total derivative basis
As an application of the discussion presented in the previous section, we determine the part of the leading-n f terms of the five-loop mixing matrix, which is independent of Riemann zeta-values. For this, the five-loop expression for γ N,N presented in [22] is also used. In total there are six terms with a priori unknown numerical prefactors, namely These can be fixed by using Eq. (2.38) for (N, k)-pairs up to N = 4, and we quickly find (4.40) We have checked that the resulting expression obeys Eq. (2.38) also for high moments. This then suggests the following pattern for the terms of the form [S α (N) − S α (k)] at l loops: whereγ denotes that we only use the denominator-structure itself, i.e. we do not include the overall numerical factor of the lower-loop expression. The five-loop N = 5 mixing matrix iŝ γ D,(4) 5 Going beyond the leading-n f -limit So far, the focus has been to apply our method to determine the mixing matrix in the limit of large n f . However, the algorithm can also be applied beyond this limit, and we will present two possible threads of this. In the first the two-loop mixing matrix is determined in the planar limit, i.e. large n c . In the second we present some low-N matrices in full QCD.

Second order mixing matrix in the leading color limit
The application of our method to the determination of the leading-n f part of the mixing matrix works particularly well because of the simple structure of these terms. This simplicity manifests itself in three ways: (1) only harmonic sums with positive indices appear, (2) the number of possible terms is restricted because of the low maximum weight and, (3) increasing the order in a s by one corresponds to increasing the maximum weight of the possible structures by one. Moving away from the leading-n f limit and towards the leading color limit, point (1) still remains valid. Hence, while the sums that need to be evaluated according to our algorithm become more numerous, they are of the same type as those in the leading-n f limit, and the method goes through in exactly the same way.
As an application of this, we present the off-diagonal piece of the second order mixing matrix in the planar limit. The two-loop anomalous dimensions can be written as The first term just represents the leading-n f piece. The leading color term, γ D,(1) N,k | lc , is proportional to C F n c /2 and the next-to-leading color one, γ D,(1) N,k | nlc , to C F /(2n c ). Using our method we then find γ D,(1) For illustration, we quote the corresponding mixing matrix for spin-5 operatorŝ The algorithm should also be applicable to the subleading color part γ D,(1) N,k | nlc . However, because of the non-planar Feynman diagrams contributing, harmonic sums with negative indices appear in the expression for the anomalous dimension. This adds an additional layer of complexity to the problem and we leave this part to future studies.

N = 5 matrix in full QCD
Here we present the mixing matrix for spin-5 operators in full QCD 8 , i.e. C F = 4/3 and C A = 3. The corresponding expressions for a general gauge group SU(n c ) can be found in Appendix A. Before presenting the explicit results, the relevant relations used in their derivation are reviewed.

Relevant relations
As explained in Sec. 2.4, the next-to-diagonal elements can be calculated directly from the forward anomalous dimensions using Eq. (2.34). This already fixes γ 1,0 , γ 2,1 , γ 3,2 and γ 4,3 . For the last column, the relation between the Gegenbauer basis and the total derivative one, Eq. (2.37), can be used to determine γ 2,0 , γ 3,0 and γ 4,0 . Explicitly At N = 5, there is not enough information to completely fix all the anomalous dimensions. Nevertheless, we can again use Eq. (2.35) to derive These relations can be used to determine the complete N = 2 matrix up to the five-loop level. To the same order we can also calculate γ 2,1 . Finally, with the three-loop information in the Gegenbauer basis from [26], the full N = 4 matrix can be calculated up to third order. For the fixed-N forward anomalous dimensions we use the four-loop expressions from [19][20][21]24] and the five-loops ones from [25]. At the three-loop level, only the leading-n f term of γ

Conclusion and outlook
We have studied the renormalization of non-singlet quark operators appearing in deep-inelastic scattering, including total derivative operators. The mixing of operators under renormalization and the explicit form of the anomalous dimension matrix depends on the choice of a basis for the operators involving total derivatives. We have discussed three choices, which have appeared in the literature, one in terms of Gegenbauer polynomials and two based on counting (total) derivatives. We have shown how to transform the anomalous dimension matrices from one basis to the other, thereby connecting as of yet unrelated results in the literature and using the results in the different bases for mutual cross-checks up to the three-loop level. This provides highly non-trivial tests of the respective computations, which have applied entirely different techniques, for one exploiting conformal symmetry and a dedicated determination of the conformal anomaly, or, otherwise, evaluating Feynman diagrams of the respective OMEs directly. The result for the three-loop evolution kernel for flavor-non-singlet operators in off-forward kinematics is thus established by independent methods. The quoted analytic expressions in moment space in the Gegenbauer and the total derivative basis in terms of harmonic sums are new.
Moreover, we have exploited consistency relations between the off-forward anomalous dimensions in the chiral limit to derive new results at the four-and five-loop level for the anomalous dimension matrices, presenting the leading-n f terms of the mixing matrix of non-singlet quark operators as well as the full QCD result for some low-N fixed moments in the MS renormalization scheme. With an additional scheme transformation to the RI scheme, these will also become useful in studies of the hadron structure with lattice QCD approaches.
The algorithms presented in this paper, i.e. consistency relations combined with direct Feynman diagram computations of the respective OMEs, are expected to work also beyond the leadingn f or the planar limit discussed here. They allow for an automation of the calculations with the help of various computer algebra programs, such as FORCER under FORM for the computation of massless two-point functions up to four loops or the symbolic summation with the help of SIGMA. Furthermore, the approach can be adapted to the calculation of mixing matrices between other types of operator in QCD, like flavor-singlet operators, or to different models all-together, like gradient operators in scalar field theories. We leave these aspects for future studies.

Acknowledgements
We thank V. Braun, J. Gracey and A. Manashov for useful discussions and for comments on the manuscript. The Feynman diagrams have been drawn with the packages AXODRAW [60] and JAXODRAW [61].
This work has been supported by Deutsche Forschungsgemeinschaft (DFG) through the Research Unit FOR 2926, "Next Generation pQCD for Hadron Structure: Preparing for the EIC", project number 40824754 and DFG grant MO 1801/4-1.
Here we present the off-diagonal elements of the N = 5 mixing matrix for a general gauge group SU(n c ). Starting from the 4-loop level, the following color factors contribute The subscripts F and A denote the fundamental and adjoint representation, and N F = n c and N A = n c 2 − 1 represent their dimensions.