Two-loop evolution equations for flavor-singlet light-ray operators

QCD in non-integer $d=4-2\epsilon$ space-time dimensions enjoys conformal invariance at the special fine-tuned value of the coupling. Counterterms for composite operators in minimal subtraction schemes do not depend on $\epsilon$ by construction, and therefore the renormalization group equations for composite operators in physical (integer) dimensions inherit conformal symmetry. This observation can be used to restore the complete evolution kernels that take into account mixing with the operators containing total derivatives from their eigenvalues (anomalous dimensions). Using this approach we calculate the two-loop (NLO) evolution kernels for the leading twist flavor-singlet operators in the position space (light-ray operator) representation. As the main result of phenomenological relevance, in this way we are able to confirm the evolution equations of flavor-singlet generalized hadron parton distributions derived earlier by Belitsky and M\"uller using a different approach.


Introduction
The exciting possibility to study the three-dimensional proton structure ("tomographic imaging") at the next generation lepton-hadron and lepton-nucleus colliders, most notably the EIC [1], poses formidable theoretical challenges. To meet these challenges precise predictions in quantum chromodynamics (QCD) are necessary and require pushing the corresponding computational tools to the highest possible accuracy. In practice, however, there exists still a considerable gap between precision and theoretical rigor applied in data analyses at the energy frontier, e.g., by the LHC community searching for New Physics beyond the Standard Model, and the QCD studies at medium energy machines, where the main goal is an improved understanding of the strong interactions.
This concerns, in particular, differences in the accuracy when considering the scale dependence of parton distributions. Whereas the next-to-next-to-leading order (NNLO), i.e., three-loop, analysis of parton distributions and fragmentation functions is becoming standard [2], the analysis of deeply-virtual Compton scattering and, e.g., hard exclusive vector meson electro-production is still often based on the leading-order (LO) evolution, despite the complete next-to-leading order (NLO) evolution kernels being available for a long time [3]. The necessity to close this gap is becoming gradually accepted so that the task to verify the results of ref. [3] by an independent calculation and to develop the techniques to push such calculations to NNLO accuracy are high on the agenda.

JHEP02(2019)191
On the technical level, the challenge is that considering off-forward matrix elements one has to take into account mixing with the operators containing total derivatives. The direct calculation of the relevant Feynman diagrams is difficult as they involve two different momenta, so that the method of choice has been to use the constraints from conformal symmetry that allow one to "save" one loop in comparison to the direct computation. This approach was pioneered by Dieter Müller [4][5][6] and further refined in refs. [3,7,8]. The existing results for the two-loop (NLO) evolution kernels [3] have been obtained using this technique.
In ref. [9] we proposed a somewhat different implementation of the same idea, based on the exact conformal symmetry of QCD in d = 4 − 2 dimensions at the critical point. Another difference to Müller's approach is our use of the position space (light-ray operator) representation that allows one to switch rather easily to momentum fraction space, avoiding the generally nontrivial problem of the restoration of the evolution kernels from the mixing matrices of local operators. Evolution kernels in position space are also interesting on their own in connection with lattice QCD calculations of Euclidean "observables" that can be factorized in terms of parton distributions, see e.g., refs. [10][11][12]. Such calculations currently attract a lot of attention.
This modified approach was tested in [9] on several examples to two-and three-loop accuracy for scalar theories. In [13][14][15] we used it to calculate the NLO (two-loop) [13] and NNLO (three-loop) [14,15] evolution kernels for the leading-twist quark-antiquark flavor-nonsinglet operators. In this work we address the flavor-singlet sector. This case is technically more complicated and also new questions arise concerning the role of gaugenoninvariant operators in conformal Ward identities. These issues will be discussed in what follows. Our main result is the derivation of the complete set of NLO evolution kernels in the position space representation for the leading twist gluon and (C-even) quark-antiquark operators. Expanding these kernels at small field separations we reproduce the results for the mixing matrices for flavor-singlet local operators given in ref. [7]. We also demonstrate how the evolution kernels in momentum fraction space can be obtained from the position space expressions by simple integration. We have verified numerically that the resulting expressions agree with the kernels that are implemented in the NLO evolution code for generalized parton distributions (GPD) by Freund and McDermott [16,17] which is based on the analytic expressions derived in ref. [3].

