NNLO anomalous dimension matrix for twist-two flavor-singlet operators

Conformal symmetry of QCD is restored at the Wilson-Fisher critical point in noninteger $4-2\epsilon$ space-time dimensions. Correlation functions of multiplicatively renormalizable operators with different anomalous dimensions at the critical point vanish identically. We show that this property allows one to calculate off-diagonal parts of the anomalous dimension matrices for leading-twist operators from a set of two-point correlation functions of gauge-invariant operators which can be evaluated using standard computer algebra techniques. As an illustration, we present the results for the NNLO anomalous dimension matrix for flavor-singlet QCD operators for spin $N\le 8$.


Introduction
The Electron-Ion Collider [1,2] will allow one to access generalized parton distributions (GPDs) [3,4,5] in a broad kinematic range.In particular the possibility to study the threedimensional gluon distributions in the longitudinal and transverse plane is new and very exciting.The scale dependence of GPDs is governed by evolution equations that are more complicated as compared to the usual parton distributions (PDFs).In the language of the operator product expansion (OPE), the added complication in this case is to take into account mixing with operators containing total derivatives.Going over to the momentum fraction space, this mixing translates to the evolution kernels involving extra variables.The complete set of the NLO (two-loop) evolution kernels is available for a long time [6] and the NNLO (three-loop) evolution kernels for flavornonsinglet operators were calculated more recently in [7].Both calculations use conformal symmetry constraints that allow one to obtain the kernels for GPDs from the known NLO and NNLO evolution kernels for PDFs and a computation of the so-called conformal anomaly from conformal Ward identities at one order less, i.e. a two-loop anomaly [8] is sufficient to obtain the NNLO kernels.The NNLO flavor-singlet kernels can, in principle, be obtained in the same way, but the calculation becomes too large to be done without using computer algebra methods.The required algorithmic implementation is, unfortunately, not available.
In this letter we suggest an alternative approach that allows one to calculate off-diagonal parts of the anomalous dimension (AD) matrices of local flavor-singlet operators from a set of two-point correlation functions which can be evaluated using standard computer algebra software packages 1 .The main advantage of this technique as compared to the direct calculation is that gauge non-invariant Equation of Motion (EOM) and BRST operators can be completely neglected.A disadvantage as compared to the approach of [6,7] is that the calculation is done for local operators with given (not very high) spin, alias for the first few moments of GPDs.The results can be used to obtain a certain approximation for the NNLO evolution kernels, but their construction is likely to be more complicated as compared to the well-studied case of PDFs.This is a separate problem that will not be considered here.
The starting point is that conformal symmetry of QCD at quantum level is restored at the Wilson-Fisher critical point [10] at noninteger space-time dimension where β 0 , β 1 ,. . .are the first few coefficients of the QCD βfunction and α s is the strong coupling.At the critical point, the two-point correlation functions of multiplicatively renormalizable operators with different ADs vanish to all orders of perturbation theory [12] [ where . . .stands for the vacuum expectation value.We will show that this condition allows one to find the eigenvectors of the renormalization group (RG) equation in the chosen operator basis from a calculation of the corresponding unrenormalized correlation functions with m ≤ n.Since the eigenvalues (ADs) are known, this information is sufficient to restore the complete mixing matrix.Last but not least, the ADs of composite operators in minimal subtraction schemes do not depend on by construction and are the same for the physical d = 4 and the critical d = 4 − 2 * space-time dimensions.Thus the calculated mixing matrix for the leading-twist operators at the critical point coincides identically with that in physical theory in four dimensions [13,8,7].
In this letter we will first explain application of this technique on a simple example in NLO, followed by a more systematic presentation for the most interesting case of flavor-singlet operators.The NNLO mixing matrix (in the Gegenbauer basis) for flavor-singlet QCD operators for spin N ≤ 8 presents our main result.As a byproduct of this calculation we re-derive and confirm the corresponding results of Ref. [7] for the flavornonsinglet operators.

Simple example
As an example, consider the twist-two operator where q 1 and q 2 are quark fields of different flavor, ∂ µ = ∂/∂x µ , 2 (y) is the Gegenbauer polynomial and ← D+, D+ are left and right covariant derivatives, respectively.The "plus" projection corresponds to a multiplication by an arbitrary light-like vector In processes involving a momentum transfer between the initial and the final states one needs to take into account mixing of O 2 (x) with the (second) total derivative of the vector current so that the renormalized operators in the MS scheme take the form It is convenient to introduce matrix notation Renormalized operators satisfy the RG equation Here is the AD-matrix and β(a) is the d-dimensional beta function where Since the vector current is conserved γ 11 = 0 and Z 11 = 1 to all orders in perturbation theory.The γ 22 entry is the usual AD of the leading-twist operator with two derivatives.It is known to five-loop order [14].For N c = 3 with etc.The advantage of using the Gegenbauer polynomial in ( 3) is that the off-diagonal ADs start at order O(a 2 ) in this basis: In what follows we describe a simple method to calculate γ (2)  21 .The mixing matrix (8) can be written in the following form with and set the space-time dimension to its critical value (1) such that the β-function (9) vanishes.With this choice, the RG equation in (7) decouples into separate equations for the "rotated" operators and conformal symmetry requires that to all orders of perturbation theory Using (5) we can rewrite this equation in terms of bare correlation functions This can be solved for A 21 or, equivalently, γ 21 , if the other entries are calculated to the sufficient accuracy.Let us note that Eq. ( 17) implies that the correlation functions This property is a consequence of conformal symmetry and is valid at the critical point only, → * .The renormalization factors in Eq. ( 18) take the form and, since γ 11 = 0, γ (2) 21 = A (1) 21 γ (1)  22 .Thus in order to find γ (2)  21 we need to calculate ), the second term on the l.h.s. of ( 18) can be omitted.The relevant Feynman diagrams are shown in Fig. 1.One obtains O where Using these expressions and the one-loop result for Z 22 (19), expanding everything to O(a) accuracy and replacing → −β 0 a one obtains from Eq. ( 18) in agreement with the known result [6,7].This calculation is much easier as compared to a direct calculation of γ (2)  21 from the two-loop Green function of O 2 and two quark fields.

