Three-loop evolution equation for flavor-nonsinglet operators in off-forward kinematics

Using the approach based on conformal symmetry we calculate the three-loop (NNLO) contribution to the evolution equation for flavor-nonsinglet leading twist operators in the $\overline{\text{MS}}$ scheme. The explicit expression for the three-loop kernel is derived for the corresponding light-ray operator in coordinate space. The expansion in local operators is performed and explicit results are given for the matrix of the anomalous dimensions for the operators up to seven covariant derivatives. The results are directly applicable to the renormalization of the pion light-cone distribution amplitude and flavor-nonsinglet generalized parton distributions.


Introduction
The remarkable progress in experimental techniques in the past two decades has provided a fresh impetus to the study of hard exclusive and semi-inclusive reactions with identified particles in the final state. Such processes are interesting as they allow one to access the hadron structure on a much more detailed level as compared to totally inclusive reactions. A (probably still distant) major goal is to understand the full three-dimensional proton structure by "holographic imaging" of quark and gluon distributions in transverse distance and momentum spaces. The related experiments have become a prominent part of the research program at all major existing and planned accelerator facilities, e.g. the Electron Ion Collider (EIC) [1].
The relevant nonperturbative input in such processes in many cases involves operator matrix elements between states with different momenta, dubbed generalized parton distributions (GPDs), or vacuum-to-hadron matrix elements related to light-front hadron wave functions at small transverse separations, the distribution amplitudes (DAs). The scale-dependence of such distributions is governed by the renormalization group (RG) equations for the corresponding operators, where, in contrast to standard parton densities, mixing with the operators involving total derivatives has to be taken into account. Going over to local operators one has to deal with a triangular mixing matrix where the diagonal entries are the anomalous dimensions, the same as in deep-inelastic scattering, but the off-diagonal elements require a separate calculation. the evolution kernel can be written as a sum of several contributions with a simpler structure. This construction is similar in spirit to the "conformal scheme" of Refs. [8][9][10]. We find that the remaining (canonically) SL(2)-invariant part of the three-loop kernel satisfies the reciprocity relation [17][18][19][20], discussed in Sect. 4. The explicit construction of the invariant kernel is presented in Sect. 5. We provide analytic expressions for the terms that correspond to the leading asymptotic behavior at small and large Bjorken x, and a simple parametrization for the remainder that has sufficient accuracy for all potential applications. In Sect. 6 we explain how our results for the renormalization of light-ray operators can be translated into anomalous dimension matrices for local operators. In this way also the formal relation to the results in [7][8][9] is established. The final Sect. 7 is reserved for a summary and outlook. The paper also contains several Appendices where we collect the analytic expressions for the kernels.