JHEP02(2019)191
Here n µ is an auxiliary light-like vector, n 2 = 0, the "plus" projection stands for F +µ = n ν F νµ , and z 1 , z 2 are real numbers. We will often omit the factor n µ in the argument specifying the field position and use a short-hand notation q(z) ≡ q(zn), etc. The gauge links in eq. (2.1) are taken in the adjoint and fundamental representations, respectively, [z 1 , z 2 ] = P exp igz 12 where A + = A a + T a and in the adjoint representation T a bb = if bab . In eq. (2.2) we have introduced another notation that will be used throughout this work: The LROs (2.1) can be viewed as the generating functions for local quark and gluon operators where χ = q, g and one can choose, for example, n (x) are Gegenbauer polynomials. For the C-parity-even operators considered here the sum goes, obviously, over odd n. The "coefficient functions" Ψ (χ) nk (z 1 , z 2 ) are homogeneous polynomials of two variables z 1 , z 2 of degree k for quarks and k − 1 for gluons. The factor 6 in the definition of the gluon operator is inserted for convenience. It ensures uniform normalization of the coefficient functions later on, see eq. (5.2).
Renormalized LROs are defined by the same expression with bare local operators replaced by the renormalized ones (we always assume dimensional regularization and minimal subtraction): (2.6) Let us stress that the polynomials Ψ nk in eq. (2.6) are exactly the same as in eq. (2.4).
Note also that the LRO on the l.h.s. of eq. (2.6) does not depend on the particular choice of the local operator basis. Going over to a different basis would yield different expressions for the polynomials Ψ nk such that the sum is unaffected. The renormalized local operators (2.5) satisfy the renormalization group equation (RGE) which has the usual form where M is the renormalization scale The anomalous dimension matrix γ χχ nn has a triangular form in n-space: its elements are nonzero only for n ≥ n . The diagonal entries γ χχ nn are 2 × 2 matrices in the χ, χ ∈ {q, g} space, they are known to three-loop accuracy [19]. The off-diagonal entries in n-space vanish to one-loop accuracy for the special choice of local operators in eq. (2.5). Beyond one loop the non-diagonal entries are non-zero and their calculation to two-loop accuracy is the main topic of this study.
The RGEs for local operators of different dimension can be combined to the RGE for their generating function, the LRO. In this representation one obtains an integrodifferential equation where H(a) (evolution kernel) is an integral operator which has a perturbative expansion The one-loop kernel H (1) χχ takes the form [18] Here we introduced the shorthand notations for the integrals The one-loop kernel commutes with the canonical generators of the (collinear) conformal transformations, where The expressions for the one-loop kernel in eq. (2.12) have been obtained in ref. [18] by direct calculation. But in fact, they are determined completely by the anomalous dimensions γ (1),χχ n , see ref. [20] and a discussion in ref. [15]. Starting from two loops this property is lost: the evolution kernel does not commute any longer with the canonical generators and cannot be restored from the anomalous dimensions alone. Nevertheless, one can simplify the calculation considerably by observing that QCD in noninteger d − 2 dimensions possesses a nontrivial critical point such that the β-function vanishes for the special value of the coupling, β(a * ) = 0, and the theory enjoys full conformal invariance. This property allows one to argue that eq. (2.14) holds true for the full kernels in arbitrary order of perturbation theory with the appropriately modified ("deformed") symmetry generators S α . A general technique for the calculation of quantum corrections to the generators of conformal transformations was developed in [9,13,14]. The present case (flavor-singlet) is more complicated as compared to the discussion of flavor-nonsinglet operators in ref. [14] because of a nontrivial mixing of gauge-invariant LROs with BRST [21] and Equation-of-Motion (EOM) operators. We consider this problem in more detail in the next section.

Light-ray operators beyond one loop
It is well known that gauge-invariant (local) operators mix under renormalization with the BRST variations and EOM operators which are not gauge-invariant. A renormalized operator [O χ,nk ] can be decomposed as [22][23][24][25] [O χ,nk (0)] = O χ,nk + B χ,nk + E χ,nk , where the three terms on the r.h.s. are the gauge-invariant, BRST and EOM operators, respectively. The last two terms do not contribute to the correlators of gauge invariant operators and in many cases can be omitted. The (renormalized) LROs can be decomposed in a similar manner and one can expect that in most cases only the gauge-invariant contributions will prove to be relevant, By definition, the variation of the LRO under collinear conformal transformations δ α , α = 0, ± (or, equivalently, dilatation and rotation in the (n,n) plane, D − M nn , translation along the light-cone P + and special conformal transformation K − , cf. appendix A) is determined by the transformation properties of local operators [9,13,14] 3)