General case
The approach sketched above can be generalized to all orders in perturbation theory and also for flavor-singlet operators.Let These operators have spin N = n + 1 and mix with each other under renormalization, Here and below [. ..] stands for a renormalization in MS scheme.Since operators containing total derivatives do not contribute to the forward matrix elements, these matrix elements satisfy the RGE of the form where α, β ∈ {q, g}.The anomalous dimensions are 2 × 2 matrices They are known to three-loop accuracy for all n [15] and to four loops for n = 1, 3, 5, 7 [16].In a theory in d = 4 − 2 dimensions the RGE (25) has the same form as in d = 4, but with the d-dimensional β-function (9).
In processes involving matrix elements with nonzero momentum transfer the RGE becomes more complicated.In this case mixing with operators containing total derivatives, has to be taken into account.For definiteness, and having in mind applications to two-photon reactions such as DVCS, we will consider C-parity-even operators n = 1, 3, 5, . . .(even spin).
which has the same form for all n, so that this subscript is essentially redundant.
It is convenient to introduce matrix notation and where each entry Z mk is a 2 × 2 matrix Z αβ mk .Note that the matrix Z m for m < n is a principal submatrix of Z n : The subscript only specifies the size of the matrix while the entries do not depend on it.The RG equation for where The diagonal entries γ nn are nothing else as the forward ADs (26), γ nn ≡ γ n , and our task is to find the off-diagonal entries γ km , k > m.In the chosen (Gegenbauer polynomial) operator basis the off-diagonal entries are O(a 2 ): The AD matrix (33) can be brought to the block-diagonal form where γ n (a) = diag{γ 1 (a), . . ., γ n (a)} and where all entries are 2 × 2 matrices in quark-gluon space, cf.(27).Next, define a rotated operator This equation means that the operators O α mm , α = q, g at = * can be written as linear combinations of two operators with certain scaling dimensions, which transform in a proper way under conformal transformations (dubbed conformal operators).Since correlation functions of conformal operators with different scaling dimensions vanish, we conclude that for m/ =n and x 0. For definiteness we assume n > m.
It proves to be convenient to write the operators O α nn in a slightly different form, The matrices A and B are related to each other as Then it follows from Eq. ( 39) Note that only one term with k = m survives in the sum on the r.h.s.We will show that this equation allows one to determine the coefficients B αγ nk .In practice, it is more convenient to do calculations in momentum representation.We consider the correlation functions of bare operators where a b is the bare coupling and s = µ 2 /(−p 2 − i0).A perturbative expansion for T αβ km can be written as where is the number of loops.The renormalized correlation functions [T αβ km ](s, a, ) are given by or, in matrix notation, [T](s, a, ) = Z T(s, Z a a, ) Z T .These functions still have a 1/ pole coming from the integration around x = 0 (recall that Eq. ( 42) holds only for x 0).This divergent contribution can be removed applying the derivative in s: This object is finite and we can put the space-time dimension to its critical value → * .In what follows we use a shorthand notation T * (s, a) = T(s, a, * ).The momentum-space version of Eq. ( 42) takes the form so that where Note that B and V (for n > m) are lower block-triangular and R is a block-diagonal matrix.The matrices V and R depend on B through A = (1l − B) −1 and implicitly through off-diagonal elements in the renormalization factors Z.It remains to expand Eq. ( 48) in powers of the coupling constant.Note that V = O(a) since the correlation functions O α n (x)O β m (0) with n m vanish in d = 4 at leading (one-loop) order.As a consequence, terms of order a k in the expansion of Eq. ( 48) only contain the B-matrix dependent terms of one order less on the r.h.s., so that it can be solved iteratively, orderby-order.Write Then where V 1 and R 0 are obtained from two-loop correlation functions (D 2 ) αβ km (44) and do not depend on B. Once B 1 is found, one can calculate the two-loop AD matrix and the two-loop renormalization factor As the next step, we obtain V 2 and R 1 with the input from threeloop correlation functions (D 3 ) αβ km .This allows one to calculate B 2 as and determine the three-loop AD matrix This procedure can be continued iteratively to any order, O(a k ), provided the correlation functions (43), (44) are calculated to the = k loops accuracy.
Finally, note that one can consider correlation functions of the operators (23) defined with two different auxiliary lightcone vectors n and n, schematically O (n) (x)O (n) (0).We have checked that this freedom does not produce new constraints while choosing n n complicates the calculations.