Evolution equations for light-ray operators
A renormalized light-ray operator, [O](x; z 1 , z 2 ) = ZO(x; z 1 , z 2 ) = Zq(x + z 1 n)/ nq(x + z 2 n), (2.1) where the Wilson line is implied between the quark fields on the light-cone, is defined as the generating function for renormalized local operators: where D µ = ∂ µ − igA µ is the covariant derivative. Here and below we use square brackets to denote renormalized composite operators (in a minimal subtraction scheme). Due to Poincare invariance in most situations one can put x = 0 without loss of generality; we will therefore often drop the x dependence and write O(z 1 , z 2 ) ≡ O(0; z 1 , z 2 ).
The renormalization factor Z is an integral operator in z 1 , z 2 which is given by a series in 1/ǫ, d = 4 − 2ǫ, with where Z = ZZ −2 q ; Z q is the quark wave function renormalization factor and γ q = M ∂ M ln Z q the quark anomalous dimension. The QCD β-function β(a) and γ q are known to O(a 5 ) [21][22][23][24][25].
It is easy to see [14] that translation-invariant polynomials (z 1 − z 2 ) N are eigenfunctions of the evolution kernel, The eigenvalues γ N (a) correspond to moments of the evolution kernel in the representation (2.9), They define the anomalous dimensions of leading-twist local operators in Eq. (2.2) where N = m+k is the total number of covariant derivatives acting either on the quark or the antiquark field. The corresponding mixing matrix in the Gegenbauer polynomial basis is constructed in Sect. 6. The leading-order (LO) result for the evolution kernel in this representation reads [26]: which satisfy the usual SL(2) algebra It can be shown that as a consequence of the commutation relations [H (1) , S α ] = 0 the corresponding kernel h (1) (α, β) is effectively a function of one variable τ called the conformal ratio [27] h (1) (α, β) =h(τ ) , τ = αβ αβ , (2.16) up to trivial terms ∼ δ(α)δ(β) that correspond to the unit operator. This function can easily be reconstructed from its moments (2.12), alias from the anomalous dimensions. Indeed, it is easy to verify that the result in Eq. (2.13) can be written in the following, remarkably simple form [27] h (1) where the regularized δ-function, δ + (τ ), is defined as Beyond the LO this property is lost. However, the evolution kernels for leading twist operators in minimal subtraction schemes retain exact conformal symmetry. Indeed, the renormalization factors for composite operators in this scheme do not depend on ǫ by construction. As a consequence, the anomalous dimension matrices in QCD in four dimensions are exactly the same as in QCD in d = 4 − 2ǫ dimensions that enjoys conformal symmetry for the specially chosen "critical" value of the coupling [14][15][16]. The precise statement is that the QCD evolution kernel H(a) commutes with three operators that satisfy the SL(2) algebra (2.15). These operators can be constructed as the generators of collinear conformal transformations in the 4 − 2ǫ-dimensional QCD at the critical point and have the following structure [14][15][16]: whereβ(a) is the QCD β-function (2.5) and S (0) α are the canonical generators (2.14). Note that the generator S − corresponds to translations along the light cone and does not receive any corrections as compared to its canonical expression, S (0) − . The generator S 0 corresponds to dilatations; its modification in interacting theory ∆S 0 = S 0 − S (0) 0 can be related to the evolution kernel H(a) from general considerations [14]. Finally S + is the generator of special conformal transformations and Eq. (2.20c) is the most general expression consistent with the commutation relations (2.15). To see that, note that ∆S + = S + − S  If H(a) is known, this equation can be used to find ∆(a) and in this way construct the SL (2) generators that commute with the evolution kernel in a theory with broken conformal symmetry. The main point is, however, that ∆(a) can be calculated independently from the analysis of the conformal Ward identity [7,13,16]. In this way Eq. (2.21) can be used to calculate the non-invariant part of the evolution kernel with respect to the canonical generators S (0) ±,0 . Indeed, expanding the kernels in a power series in coupling constant one obtains from (2.21) a nested set of equations [14] [S (0) is expressed in terms of the lower order kernels, H (k) and ∆ (k) with k ≤ ℓ − 1.
The first of them, Eq. (2.23a), is the usual statement that the LO evolution kernel commutes with canonical generators of the conformal transformation. As a consequence, the corresponding kernel h (1) (α, β) can be written as a function of the conformal ratio (2.16) and restored from the spectrum of LO anomalous dimensions. The result is presented in Eqs. (2.13), (2.17).
The second equation, Eq. (2.23b), is, technically, a first-order inhomogeneous differential equation on the NLO evolution kernel H (2) . To solve this equation one needs to find a particular solution with the given inhomogeneity (the expression on the r.h.s.), and add a solution of the homogeneous equation [S (0) + , H (2) ] such that the sum reproduces the known NLO anomalous dimensions. This calculation was done in Ref. [15] and the final expression for H (2) is reproduced in a somewhat different form below in App. A.
In this work we solve Eq. (2.23c) and in this way calculate the three-loop (NNLO) evolution kernel H (3) . The main input in this calculation is provided by the recent result for the two-loop conformal anomaly ∆ (2) [16]. Since the algebraic structure of the expressions at the three-loop level is quite complicated, we separate the calculation in several steps in order to disentangle contributions of different origin. The basic idea is to simplify the structure as much as possible by separating parts of the three-loop kernel that can be written as a product of simpler kernels.

Similarity transformation
The symmetry generators S α in a generic interacting theory (2.20) involve the evolution kernel H(a) and additional contributions ∆(a) due to the conformal anomaly. These two terms can be separated by a similarity transformation Note that H and H obviously have the same eigenvalues (anomalous dimensions). Going over to the "boldface" operators can be thought of as a change of the renormalization scheme, The "rotated" light-ray operator [O(z 1 , z 2 )] U satisfies the RG equation Looking for the operator U in the form U = e X , X(a) = aX (1) + a 2 X (2) + . . . , (3.4) we require that the "boldface" generators do not include conformal anomaly terms, 0 we can assume that X (k) commute with the same generators as well, whereas comparing Eqs. (2.20c) and (3.5c) yields the following set of equations: These equations can be used to fix X (1) and X (2) from the known one-and two-loop expressions for the conformal anomaly, ∆ (1) and ∆ (2) [15,16], up to SL(2) (canonically) invariant terms. In other words, the transformation U that brings the conformal generators to the form (3.5) is not unique; there is some freedom and we specify our choice later on. The one-loop result, X (1) , turns out to be rather simple whereas the two-loop operator, X (2) , is considerably more involved. Explicit expressions are presented in App. B. The rotated, "boldface" evolution kernels satisfy a simpler set of equations as compared to Eqs. (2.23), as the terms involving the conformal anomaly are removed, These equations are solved by inv , inv + T (1) β 1 + 1 2 H inv + T inv are (canonically) SL(2)-invariant operators with kernels that are functions of the conformal ratio (2.16) and the operators T (i) commute with S Similar to the X kernels defined as the solutions to Eqs. (3.7), the T kernels are fixed by Eqs. (3.10) up to SL(2) (canonically) invariant terms. Explicit expressions are collected in App. C. Note that the expressions for the perturbative expansion of the evolution kernel in Eq. (3.9) can be assembled in the following single expression: (3.11) Finally, adding the contributions from the rotation matrix U = exp{aX (1) + a 2 X (2) + . . .} we obtain the following results for the first three orders of the evolution kernel in the MS scheme: inv , I ] where all entries are known except for the SL(2)-invariant part of the three-loop kernel H inv that has yet to be determined. Explicit expressions for the X and T kernels are given in App. B and App. C, respectively. 1 The SL(2)-invariant kernels H (k) inv can be written in the following general form Here z α 12 is defined in Eq. (2.10), τ = αβ/(ᾱβ) and P 12 is the permutation operator, P 12 f (z 1 , z 2 ) = f (z 2 , z 1 ). Γ cusp is the cusp anomalous dimension which is known to the required accuracy [6]: 14) The LO expression (2.13) corresponds to The two-loop constant term χ (2) 0 and the functions χ inv (τ ) are given in App. A. The three-loop expressions will be derived below.
Note that the expression for the two-loop kernel in Eq. (3.12) differs from the one derived in Ref. [15] where the non-invariant part is written in form of a single expression. The representation of the non-invariant part of the two-and three-loop kernel as a product of simpler operators suggested here seems to be sufficient and probably more convenient for most applications.

Reciprocity relation
The eigenvalues of H (k) correspond to the flavor-nonsinglet anomalous dimensions in the MS scheme that are known to three-loop accuracy [6]. Note that our definition of the anomalous dimension γ (k) (N ) differs from the one used in Ref. [6] by an overall factor of two. In addition, in our work N refers to the number of derivatives whereas in [6] the anomalous dimensions are given as functions of Lorentz spin of the operator. Thus Ref. [6] .

(4.2)
One can show that the eigenvalues of the invariant kernels H with a "natural" choice of T operators specified in App. C satisfy the following symmetry relation: Let j N = N +2 (conformal spin); the asymptotic expansion of γ (k) inv (N ) at large j → ∞ only contains terms symmetric under the replacement j N → 1 − j N . Note that this symmetry does not hold for the anomalous dimensions γ (k) (N ) themselves.
The argument goes as follows. As well known [28], conformal symmetry implies that the evolution kernels can be expressed in terms of the quadratic Casimir operator of the symmetry group. As a consequence it is natural to parameterize the anomalous dimensions in the form It has been shown [18][19][20] that the asymptotic expansion of the function f (j) = af (1) (j)+a 2 f (2) (j)+ . . . at large j only contains terms that are symmetric under reflection j → 1 − j. In all known examples the f -function also proves to be simpler than the anomalous dimension itself. Expanding both sides of Eq. (4.4) in a power series in the coupling one obtains 2 so that the values of f (k) (j N ) are related to the anomalous dimensions at the same order of perturbation theory up to subtractions of certain lower-order terms. Let us compare this expansion with the relations between the eigenvalues of H (k) in Eq. (4.1) and of the invariant kernel H (k) inv in Eq. (4.3). Using Eq. (3.12) and explicit expressions for the eigenvalues of the T kernels in Eq. (C.6) one obtains (note that the commutator terms do not contribute to the spectrum) Comparing this expansion to the one in Eq. (4.5) we see that In other words, the QCD anomalous dimension in the MS scheme can be written in terms of the eigenvalues of the invariant kernel as 8) and the asymptotic expansion of γ inv (N ) at large N only contains terms invariant under the reci- [18][19][20]. This relation implies a certain condition for the choice of T kernels which appears to be natural in the first few orders of perturbation theory, see App. C.

Three-loop invariant kernel H (3) inv
The three-loop invariant kernel H (3) inv takes the form The three-loop cusp anomalous dimension Γ cusp is known (3.14) so that our task is to determine the constant χ  inv (τ ). This can be done by using the information on the spectrum of H can easily be calculated from the known one-, two-, and three-loop flavor-nonsinglet anomalous dimensions [6]. The constant term χ in the invariant kernel (5.1) corresponds to the constant term in the large-N asymptotic of the anomalous dimension and is straightforward to obtain. Using the expressions from Ref. [6] we find The calculation of the functions χ inv (τ ) and χ P (3) inv (τ ) is much more involved. As usual one has to consider even and odd values of N separately. We write inv (N ) correspond to the eigenvalues of the invariant kernel for even (odd) N . Using the representation in (5.1) one obtains These relations can be inverted to express the kernels as functions of the anomalous dimensions [15] 9) and P N +1 is the Legendre function. All singularities of the anomalous dimensions have to lie to the left of the integration contour. The algebraic structure of the three-loop anomalous dimension γ inv (N ) [6] is, unfortunately, too complicated to do the integrals in (5.8) analytically. We, therefore, adopted the following strategy: We will provide analytic expressions for the terms that correspond to 1. the leading contributions ∼ (ψ(j) − ψ(2)) k j(j − 1) , j = N + 2, to the large-N expansion of the anomalous dimensions, 2. the contribution of the leading singularity (pole at N = −1) in the complex plane.
The remainder will be parameterized by a sufficiently simple function with a few fit parameters.

Splitting functions
It turns out to be advantageous to use the representation for the anomalous dimensions in terms of the splitting functions. The three-loop splitting functions are available from Ref. [6] and involve harmonic polylogarithms (HPL) up to weight five, see Ref. [29]. Using this result and Eq. (5.3) it is straightforward to calculate the splitting functions for the invariant kernels such that The "plus" function can be written as We want to find a parametrization for the splitting functions consistent with the reciprocity relations (5.12) and separating the leading contributions at x → 0 and x → 1.
To this end, let us define the set of functions φ k (x) such that (5.13) They can be constructed recursively, so that where here and belowx We write the splitting functions as where the addenda, δH inv (x) and δH (3−) inv (x), do not include, by construction, terms ln k x, k ≥ 1 (for x → 0) andx ln kx , k ≥ 0 (for x → 1). The maximum powers of the logarithms are found by inspection of the known analytic expression. Thus, with normalization constants or, equivalently, in moment space,  Table 1 where, for completeness, we also repeat the expressions for χ cusp (3.14). In all cases we show the coefficients for the following color decomposition: cusp . The remaining terms δH inv (x) contributes at most 6% to the full splitting function in the whole range 0 < x < 1 so that for all practical purposes it can be approximated by a simple expression with a few parameters.
Due to the reciprocity property (5.12) the functions δH inv (x) can be parameterized in the form We choose the following ansatz where a ± and b ± are fit parameters and the normalization constants H ± 0 are determined analytically from the condition δH The fitted values of the parameters a ± and b ± for the different color structures can be found in Table 2. With this simple parametrization we reduce the deviation from the exact splitting functions to less than 0.5%, see Fig. 1.

Mellin transformation
The following Mellin representation of the kernels χ inv (τ ) and χ P inv (τ ) proves to be useful in order to restore them from the splitting functions and allows one to write all terms in the form that automatically respects the reciprocity relation: inv (ρ)Γ 2 (1 + ρ) , the asymptotic expansion of the integrals in Eqs. (5.24) at j N → ∞ can be obtained by moving the integration contour to the right and picking up the corresponding residues. It is easy to check that if the only singularities of the Mellin-transformed kernels χ(ρ) in the right half-plane are poles (of arbitrary order) at real integer values of ρ, then a generic term of the asymptotic expansion of the anomalous dimensions has the form which is required by the reciprocity symmetry under the j N → 1 − j N transformation [18][19][20]. Under the same condition ( χ(ρ) only has poles at integer ρ values in the right half-plane), the kernels χ(τ ) at small values of the conformal ratio have the expansion The Mellin transform of the one-loop kernel is very simple: More examples are collected in Table A and similar for χ inv (x). The kernels in τ space can finally be obtained by the inverse Mellin transformation (5.23a). We found this two-step approach to be the most effective for the three-loop case.
For the remainder function δH and for the simple ansatz in Eq. (5.22) Using Eq. (5.23a) we finally obtain the following expression for the corresponding contribution to the invariant kernel We also need the expressions for the functions φ k (5.14) in ρ-and τ -space defined as One obtains and where H p...q (z) are harmonic polylogarithms, see Ref. [29].
With these expressions at hand, our result for the invariant kernel (5.1) is complete. We obtain where the functions ϕ k (τ ) are defined in Eq. (5.33), the coefficients B can be found in Table 1 and the parameters for δχ

From light-ray to local operators
Light-ray operators are nothing but the generating functions for the renormalized local operators so that the mixing matrices for flavor-nonsinglet local operators can be calculated, in principle, by evaluating the evolution kernels on the test functions of the form f (z 1 , z 2 ) = z m 1 z k 2 , cf. (2.2). The results in this form are required for several applications, e.g. the calculation of moments of the distribution amplitudes and generalized nucleon parton distributions using lattice QCD techniques where the precision is increasing steadily and in some case already now requires NNLO accuracy [5].
Instead of using mixing matrices for the operators with a given number of left and right derivatives, as in Eq. (2.2), it proves to be more convenient to go over to the Gegenbauer polynomial basis. To the leading-order accuracy these operators diagonalize the evolution equations. Apart from convenience, writing the results in this form will allow us to make explicit connection to the formalism and notations used in [13] where the NLO expressions have been presented in this basis. We will see that the local operator formalism has its own advantages, e.g., solving the conformal constraint (2.21) is significantly easier in this language. Also the final step, reconstructing the invariant kernels from the eigenvalues, can be completely avoided here, as they directly enter the anomalous dimension matrices as diagonal elements.
Our goal is to translate the evolution kernels for light-ray operators into the anomalous dimension matrices for local operators of the form Here k is the total number of derivatives and the operator of the lowest dimension for given n, O n ≡ O nn , is a conformal operator (lowest weight of the representation of the SL(2) group). Increasing k for fixed n corresponds to adding total derivatives. The operators O nk mix under the evolution The mixing matrix γ nn ′ is triangular and its diagonal elements are equal to the anomalous dimensions Since γ nn ′ does not depend on k, the second subscript k is essentially redundant. In what follows we will use a "hat" for the anomalous dimensions and other quantities in matrix notation The light-ray operator (2.2) can be expanded in terms of the local operators defined in Eq. (6.1) where the coefficients Φ nk (z 1 , z 2 ) are homogeneous polynomials of two variables of degree k [30]: + ) k−n z n 12 , ω nk = 2 2n + 3 (k − n)! Γ(n + 2) Γ(n + k + 4) .
These polynomials are mutually orthogonal and form a complete set of functions 3 w.r.t. the canonical SL(2) scalar product (see, e.g., [30]) Φ nk |Φ n ′ k ′ = δ kk ′ δ nn ′ ||Φ nk || 2 = δ kk ′ δ nn ′ ω nk ρ −1 n , ρ n = 1 2 (n + 1)(n + 2)! (6.7) so that the local operators (6.1) can be obtained by the projection on the corresponding "state" The (canonical) conformal spin generators S (0) ± act as rising and lowering operators on this space whereas S (0 Thus the set of coefficient functions Φ n,k (z 1 , z 2 ) for k = {n . . . ∞} forms an irreducible representation of the SL(2) algebra, which is usually referred to as the conformal tower. Let A be a certain operator A acting on quantum fields. Its action can be realized by the expansion in terms of local operators with "matrix elements" serving as expansion coefficients Alternatively one can represent A by some integro-differential operator A acting on the arguments z 1 , z 2 of the light-ray operator O(z 1 , z 2 ) and, by means of the expansion (6.5), on the coefficient functions, Comparing the representations in Eqs. (6.10) and (6.11) we see that the action of A on the coefficient functions of local operators is given by the transposed matrix Using the orthogonality relation (6.7) one obtains which is the desired conversion, for a generic operator, between the light-ray and local operator representations.
The "matrix elements" A kk ′ nn ′ depend in general on four indices. However, if the operator A has a certain (canonical) dimension, i.e [S To translate this equation into the local operator representation, we define the matrices a mn (k) = m, k|S γ mn = m, k|H|n, k , w mn = m, k|z 12 ∆|n, k − 1 . (6.15) The first two matrix elements are easily computed, where we introduced a discrete step function ϑ mn = 1 if m − n > 0 and even 0 else.
The remaining two are nontrivial and can be written as a perturbative expansion Eq. (6.14) becomes in matrix notation Note that the matrices a(k) and b(k) (6.16) depend in principle on the total number of derivatives k. However, due to the fact that only diagonal elements depend on this parameter, the dependence on k drops out in the commutator. Hence we can safely omit it.
In complete analogy to the light-ray operator formulation, this equation fixes the non-diagonal (i.e. canonically non-invariant) part of the anomalous dimension matrix. Indeed, the commutator on the l.h.s. of Eq. (6.18) takes the form [ a, γ(a)] mn = (−a(m, k) + a(n, k))γ mn = −a(m, n)γ mn , (6.19) so that the non-diagonal elements of the mixing matrix are given by [8] γ In particular to the two-loop accuracy Here γ (1) is the well-known (diagonal) matrix of one-loop anomalous dimensions γ (1) mn = γ (1) n δ mn = 2δ mn C F 4S 1 (n + 1) − 2 (n + 1)(n + 2) − 3 (6.23) and w (1) is the one-loop conformal anomaly where Collecting everything one obtains the two-loop anomalous dimension matrix: The first few elements (0 ≤ n ≤ 7, 0 ≤ m ≤ 7) for N c = 3 are (6.27) Our expressions for the one-loop conformal anomaly and the two-loop anomalous dimension matrix coincide identically with the results in [13]. 4 Expanding Eq. (6.20) to the third order, we obtain the three-loop nondiagonal anomalous dimension matrix in the form In addition to the already known quantities, this expression involves the matrix element of the two-loop conformal anomaly (B.6), where ∆w (2) mn = m, k|z 12 ∆ (2) The explicit expression for the operator ∆ (2) + (B.7) can be found in [16]. We have not found a closed analytic expression for the matrix ∆w (2) mn , but the values for given m, n can be evaluated in a straightforward way.
Splitting the result into different color structures, we get for the first few elements (0 ≤ n ≤ 5, 1 ≤ m ≤ 7)  Using these expressions and the diagonal matrix elements from [6] we obtain the full three-loop anomalous dimension matrix 1 , . . .} + γ 1 + n f γ where the off-diagonal matrices for N c = 3 and different powers of n f in the range 0 ≤ n ≤ 7, 0 ≤ m ≤ 7 are given by the following expressions: (6.33) For completeness we also write the first few anomalous dimensions [6]: the anomalous dimension matrix for n f = 4 in the same range (0 ≤ n ≤ 7, 0 ≤ m ≤ 7):

Conclusions
Using the two-loop result for the conformal anomaly obtained in Ref. [16] we have completed here the calculation of the three-loop evolution kernel for the flavor-nonsinglet leading-twist operators in off-forward kinematics. The result is presented in the form of the evolution equation for the relevant nonlocal light-ray operator. In addition we derive the explicit expression for the three-loop anomalous dimension matrix for the local operators of dimension D ≤ 10, i.e., containing up to seven covariant derivatives. In the latter form, our result is directly applicable to the renormalization of meson distribution amplitudes and will be useful for lattice calculations of their first few moments. Practical methods for the solution of the three-loop evolution equations in more general GPD kinematics still have to be developed. Our results can, in principle, be translated to the evolution equation for GPDs by a Fourier transformation, although algebraic complexities of this transformation may prohibit obtaining the analytic expressions. An alternative method using the Mellin transformation in the conformal spin [31,32] is very attractive and it has become the standard tool in the NLO analysis of the DVCS [33][34][35] and deeply-virtual meson production [36,37]. The extension of this technique to NNLO was considered in [33][34][35] in a special "conformal" renormalization scheme. The transformation to the conventional MS scheme can be done using the results of our paper, but it remains to be seen whether this works in practice without major complications. The NNLO analysis of the DVCS data will, of course, require the extension of our results to flavor-singlet operators.

Acknowledgments
This study was supported by Deutsche Forschungsgemeinschaft (DFG) with the grants MO 1801/1-2 and SFB/TRR 55.

B X kernels
In this appendix we present explicit expressions for the kernels X (k) appearing in the operator of similarity transformation (3.1). The one-loop kernel X (1) is defined as a solution to the differential equation The operators X IA , X IB which originate from the two-loop anomaly ∆   It is easy to check that eigenvalues of the kernels defined in (C.2) and (C.4) are given by the following expressions: 1 z N 12 = T  This choice is convenient for our present purposes as it leads to a certain symmetry of the three-loop invariant kernel H inv that allows one to obtain somewhat simpler expressions, see Sec. 5.