JHEP02(2019)191
where δ α O χ,nk (0) can be expanded over the basis of gauge-invariant local operators with finite coefficients 2 δ α O χ,nk (0) = χ ,n (c α ) χχ nn O χ ,n ,k+α (0) + . . . , (3.4) where the ellipses stand for contributions of gauge non-invariant operators. For dilatations, the structure of such terms repeats eq. (3.1) -one obtains local operators that can be presented as a BRST variation, and EOM operators. For the special conformal transformations the structure of such gauge non-invariant contributions is more complicated, and they are in general non-local. A detailed discussion and explicit expressions can be found in ref. [26] where it is proven that such extra terms do not contribute to correlation functions with gauge-invariant operators. For the symmetry transformations from the Poincaré group there are no such contributions, i.e., they involve only the sum over gauge-invariant operators.
Note that the expressions for the symmetry transformations depend on the choice of the local operators but this dependence is compensated by the modification of the coefficient functions. As the result, the symmetry generators in the LRO representation are defined unambiguously and do not depend on the choice of the operator basis.
Omitting gauge non-invariant contributions and substituting eq. (3.4) in eq. (3.3) one can represent the result as an action on the LRO of a certain linear operator The generators S α,χχ , α = ±, 0 are integro-differential operators in z 1 , z 2 , which are 2 × 2 matrices in the quark-gluon space, χ, χ ∈ {q, g}. From general considerations [9,14] it follows that they can be written in the form where the canonical generators, S α , are defined in eq. (2.15). Thus quantum corrections to the classical symmetry generators involve the evolution kernel H and the operator ∆ which we will refer to as conformal anomaly. Both of them are matrices in the quark-gluon space.

Flavor-singlet conformal anomaly ∆
The leading-order (one-loop) conformal anomaly for the case under consideration

JHEP02(2019)191
has been calculated in ref. [7] from the analysis of the scale and conformal Ward identities. Our approach is in principle similar but differs from [7] in details. The calculation is described in appendix B. We obtain .

(3.8)
This expression is in agreement with ref. [7]. The general procedure how one can restore the evolution kernels from the spectrum of anomalous dimensions at the same order in perturbation theory and the conformal anomaly at one order less is explained in detail in [15] and it can be used for the flavor-singlet operators as well.
The central observation is that the evolution kernel at critical coupling must commute with the generators of conformal transformations, , (3.11) where h inv (τ ) are functions of τ = αβ αβ dubbed "the conformal ratio". Indeed, the one-loop evolution kernel eq. (2.12) can be rewritten in this form [13]. Eq. (3.10b), in turn, can be viewed as an inhomogeneous differential equation for the two-loop evolution kernel H (2) and is solved by [15] where H inv can be any invariant operator that takes the form eq. (3.11) and the operators X (1) and T (1) are 2 × 2-matrices which must satisfy the constraints

JHEP02(2019)191
Solving these two equations we get 14) The operator T (1) takes the form .
inv which still has to be determined. The expressions given above define the non-invariant part of the two-loop evolution kernels in the factorized form as a product of relatively simple integral operators (3.14), (3.15). For certain applications and in particular for numerical studies of the scale dependence of the GPDs, it can be advantageous to have explicit expressions for these products.
The results can be written in the form For the pure singlet-contribution to the quark kernel we get ∆r PS qq = ϑ PS qq = 0 and

JHEP02(2019)191
where the "plus" distribution is defined as Finally, for the gluon-gluon kernel we obtain The expressions in eqs. (3.19)-(3.24) supplemented by the invariant kernels which are calculated in the next section provide one with the complete two-loop flavor-singlet evolution kernels in the LRO representation. These are our main results.

Anomalous dimensions vs invariant kernels
Substituting the LRO in eq. (2.10) by its expansion in terms of local operators and comparing the resulting expression with the RGE (2.7) one gets Conformal operators of the lowest dimension for the given spin correspond to the highest weights of the representation and are annihilated by S − = S (0) − . As a consequence, the coefficient functions of the operators with k = n are translation-invariant and by dimension counting Ψ q nn ∼ z n 12 and Ψ g nn ∼ z n−1 12 . It follows that where the coefficients h χχ are related to the matrix of anomalous dimensions γ χχ (n) in the normalization chosen in ref. [19] as h χχ (n) = 2γ ref. [19] χχ (n + 1) , h qg (n) = (2/n)γ ref. [19] qg (n + 1) , h gq (n) = 2n γ ref. [19] gq (n + 1) .

JHEP02(2019)191
Eq. (3.10b) allows one to find the spectrum of the invariant kernel H (2) inv from the known results for the two-loop anomalous dimensions. Following the notation in ref. [19] we split the two-loop flavor-singlet quark kernel in two parts -flavor non-singlet and pure-singlet, PS . The expression for the non-singlet part can be found in ref. [15], while for the pure singlet one finds where h X χχ e χ n = X χχ (n)e χ n , .
and the two-loop expression for the pure-singlet anomalous dimension from ref. [19], we get (4.7) The corresponding kernel as a function of the conformal ratio takes the form In a similar manner one obtains from eq. (3.10b) the equations for the remaining entries: The relevant eigenvalues of the X kernels are equal to cf. eq. (4.5) for X (1) gq , and the eigenvalues for the diagonal entries of the one-loop evolution kernels take the familiar form Their derivatives can be written as Combining all these expressions we obtain, after some algebra, the following results for the eigenvalues of the invariant kernels

JHEP02(2019)191
and h inv gg (n) = 16 3 (4.16) It remains to restore the invariant kernels in the LRO representation from these results for the spectrum. The general expression for H inv can be written as a 2 × 2 matrix The entries are invariant operators which we parameterize as follows , where τ = αβ/ᾱβ. Adding the pure-singlet kernel (4.9) to the expression for the flavor non-singlet kernel given in refs. [13,15] one obtains For the quark-gluon kernel we find Γ qg = r qg = 0 and Finally, for the gluon-gluon kernel we obtain and ω gg (τ ) = 16 3 C A n f τ lnτ − 50 3 + 23 + 40C F n f τ lnτ + 2 5 + 2τ We close this discussion with a remark on the so-called reciprocity symmetry of the invariant kernels [27]. This symmetry arises, technically, from the observation [27] that invariant kernels can in general be presented in terms of the quadratic Casimir operator of the collinear conformal group. For an operator with conformal spin j (in our case j = n + 2) and its anomalous dimension given by γ(j) = f (j + 1 2 γ(j)), the asymptotic expansion of the function f (j) for large j consists of terms invariant under j → 1 − j. It was shown in ref. [15] that the function f the situation that there are two and or more operators of the same conformal spin, the reciprocity cannot be expected in general for the off-diagonal elements of the anomalous dimension matrix, because they depend explicitly on the assumed normalization for the operators [27]. Since the evolution kernels for the LROs do not depend on the basis of local operators, one should expect, however, that in this representation the reciprocity holds for the off-diagonal elements (kernels) as well. Indeed, one can verify that the eigenvalues of the off-diagonal invariant kernels, eqs. (4.14) and (4.15) are invariant under j → 1 − j with j = n + 2. For the invariant kernels themselves this condition implies an expansion k h k (τ ) ln k τ for τ → 0, where h k (τ ) is an analytic function in the vicinity of τ = 0.

Anomalous dimension matrix for local operators
For a comparison with ref. [7] and also for the application to the scale-dependence of the meson light-cone distribution amplitudes, it is desirable to have the results also in a different form, as an anomalous dimension matrix for local operators in the Gegenbauer basis (2.5).
These polynomials are mutually orthogonal and form a complete set of functions w.r.t. the canonical SL(2) scalar product (see, e.g., ref. [28]) The local operators (2.5) can be obtained by the projection of the LROs on the corresponding coefficient function In order to translate the action of the evolution kernel H in the LRO representation in eq. (3.12) in terms of local operators we adopt the following dictionary: where A ∈ {H, X, T}. Note that for all operators from this set [A, S − ] = 0 and, therefore, we can choose k = n in eq. (5.7) for convenience without loss of generality. Moreover, it is worth to be mentioned that for any invariant operator A χχ inv the matrix A χχ inv is diagonal: Making use of these properties we find the following results for the local matrix elements of the X operators (3.14) with the eigenvalues X where For the operator T (1) defined in eq. (3.15) we obtain with the one-loop anomalous dimensions in eqs. (4.6), (4.12), (4.13).
The result in this form can directly be compared to the two-loop anomalous dimension matrix calculated in ref. [7]. The diagonal elements H  ref. [7] g,nk , (5.15) and also that the anomalous dimensions in ref. [7] are written as an expansion in α s /(2π). Taking these differences into account we find a perfect agreement.

Evolution kernels in the momentum fraction representation
Evolution equations in the momentum fraction representation can be derived straightforwardly from the kernels in the coordinate representation. Let f (x 1 , x 2 ) be a Fourier transform of the position-space distribution f (z 1 , z 2 ) defined as The general expression for an evolution kernel in the momentum fraction representation takes the form (Hf )(x 1 , x 2 ) = DuH(x 1 , x 2 |u 1 , u 2 )f (u 1 , u 2 ) , (6.2)

JHEP02(2019)191
where Du = du 1 du 2 δ(x 1 + x 2 − u 1 − u 2 ) and the δ-function is due to the momentum conservation. Following the decomposition in eq. (4.18) we split the kernel H into three parts where H ϑ , H ω originate from the contributions which have the form in the coordinate space (LRO) representation , It is easy to show that the corresponding kernels in the momentum fraction representation are given by the following expressions: ( The Θ function is defined as follows Θ(a 1 , . . . , a n ) = and functions A ω , B ω , C ω are given by the following integrals where it is implied that x 1 + x 2 = u 1 + u 2 and The kernel ω(α, β) is a single-valued function in the simplex 0 < α + β < 1. It is easy to see that if the variables x i , u i belong to the regions determined by the corresponding Θfunctions, then the arguments α x , β x lie inside the simplex and the corresponding integrals are unambiguously defined. In the present case, ω(α, β) is symmetric under the interchange of its two arguments ω(α, β) = ω(β, α) (this is not always the case, see ref. [29]), and as a

JHEP02(2019)191
consequence C ω (x 1 , x 2 , u 1 , u 2 ) = B(x 2 , x 1 , u 2 , u 1 ). Provided that the functions B ω , C ω are continued in an appropriate way beyond their analyticity domain one gets A ω = B ω − C ω . Note also that all these kernels are effectively functions of two variables, x = x 1 /(x 1 + x 2 ) and y = u 1 /(x 1 + x 2 ), e.g. A ω = a(x, y)/(x 1 + x 2 ), etc. All integrals that arise in the transition from the position to the momentum fraction representation can be taken analytically in terms of elementary functions and Li 2 polylogarithms. The resulting expressions are collected in ref. [3] and are implemented in the NLO evolution FORTRAN code by Freund and McDermott [16,17]. 3 We have checked numerically that our results for the kernels in the momentum fraction representations agree with the kernels implemented in this code for all color structures. In this way we can confirm the results of ref. [3], where these kernels are given in analytic form. 4

Summary
We have presented a re-derivation of the two-loop flavor-singlet evolution equations for the leading-twist operators in off-forward kinematics, based on using conformal symmetry of QCD in non-integer d − 2 space-time dimensions at the critical point. This case is more complicated as compared to the flavor-nonsinglet evolution, both, technically and conceptually, due to potentially dangerous contributions of gauge-noninvariant operators. Our analysis is based on studies of conformal Ward identities in ref. [26], where it is proven that such extra terms do not contribute to correlation functions with gauge-invariant operators.
The results in this work are given in the position-space or light-ray operator representation. This form has some technical advantages and can also be interesting in applications to lattice QCD calculations. Expanding our results in powers of the field separation we reproduce the results for the mixing matrices for flavor-singlet local operators derived in ref. [7]. The evolution kernels in momentum fraction space can be obtained from our expressions by simple integration. We have verified numerically that the resulting kernels agree with the kernels that are implemented in the NLO GPD evolution code by Freund and McDermott [16,17] which is based on the analytic expressions derived in ref. [3].
To summarize, the two-loop evolution equations for generalized parton distributions have now been derived independently by two groups, and perfect agreement is found. In view of the projected high statistics of the relevant experiments at the JLAB 12 GeV upgrade and, in future, the EIC, it is necessary to implement and use these results now in NLO analyses of deeply-virtual Compton scattering and similar reactions.

JHEP02(2019)191
A Scale and conformal transformations Scale (D) and conformal (K) field transformations for the fundamental fields take the generic form where ∆ Φ = dim Φ are the canonical dimensions of the fields. It is convenient to choose them in 4 − 2 dimensions to be exactly the same as in the four-dimensional theory, For this choice the gluon strength tensor F σρ transforms in a covariant way and the covariant derivative of the ghost field D ν c transform as a vector field, The variation of the different parts of the QCD action under the special conformal transformation takes the form Note that the ghost and the gauge-fixing terms break conformal symmetry explicitly even in d = 4 dimensions. Summing up all contributions one obtains

B One-loop conformal anomaly
Our analysis follows closely the lines of ref. [14] so that in this appendix we concentrate mainly on specific problems that arise for the flavor-singlet operators. Consider [14] the correlation function of two (renormalized) LROs stretched in different directions where z = {z 1 , z 2 }, w = {w 1 , w 2 } and n,n are two auxiliary light-like vectors, n 2 =n 2 = 0.
In what follows we assume that (nn) = 1 and (x · n) = (x ·n) = 0. Starting from the representation in terms of local conformal operators it can be shown, see ref. [14], that the correlation function (B.1) satisfies the constraint where the operator S +,σχ is the generator of special conformal transformations defined in eq. (3.5). The explicit expression for S +,σχ can be found from the analysis of the conformal Ward identity (CWI) for the correlation function (B.1). Making the corresponding change of variables in the path integral representation for this correlation function one obtains where δ stands for the conformal variation along then direction, δ =n µ δ K µ . The second term on the l.h.s. is solely responsible for the r.h.s. of eq. (B.2) while the two others contribute to the l.h.s. of this relation. The first of them, the variation of the LRO operator, takes the form The renormalization factor Z in this equation is an integral operator in z 1 , z 2 and it does not commute with S Since all singular terms must cancel in the final result we can drop them. The second term that contributes to the l.h.s. of eq. (B.2) is due to the conformal variation of action as given in eq. (A.8). As the first step, it has to be re-expanded in JHEP02(2019)191 terms of renormalized operators. In Landau gauge that we will assume from now on, the corresponding expression takes the form , Ω Φ = ΦδS R /δΦ, Ω qq = Ω q +Ωq. Note that the anti-ghost EOM term does not contribute to the correlation function in question since the LROs O χ do not contain (anti)ghost fields. The terms Ω A and ∂ µ [B µ ] contain contributions of the gauge parameter proportional to 1/ξ, which, however, cancel each other. For nonzero ξ the correlation function of B µ with any gauge invariant functional F vanishes, implying that Integrating by parts the contributions of the EOM terms in δS R O (n) χ (x, w) we can write, e.g., for the first operator insertion, . Let us consider first the quark operator, O q (0, z). In order to streamline the notation, hereafter, we do not display the dependence on the auxiliary vector n and the space-time coordinates, if this is clear from the context. Our goal is to calculate the one-loop correction to the generators. It can be shown that there are no BRST or EOM counterterms to the quark LRO at one loop, i.e. O q (z) = Z qq O q (z) + Z qg O g (z) + O(a 2 ). It is also obvious that the one-loop correction to the quark part of the generator, S qq , is the same as in the non-singlet case. We consider therefore only the off-diagonal entry S qg , i.e. we are looking for the contributions of the gluon LRO O g on the r.h.s. of eq. (B.8).
First, note that the term 2γ c D A (x) O q (z) gives rise to a O(a 2 ) contribution only and can be omitted. Second, there are no gluon pair counterterms for the product L Y M (x) O q (z) at one loop. Thus, the only relevant contribution comes from the quark EOM term. We write The gluon contribution of interest corresponds to the last term in this expression. Taking into account that the renormalization factor is related to the evolution kernel as JHEP02(2019)191 qg + O(a 2 ) and multiplying this relation by − γ q , one can bring this contribution to the form This result implies that ∆ (1) qg = 0 in agreement with ref. [7]. Next, consider the S gq generator. To this end we again omit the D A term in eq. (B.8) since it is O(a 2 ). The quark EOM contribution in eq. (B.8) can be handled as above yielding In addition we have to consider the contribution from the product L Y M (x) O g (z). This product is divergent when x → 0 so that one has to add and subtract the corresponding pair counterterms. The counterterm ∼ O q corresponds to the Feynman diagram shown in figure 1 where the crossed circle stands for the insertion of the vertex The terms involving (∂A) vanish since the propagators are transverse in Landau gauge. The gluon line with the crossed-circle insertion (modified propagator) can be written in the form where G µν is the usual gluon propagator in Landau gauge and (B.14) Taking into account the contribution of the crossing-symmetric diagram one obtains (divergent part only) 5

JHEP02(2019)191
where {. . . , . . .} + stands for the anti-commutator. Multiplying this expression by (−β(a)/2a) and adding eq. (B.11) we obtain the result for ∆ (1) gq given in eq. (3.8). The calculation of the pure gluon contribution, ∆ gg , is more involved. It was shown in [26] that the r.h.s. of eq. (B.8) can be represented as sum of the fully renormalized product [F 2 µν (x)O g (z)] which enters with the coefficient β(a)/a, contributions of gauge-invariant operators and gauge-noninvariant contributions that are BRST variations of nonlocal operators or nonlocal EOM operators which drop out of all correlation functions with gaugeinvariant operators, hence also of eq. (B.3). Nevertheless, we need to know the structure of the gauge-noninvariant terms in order to correctly separate the gauge-invariant contributions. Consider first the BRST and EOM counterterms to the gluon LRO O g . Since they should be twist-two operators, there are two possible structures: Here we tacitly assume that the color indices are contracted in some way. Both operators contain singular terms in ξ that have to cancel in the sum. Therefore they can appear only in the combination, B k (z 1 , . . . , z k ) − E k (z 1 , . . . , z k ). The lowest, k = 2, counterterm can be restored by calculating the one-loop F +µ F +µ →cc diagram. The result reads (B.17) The BRST variation of this operator is given by a sum of EOM contributions which, in agreement with the general statement, do not contain the anti-ghost EOM. The contributions for k = 3, 4, . . . can be neglected to our accuracy as they are higher order in the coupling. Next we notice that for the calculation of ∆ gg to one-loop accuracy one can omit the quark EOM term ∼ D qq in eq. (B.8), and the remaining contribution ∼ D A (x) can be written as where (B − E) 2 is given in eq. (B.17) and the ellipses stand for k > 2 (gauge-noninvariant) counterterms which contribute only starting from O(a 2 ). The term with the BRST variation B 2 also drops out from the correlation function at order O(a). The remaining EOM term can be rewritten (inside the correlation function) with the required accuracy as .

(B.22)
This is the only contribution from the BRST and EOM counterterms at this order. Note that the operator in eq. (B.22) is not a BRST variation so that, according to the general statement of ref. [26], this contribution must be complemented by contributions of other diagrams to produce a BRST variation or EOM operator in the sum, which is not easy to see on a diagrammatic level. For our purposes it is, however, sufficient to note that all such "unwanted" contributions contain either A + or A − fields, and can be avoided if one looks for the contributions of transverse gluons.
In order to calculate the pair counterterm to the operator product L Y M (x) [O g (z)] it is convenient to make some rearrangements. Namely, we rewrite The integral of the first piece can be written as (up to (∂A) contributions) where L Y M 3,4 stand for the three-and four-gluon vertices. The ghost terms can be omitted as they do not contribute to the gauge-invariant counterterms. In the second piece, the BRST operator does not contribute to the correlation function and the Ω A -term gives All these contributions are multiplied by (−β(a)/2a) = 2( − γ g ).
The main contribution to the CWI is due to the pair counterterms to The corresponding diagrams are shown in figure 2 (adding the symmetric contributions is implied). They can be related to the counterterm diagrams for the gluon LRO with some extra terms.
Since the four-gluon vertex does not contain derivatives, the counterterm corresponding to the last diagram in figure 2 can be rewritten as where G ab µν (y, z 1 , z 2 ) is the usual QCD diagram with the four-gluon vertex and G ab µν (y, z 1 , z 2 ) is the diagram with the modified vertex, G ab µν (y, z 1 , z 2 ) ≡ 2(ny)G ab µν (y, z 1 , z 2 ). For the rest of the contributions in figure 2 that involve a three-gluon vertex (which contains a derivative), this relation has to be modified as follows. Let y 1 be the coordinate of the modified vertex (insertion of (n · y 1 )). Then where we suppressed color and Lorentz indices on the gauge field. In this expression, as above, G(y 1 , y 2 , z 1 , z 2 ) is the diagram with the standard QCD three-gluon vertex, G(y 1 , y 2 , z 1 , z 2 ) = (ny 1 )G(y 1 , y 2 , z 1 , z 2 ) and G 1 (y 1 , y 2 , , z 1 , z 2 ) is the Feynman diagram obtained from G(y 1 , y 2 , z 1 , z 2 ) by the replacement of the triple vertex V (q, r, k) at the position y 1 by the new vertex V (q, r, k) defined as V (q, r, k) = −i(n, ∂ q )V (q, r, k) q=0 . (B.29) Using these relations the divergent parts of the contributions in figure 2 can be presented in the form Here CT(O g ) stands for the counterterms to the operator O g , D i are the first four diagrams in figure 2 (with the triple-gluon vertex) with the crossed vertex replaced by eq. (B.29) and D i are the diagrams shown in figure 3 (upper row), where the big black circle stands for the vertices generated by the operator (nn)A µ (z 2 ) −n µ A + (z 2 ) + (nn)z 2 F +µ (z 2 ) . The first contribution in eq. (B.30) absorbs all terms with the derivative D A in eqs. (B.28) and (B.27) but also generates additional contributions, the diagrams D i shown in the upper row in figure 3, that have, therefore, to be subtracted. We recall that eq. (B.30) has to be multiplied by factor (−β(a)/2a) = − γ g and also the contribution due to eq. (B.25) has to be added, which reads  This expression has still to be rewritten as the sum of renormalized operators minus the corresponding counterterms. Since γ g (a * ) = the contributions of renormalized operators vanish at the critical point. For the first term we can write  but their sum can be simplified to the diagrams shown in lowest row. Strictly speaking, there is also a contribution of the form (z 1 − z 2 )O g (z 1 , z 2 ) that cancels in the sum with the symmetric diagram (z 1 ↔ z 2 ). Finally, the counterterm diagrams for the operator O 2g are shown in figure 4. The third and the fourth diagram in figure 4 cancel with the first and the second diagrams in the lowest row in figure 3, respectively. The remaining ones give rise to the conformal anomaly ∆ gg . These are: • The first four diagrams in figure 2 with the triple gluon vertex (B.29).
• The last diagram in the lowest row in figure 3.
• The first, second, fifth and sixth diagram in figure 4.
Technically, it is convenient to combine the first two diagrams in figure 2 with the first two diagrams in figure 4 as there are strong cancellations. The calculation of all other diagrams is rather straightforward. Note, that the final answer contains not only the term F µ+ F µ+ which we are interested in, but also gauge non-invariant contributions. These terms contain the field A + or A − or both of them and must be part of the BRST and EOM operators which can appear on the r.h.s. of eq. (B.8). We did not keep track of such contributions. But we have checked that the gauge-noninvariant terms, (nn)A µ (z 1 n)F +µ (z 2 n), which are neither BRST nor EOM operators, cancel out in the sum of all diagrams. The resulting expression for the conformal anomaly is presented in eq. (3.8). It has a very simply form and coincides with the result of Belitsky and Müller [7], that is the χχ K ω kernels in their notation.

C RG identities in Landau gauge
The analysis of the Ward identities in this work is based on the re-expansion of the renormalized Lagrangian 2 L R in terms of renormalized operators. The corresponding general expression in arbitrary covariant gauge reads [7,14,30] 2 Here L R = L R − 1 2 ∂ ρ (q 0 γ µ q 0 ), γ Φ are the anomalous dimensions of the fields, Ω Φ (x) = Φ(x)δS R /δΦ(x) and Ω µ =cD µ c − ∂ µc c is a conserved current, ∂ µ [Ω µ ] = Ω c − Ωc. The coefficients z b (g, ξ) and z c (g, ξ) are at this stage unknown functions.

JHEP02(2019)191
The gauge fixing term, [L gf ] = 1 2ξ (∂A) 2 , can be written as a combination of BRST and EOM operators, An arbitrary diagram with an insertion of this operator vanishes in Landau gauge. Hence, this contribution can safely be omitted. Next, whereas all counterterms are polynomials in the gauge parameter ξ, the terms Ω A and ∂ µ [B µ ] on the r.h.s. of eq. (C.1) contain contributions ∼ 1/ξ which have to cancel in the sum. This implies that z g (g, ξ) = γ A + γ g + O(ξ). The last coefficient z c (g, ξ) can be fixed in Landau gauge by observing that in this gauge the QCD action has an additional discrete symmetry corresponding to the interchange of ghost and antighost fields, c → −c,c → c. Indeed it is easy to check that Since the last term can be omitted in Landau gauge, the action is invariant with respect to this transformation. Clearly eq. (C.1) should respect this property. Taking into account that Ω c + Ωc is invariant under this transformation, Ω µ → −Ω µ , and one obtains z c (g, ξ → 0) = γ c = −(γ A + γ g )/2. Using these expressions we obtain the following result which holds in Landau gauge. Note, that here we do not keep the gauge fixing term L gf , since it does not contribute in this gauge.

D Evolution kernels vs. anomalous dimensions
In this appendix we collect the expressions for the invariant kernels that appear as building blocks in two-loop evolution equations and the anomalous dimensions. The anomalous dimension γ(j) corresponding to the kernel h(τ ) is defined as follows: The inverse relation reads [13] h(τ ) = 1 2πi C dj (2j + 1)γ(j) P j 1 + τ 1 − τ , (D.2)

JHEP02(2019)191
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.