Flavor-singlet operators with spin N ≤ 8
We have calculated the correlation functions (D ) qq km , (D ) qg km , (D ) gq km , (D ) gg km as defined in Eqs. ( 43), (44) for k, m = 1, 3, 5, 7 to three-loop accuracy in d − 2 dimension for a generic gauge group.All the diagrams were generated with the help of QGRAF [17] and evaluated with FORM [18] programs MINCER [19] and COLOR [20].
The results are collected in the ancillary file.Using these expressions we determined the off-diagonal parts of the AD (mixing) matrices for C-parity even flavor-singlet operators.The complete expressions with all color structures are lengthy and are given in the second ancillary file.Here we present the results for N c = 3, separating contributions with the different n f dependence γ (2) = γ (2,0) + n f γ (2,1) , γ (3) = γ (3,0) + n f γ (3,1) + n 2 f γ (3,2) .
For the two-loop ADs we obtain These expressions coincide with those obtained in [6,21].The three-loop mixing matrix presents our main result:  As a byproduct of this calculation we have considered flavornonsinglet operators as well, and confirm the corresponding results of Ref. [7].
The size of the three-loop corrections for a ∼ 1/40 and n f = 4 is typically of the order of 20% of the two-loop results, with a few exceptions.The γ gq and γ gg entries are in all cases much larger than γ qg and γ qq .

Conclusions
We have presented a method to calculate off-diagonal parts of the mixing matrices of leading-twist operators with the operators including total derivatives based on conformal symmetry of QCD at the Wilson-Fisher critical point in noninteger dimensions.In this approach, the calculation of the ADs to -loop accuracy, O(a ), is reduced to a calculation of -loop gaugeinvariant correlation functions of leading-twist operators.As an illustration, we have calculated three-loop ADs of flavor-singlet operators for spin N ≤ 8 which contribute, e.g., to the moments of generalized parton distributions.The main advantage of this technique is that mixing with non-gauge-invariant operators can be ignored altogether and also the number of Feynman diagrams is much smaller as compared to the standard approach.An extension to higher moments and to four loops is straightforward but will require significant computer resources.Restoration of the off-forward evolution kernels in momentum fraction space from the results for a given set of moments is a nontrivial problem which goes beyond the task of this letter.
) and set the space-time dimension to the critical value → * , β(a)| = * = 0.The RGE for O n (at the critical point) takes the form µ∂ µ + γ n O n = 0 and decouples in n independent equations, µ∂ µ + γ m O mn = 0. Further, since O mn = ∂ n−m + O mm , the dependence on n is trivial and it is sufficient to consider the case n = m: