Two-loop conformal generators for leading-twist operators in QCD

QCD evolution equations in minimal subtraction schemes have a hidden symmetry: One can construct three operators that commute with the evolution kernel and form an $SL(2)$ algebra, i.e. they satisfy (exactly) the $SL(2)$ commutation relations. In this paper we find explicit expressions for these operators to two-loop accuracy going over to QCD in non-integer $d=4-2\epsilon$ space-time dimensions at the intermediate stage. In this way conformal symmetry of QCD is restored on quantum level at the specially chosen (critical) value of the coupling, and at the same time the theory is regularized allowing one to use the standard renormalization procedure for the relevant Feynman diagrams. Quantum corrections to conformal generators in $d=4-2\epsilon$ effectively correspond to the conformal symmetry breaking in the physical theory in four dimensions and the $SL(2)$ commutation relations lead to nontrivial constraints on the renormalization group equations for composite operators. This approach is valid to all orders in perturbation theory and the result includes automatically all terms that can be identified as due to a nonvanishing QCD $\beta$-function (in the physical theory in four dimensions). Our result can be used to derive three-loop evolution equations for flavor-nonsinglet quark-antiquark operators including mixing with the operators containing total derivatives. These equations govern, e.g., the scale dependence of generalized hadron parton distributions and light-cone meson distribution amplitudes.


Introduction
Scale dependence of physical observables in strong interactions involving a large momentum transfer is governed by the renormalization group (RG) equations for the corresponding (composite) operators. They have to be calculated to a sufficiently high order in perturbation theory in order to make the theory description fully quantitative. The anomalous dimensions of the leading twist-two operators are known to NNLO accuracy (three loops), and these results have been converted to the state-of-the-art NNLO evolution equations [1,2] for parton distributions that are used in modern description of inclusive reactions, e.g., at the LHC.
A remarkable progress in accelerator and detector technologies in the last decades has made possible the study of hard exclusive reactions with identified particles in the final state. Such studies have become a prominent part of the research program at all major existing and planned accelerator facilities. The relevant nonperturbative input in such processes 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 different momenta in the initial and the final state complicates the RG equations since mixing with the operators involving total derivatives has to be taken into account. The arising mixing matrix (for a given moment, or operator dimension) is triangular so that the diagonal entries correspond to the anomalous dimensions that are known to NNLO accuracy, but the nondiagonal contributions require a dedicated calculation.
A direct calculation in higher orders is quite challenging, however, it has been known for some time [3] that conformal symmetry of the QCD Lagrangian allows one to restore nondiagonal entries in the mixing matrix and, hence, full evolution kernels at given order of perturbation theory from the calculation of the special conformal anomaly at one order less. This result was used to calculate the complete two-loop mixing matrix for twist-two operators in QCD [4][5][6], and derive the two-loop evolution kernels for the GPDs [7][8][9].
In ref. [10] we have suggested an alternative technique, the difference being that instead of studying conformal symmetry breaking in the physical theory [4][5][6] we make use of the exact conformal symmetry of a modified theory -QCD in d = 4−2ǫ dimensions at critical coupling. Exact conformal symmetry allows one to use algebraic group-theory methods to resolve the constraints on the operator mixing and also suggests the optimal representation for the results in terms of light-ray operators. In this way a delicate procedure of the restoration of the evolution kernels as functions of two variables, e.g. momentum fractions, from the results for local operators can be avoided.
Utility of this modified approach was illustrated in [10] on several examples to twoand three-loop accuracy for scalar theories, and in [11] on the example of the two-loop evolution equation for flavor-nonsinglet operators in QCD. The present work is the first step towards the three-loop calculation in QCD. Our main result is the calculation of the two-loop contribution to the generator of special conformal transformations for flavornonsinglet leading-twist operators in QCD in non-integer d = 4 − 2ǫ space-time dimensions at critical coupling.
The presentation is organized as follows. Section 2 is introductory. We explain there the general strategy of our approach and introduce the necessary formalism and notations. Section 3 and related appendices A, B contain a detailed analysis of the scale and special conformal Ward Identities (WI) in d-dimensional QCD. The expression in eq. (3.47) for the ℓ-loop quantum correction to the generator of special conformal transformations is the main outcome of this analysis. In section 4 we explain some technical issues that one encounters in the calculation. The results for separate Feynman diagrams are collected in appendix C. Section 5 contains our principal result: the two-loop expression for the generator of special conformal transformations. The two-loop expression for the evolution kernel in the light-ray operator representation [11] is given as well. The final section 6 contains a short summary and outlook.

QCD in 4 − 2ǫ dimensions at the critical point
We consider QCD in the d = 4 − 2ǫ Euclidean space. The action reads where D µ = ∂ µ − ig B A a µ T a with T a being the SU(N ) generators in the fundamental (adjoint) representation for quarks (ghosts). The bare coupling constant is g B = gM ǫ where M is the scale parameter, and the strength tensor is defined as usual The renormalized action is obtained from (2.1) by the replacement where Z ξ = Z 2 A and the renormalization factors are defined using minimal subtraction where z jk are ǫ-independent constants. Note that we do not send ǫ → 0 in the action and the renormalized correlation functions so that they explicitly depend on ǫ.
Formally the theory has two charges -g and ξ. The corresponding β-functions are where with The anomalous dimensions of the fields Φ = {q,q, A, c,c} are defined as They are known to a high order, O(α 5 s ) for the quark anomalous dimension [12]. In what follows we also use a notation For a sufficiently large number of flavors, N f , one obtains β 0 < 0. Therefore, there exists a special (critical) value of the coupling, g = g * (ǫ) such that β g (g * ) = 0, alias ǫ = −γ g (a * ) or, equivalently, The β-function associated with the gauge parameter ξ vanishes identically in the Landau gauge ξ = 0. As a consequence Green functions of the quark and gluon fields in Landau gauge at critical coupling enjoy scale invariance [13][14][15]. Scale invariance usually implies conformal invariance of the theory: it is believed that "physically reasonable" scale invariant theories are also conformally invariant, see ref. [16] for a discussion. In non-gauge theories conformal invariance for the Green functions of basic fields can be checked in perturbative expansions [17,18]. For local composite operators a proof of conformal invariance is based on the analysis of pair counterterms for the product of the trace of energy-momentum tensor and local operators [19]. In gauge theories, including QCD, conformal invariance does not hold for the correlators of basic fields and can be expected only for the Green functions of gauge-invariant operators. Extra complications are due to mixing of gauge-invariant operators with BRST variations and equation-ofmotion (EOM) operators. We will discuss these issues briefly in what follows.
Renormalization ensures finiteness of the correlation functions of the basic fields that are encoded in the QCD partition function. Correlation functions with an insertion of a composite operator, O k , possess additional divergences that are removed by the operator renormalization, where the sum goes over all operators with the same quantum numbers that get mixed; Z kj are the renormalization factors that have a similar expansion in inverse powers of ǫ as in eq. (2.3). Here and below we use square brackets to denote renormalized composite operators (in a minimal subtraction scheme). Renormalized operators satisfy a RG equation with the anomalous dimension matrix (or evolution kernel, in a different representation) H ∼ (−M ∂ M Z)Z −1 (up to field renormalization) which has a perturbative expansion with the coefficients that in minimal subtraction schemes do not depend on ǫ by construction. As a consequence, the anomalous dimension matrices are exactly the same for QCD in d dimensions that we consider at the intermediate stage, and physical QCD in integer dimensions that is our final goal. Namely, if in d-dimensional QCD at the critical point with the same matrices H (k) . All what one has to do in going over to the four-dimensional world is to reexpress consistently all occurrences of ǫ = (4 − d)/2 in terms of the critical coupling ǫ = β 0 a * + . . . and replace a * → a in the resulting expressions. The requirement of large N f for the existence of the critical point is not principal since, staying within perturbation theory, the dependence on N f is polynomial. In this sense the above connection holds for an arbitrary number of flavors. Conformal symmetry of QCD in d-dimensions at the critical point means that evolution equations in physical QCD in minimal subtraction schemes to all orders in perturbation theory have a hidden symmetry: one can construct three operators that commute with H and form an SL(2) algebra, i.e. they satisfy (exactly) the SL(2) commutation relations. As we will see below, perturbative expansion of these commutation relations produces a nested set of equations that allow one to determine the non-diagonal parts of the anomalous dimension matrices with a relatively small effort. A digression to the 4 − 2ǫ dimensional world, from this point of view, is just a technical trick in order to obtain the explicit expression for one of these operators, the generator of special conformal transformations. To avoid misunderstanding, we stress that QCD in d = 4 dimensions is certainly not a conformal theory. The symmetry that we are going to exploit is the symmetry of RG equations in QCD in a specially chosen regularization scheme based on dimensional regularization with minimal subtraction. The whole construction becomes simpler and more transparent going over from local operators to the corresponding generating functions that are usually referred to as light-ray operators. This representation is introduced in the next section.

Leading-twist operators
Poincare symmetry of the theory is enhanced at the critical point a = a * , β(a * ) = 0 by the dilatation symmetry (scale invariance) and symmetry under space-time inversion. The subject of this work are flavor-nonsinglet twist-two (symmetric and traceless) operators where q andq are quark (antiquark) field operators that we tacitly assume to be of different flavor, D µ is a covariant derivative, and n µ is an auxiliary light-like vector, n 2 = 0. Symmetry transformations that act nontrivially on these operators form the so-called collinear SL(2, R) subgroup of the full conformal group that leaves the light-ray x µ = zn µ invariant, see ref. [20] for a review. Collinear conformal transformations are generated by translations along the light-ray direction n µ , special conformal transformations in the alternative light-like directionn, n 2 = 0, (nn) = 1, and the combination of the dilatation and rotation in the (n,n) plane Explicit expressions for the generators of translations P µ , dilatations D, special conformal transformations K µ and Lorentz rotations M µν can be found, e.g., in ref. [20]. Here and below we use a shorthand notation P n = n µ P µ etc. The generators defined in this way satisfy standard SL(2) commutation relations Local composite operators can be classified according to irreducible representations of the SL(2) algebra. A (renormalized) operator [O N ](x) is called conformal if it transforms covariantly under the special conformal transformation: Here ∆ * N is the scaling dimension of the operator (at the critical point): where γ * N is the anomalous dimension at the critical point, γ * N = γ N (a * ). The scaling dimension is given by the sum of the canonical and anomalous dimensions, For the operators under consideration ∆ N = 2∆ q +N where ∆ q = d/2−1/2 is the canonical dimension of the quark field. In a conformal theory, the correlation function of conformal operators is annihilated by the generator of special conformal transformations, where we added the superscripts K All operators O N k in a conformal tower have, obviously, the same anomalous dimension γ * N .

Light-ray operators
A renormalized light-ray operator, where the Wilson line is implied between the quark fields on the light-cone, is defined as the generating function for renormalized local operators: Due to Poincare invariance in most situations one can use x = 0 in the definition of the light-ray operator (2.25) without loss of generality; we will often use a shorthand notation The renormalization factor Z is an integral operator in z 1 , z 2 which is given by a series in 1/ǫ where H is an integral operator (evolution kernel) acting on the light-cone coordinates of the fields. It is related to the renormalization factor (2.26) as follows where Z = ZZ −2 q . The evolution kernel can be written as [21] H where h(α, β) is a certain weight function (evolution kernel). Here and below we use the notation In perturbation theory h(α, β) is given by a series in the coupling constant It is important to note that the fixed-order kernels h (k) (α, β) in the MS scheme do not depend on the space-time dimension by construction. Thus the dependence of h(α, β) on ǫ in QCD in 4 − 2ǫ dimensions at the critical point a * = a * (ǫ) comes exclusively through the coupling constant.
Going over from the description in terms of conformal towers of local operators to the light-ray operators essentially corresponds to going over to a different realization of conformal symmetry. The light-ray operator [O(x; z 1 , z 2 )] can be expanded in terms of that we will refer to as coefficient functions. The action of the generators of conformal transformations L ±,0 on the light-ray operator is defined via their relation to local operators, where α = ±, 0 and the coefficients ℓ N km α can be read off eq. (2.22). 1 In this way the action of the generators L ±,0 on the quantum fields in the light-ray operator can be traded for the operators S ±,0 acting on the coefficient functions Ψ N k (z 1 , z 2 ): and they can be represented as certain integro-differential operators S ±,0 acting on the quark coordinates in the light-ray operator itself, in particular

JHEP03(2016)142
The generators S ±,0 in this (position-space) representation obey the SL(2) commutation relations (note a different sign as compared to the algebra of quantum operators (2.16)), and commute with the evolution kernel Explicit expressions for the generators in the interacting theory (at the critical point) are nontrivial as, with the exception of S − , they are modified by quantum corrections. One can write them in the following form: Note that quantum corrections to S 0 are completely determined by the evolution kernel H, whereas the generator of special conformal transformation along then direction, S + , contains an additional contribution ∆ + that can be calculated order by order in perturbation theory, By construction, the evolution kernel H and the operator ∆ + commute with the canonical generator S

Conformal constraints for the evolution equation
One can show that the coefficient functions of the operators from the conformal tower are eigenfunctions of the evolution kernel for the light-ray operator and have the following form Thus, the coefficient function of the conformal operator [O N ] is ∼ z N 12 , and the coefficient functions of the operators with extra total derivatives, (n∂) k [O N ], are obtained by the repeated application of the "step-up" operator S + . As a consequence, anomalous dimensions of local operators correspond to the moments of the evolution kernel for the light-ray operator (2.45) Using the representation for the generators in (2.40) and expanding the commutation relation [S + , H] = 0 in a power series in the critical coupling a * one obtains a nested set of equations [10] [S (0) The first equation (2.46a) expresses the usual wisdom that one-loop QCD evolution equations (in four dimensions) respect conformal symmetry of the QCD Lagrangian [22]. In this case it can be shown that the corresponding kernel h (1) (α, β) (up to trivial terms ∼ δ(α)δ(β) that correspond to the unit operator) takes the form [23] and is effectively a function of one variable τ called the conformal ratio. This function can easily be reconstructed from its moments (2.45), alias from the anomalous dimensions. This prediction is confirmed by explicit calculation [21]: The corresponding one-loop kernel h (1) (α, β) can be written in the following, remarkably simple form [23] where the regularized δ-function, δ + (τ ), is defined as =0, is given by the commutator of the one-loop kernel and the one-loop modification of the generator of special conformal transformation [9,11] ∆S (1) Since the canonical generator S (0) + is nothing but the first-order differential operator, eq. (2.46b) can be viewed as the first-order inhomogeneous differential equation for the two-loop kernel H (2) . The general solution of this equation can be found as a special solution of the inhomogeneous equation, corresponding to the symmetry breaking part of the evolution kernel, complemented by a general solution of the corresponding homogeneous equation [S (0) + , H (2) ] = 0 which has to be fixed by the requirement that the moments (2.45) reproduce the known two-loop anomalous dimensions, see ref. [11] for the details. The explicit expression for the two-loop kernel h (2) (α, β) is given in appendix B. It is equivalent to the result for the two-loop splitting functions (in a different representation) for flavor-nonsinglet GPDs derived in [9] by a somewhat different method.
It is easy to see that this hierarchy continues to all orders in perturbation theory: the evolution kernels at a given order of perturbation theory can be obtained from the spectrum of anomalous dimensions at the same order and an additional calculation of the modification of the generator of special conformal transformations at one order less. In particular the three-loop evolution kernels require the knowledge of S + to two-loop accuracy, see eq. (2.46c). The corresponding calculation is the subject of this paper.

Scale and conformal Ward identities
Ward Identities (WI) follow, in general, from invariance of suitable correlations functions under the change of variables in their path-integral representation, corresponding to a symmetry transformation. The standard choice is the correlation function of the composite operator in question with the set of fundamental fields. In gauge theories and in particular in QCD it is more convenient to consider for the same purpose the correlation functions of light-ray operators, which are gauge-invariant.
As mentioned above, the operator S + in the light-ray operator representation is defined as the generator of special conformal transformations in then direction acting on the lightray operator aligned in the opposite n-direction and centered at the origin, x = 0: (3.1) (Here we display explicitly the dependence on the auxiliary vector n in the definition of the light-ray operator). On the other hand, taking instead the n-projection and for arbitrary x such that (x · n) = 0 one gets or, changing n →n, Consider the correlation function of two light-ray operators aligned in opposite light-like directions and separated by a transverse distance (x · n) = (x ·n) = 0: where we use a shorthand notation where the superscript S (z) + reminds that it is a differential operator acting on the z 1 , z 2 coordinates. The explicit expression for S (z) + can be derived starting from the path-integral representation Here N is a normalization factor, S R (Φ) is the renormalized QCD action, Φ = {A, q,q, c,c} and the functional integration goes over all fields. Let us make a change of variables in the path-integral corresponding to the dilatation and special conformal transformations, respectively, see e.g. ref. [20]. Σ µν in (3.8) is the generator of spin rotations, and ∆ Φ are the scaling dimensions of the QCD fundamental fields, which are conveniently chosen as follows [9]:

JHEP03(2016)142
The choice ∆ A = 1 ensures that the nonabelian field strength tensor transforms covariantly under conformal transformations 10) and the rationale for ∆ c = 0 is that for this choice a covariant derivative of the ghost field D ρ c(x) transforms as a vector field of dimension one, i.e. in the same way as the gluon field A ρ . Invariance of the path-integral representation of the correlation function of two lightray operators G(x; z, w) under the change of variables implies the identity where δ = δ D and δ = δ K =n µ δ µ K for scale and conformal transformations, (3.7) and (3.8), respectively, and δS R is the corresponding variation of the QCD action Note that the coefficient of ∂ ρ B ρ (x) in the conformal variation does not vanish for ǫ → 0. Hence the QCD action is not invariant under conformal transformations even for integer d = 4 dimensions. The operator B µ (x) can, however, be written as a BRST variation of c a A a µ [9], see appendix A. Thus this term does not contribute to correlation functions of gauge-invariant operators [24] and can be dropped in most cases, which greatly simplifies the analysis.
In what follows we analyze the structure of the Ward Identities (3.11) in detail.

Scale Ward identity
Let us first consider the scale, or dilatation, WI (SWI). The variation of the renormalized operators on the l.h.s. of eq. (3.11) is given by where we used that the renormalization Z-factor commutes with the operator i z i ∂ z i counting the canonical dimension. (This is nothing but a usual observation that only the operators of the same canonical dimension mix under renormalization.) We obtain Since the l.h.s. can also be written as a derivative over the scale parameter M ∂ M G(x; z, w) and the RG equations for the light-ray operators take the form (2.27), the expression on the r.h.s. of (3.15) that contains the variation of the action δ D S R (3.12) can be written as It is instructive to derive this result by a direct calculation using a method that can be generalized to the more complicated case of the conformal WI (CWI) (see also refs. [7,9]).
The starting observation is that correlation functions of the basic fields with an insertion of the operator N (x) (3.13) are finite, as follows from the structure of the corresponding scale and conformal WIs. Thus N (x) can be expanded in terms of renormalized operators and the coefficients in this expansion can be fixed (apart from certain terms involving total derivatives) from the renormalization group analysis. The result reads [7,9,19,20,24,25] where Ω Φ is an EOM operator, The constants z c (g, ξ) and z b (g, ξ) cannot be determined in this method. In order to make the presentation self-contained we explain their derivation in appendix B. The last term in (3.17), being a BRST variation, does not contribute to the correlation function in (3.15). The ghost EOM terms Ωc and Ω c also do not contribute since the lightray operators do not contain ghost fields, e.g., Further, the gauge fixing term can be replaced by the sum of EOM terms using eqs. (B.17) and (A.5), Note that the coefficients ξ∂ ξ ln Z Φ are given by a series in 1/ǫ without a constant term. Since the WI is finite, all such singular terms must cancel in the final answer and it is sufficient, in principle, to trace the nonsingular terms only. In other words, although the terms in (3.18) contribute to the WI, their only role is to cancel some other singular contributions. It is instructive, nevertheless, to trace these cancellations explicitly. Thus we can replace N (y) by a somewhat simpler expression The EOM contributions give rise to contact terms that can be evaluated integrating by parts in the path-integral (3.20) The quark and the antiquark EOM operators, Ω qq = Ω q + Ωq, give together Gluon EOM contributions are more complicated (because light-ray operators contain terms with an arbitrary number of gluon fields), but we will show that they cancel. The main contribution comes from the insertion of the renormalized Lagrangian (3.23) Since this correlation function involves three renormalized operators, the counterterms corresponding to operator renormalization are already subtracted. All remaining divergences correspond to pair counterterms for the contraction of [L Y M +gf ](y) and one of the light-ray operators, y → 0 or y → x. We can write, schematically, where PCt(A(x)B(y)) denotes the pair counterterm for the contraction of the operators A and B. The first term on the r.h.s. of eq. (3.24) is finite, by definition, so that it does not contribute to (3.23) at the critical point β(a * ) = 0 [18,19]. Pair counterterms for the product of two arbitrary operators A(x) and B(y) have the following general structure [19] PCt where C i (x), C i µ (x) are local operators and Z i , Z i are singular coefficients. The ellipses stand for the contributions with more than one derivative acting on the δ-function. For the case at hand only the terms without derivatives are relevant, which can be found without explicit calculation.
To this end let us compare the structure of divergent contributions in the correlation function of the two light-ray operators with and without the d d y[L] insertion, vs.
where the sum goes over all one-particle-irreducible (1PI) Feynman diagrams, D{. . .} stands for the expression for a given diagram, and the KR ′ operation corresponds to taking singular part after subtraction of divergences in all subgraphs.
An insertion of d d y L Y M +gf in a generic Feynman diagram for the correlation function generates two types of contributions: the kinetic term gives rise to an insertion (of unity) in gluon propagators, and the interaction terms correspond to a replacement of one of the "usual" QCD vertices (three-gluon or four-gluon) by the "special" vertex which is in fact identically the same as the "usual" one. The only effect of these substitutions is an extra combinatorial factor: e.g. the kinetic term can be inserted in any gluon line, thus the original diagram is effectively multiplied by the number of gluon lines, and similar for the vertices. It is easy to convince oneself that the combined effect of all insertions is a multiplication of the diagram by the number of loops (minus one, because the leading-order diagram for the correlation function already contains a loop).
Divergent contributions to the correlation function of the light-ray operators (3.27) obviously correspond to their renormalization. Note that a single light-ray operator contains contributions with arbitrary many gluon fields, ℓ = 0, 1, . . . so that the renormalized light-ray operator takes the form, schematically (3.28) The pair counterterms of interest are given by the same set of diagrams that give rise to the above product of renormalization factors, correcting for their combinatorial factors. Note that a multiplication by the number of loops (in a particular divergent subgraph) amounts to taking a derivative a∂ a . Hence we can write e.g. for the operator counterterm corresponding to the contraction of L Y M +gf (y) and O (n) (x, z) (cf. ref. [9]) Adding the second pair counterterm, and taking into account that We stress that this term does not contribute to the SWI at the critical point. Also, since all singular terms in ǫ have to cancel, we could drop them from the beginning and only consider finite contributions. This cancellation is rather nontrivial on a diagrammatic level. We have demonstrated how it works for the SWI, but we will simply assume of in the analysis of CWI in the next section.

Conformal Ward identity
The two terms on the l.h.s. of the conformal Ward identity (CWI), eq. (3.11), correspond to the variation of the light-ray operators. The first one can be expressed in terms of S + , where S (ǫ) , the term −ǫ(z 1 + z 2 ) is due to the modification of the quark scaling dimension ∆ q = 3 2 − ǫ, cf. (3.9). The product Z S where the ellipses stand for the singular 1/ǫ terms. As discussed above, the explicit expression for the singular contributions is not needed since they must cancel in the final result. It is easy to show that the conformal variation of the second light-ray operator retains its leading order form (for our choice (x ·n) = 0) Thus the CWI takes the form 2(nn)ZS (ǫ) where we have discarded the term due to the BRST operator ∂ ρ B ρ (3.12) as it does not contribute to gauge-invariant correlation functions.
The contribution due to the quark EOM reads Again, the singular terms can be dropped since they must cancel. The next contribution is due to Similar to the case of the SWI, the correlation function can be written as the finite part, plus contributions of pair counterterms corresponding to the contraction of [L Y M +gf (y)] and one of the light-ray operators. The principal difference is that now we need terms involving the first derivative of the delta-function in (3.25). Such terms cannot be written in terms of the evolution kernel H and require a separate calculation. It is easy to see that for the case of O (n) aligned in the same direction as the parameter in the conformal transformation, δ K =n·K, all pair counterterms vanish as the factor (n·y) under the integral inevitably producesn 2 = 0. Thus we only need pair counterterms for the product [L Y M +gf (y)][O (n) ](0, z) which can be calculated considering the Green function O (n) (0, z)q(p)q(p ′ )) with an insertion of the additional vertex d d y 2(n · y)L Y M +gf (y). The corresponding contribution to the correlation function of the two light-ray operators is, for a given Feynman diagram D,   One should expect that in the final answer G 0 (x; z, w) will be substituted by the complete correlation function G(x; z, w) = G 0 (x; z, w) + G 1 (x; z, w) + . . .

JHEP03(2016)142
so that the same integral operator ∆S + (a) appears for the correlation function at any order of perturbation theory. This property does not hold for the pair counterterm contributions alone where, in general, different operators δS + , δS ′ + can appear, δS + G 0 (x; z, w) + δS ′ + G 1 (x; z, w) + . . . , (3.42) etc., and it has to be restored by adding the gluon EOM contributions ∝ (γ A + γ g )Ω A . 2 Summing all contributions and taking into account that the first contribution in (3.38) vanishes at the critical point and all singular terms 1/ǫ must cancel, the CWI (3.11) takes the expected form, (3.5), where the operator of special conformal transformation S + is given by the following expression: Here γ * q = γ q (a * ) and it is understood that the shift in the space-time dimension ǫ = (4 − d)/2 is written as an expansion in terms of the critical coupling, ǫ = ǫ(a * ). The role of the term γ * q − ǫ (z 1 + z 2 ) is to shift the conformal spin of the quark field to its correct value at the critical point, where the operator ∆ + commutes with S − and anticommutes with the permutation operator of quark coordinates P 12 f (z 1 , z 2 ) = f (z 2 , z 1 ), The role of the term −γ * q (z 1 + z 2 ) is to cancel the corresponding term in (3.43), (3.44) such that the (gauge-dependent) quark anomalous dimension falls out of the final answer. We will see that the structure (3.45) indeed arises naturally in the calculation.

Technical details
In this section we present technical details of the calculation. The problem reduces to the calculation of singular contributions to the 1PI Feynman diagrams for the Green function Note that the counterterms corresponding to the renormalization of the light-ray operator and the Lagrangian insertion (that we do not need) are disposed of by the usual R-operation. It is convenient to rewrite Here L Y M int (y) contains the three-and four-gluon interaction vertices and The last term in (4.1) can be represented in the form A µ (∂A) = ξ(cD µ c − B µ ). The operator B µ is a BRST variation and, hence, does not contribute to the correlation function.
Omitting this term one gets for (4.1) where the gray boxes at the endpoints stand for the multiplication by the corresponding coordinates, (n · x) and (n · y), respectively. Such insertions violate translation invariance and the main trick is to move them either to the external quark lines or to the quark positions in the light-ray operator in which case the corresponding singular contributions can be related to the evolution kernel. Examples will be given below. A shift x → y corresponds to the simple rewritinḡ n · x =n · (x − y) +n · y where the gluon propagator with a thick arrow (in Feynman gauge) is defined as 3 It is sometimes convenient to move the coordinate insertion along the quark lines and/or along the gauge link so we also introduce notations

One-loop calculation
Only the two-gluon effective vertex, − 1 2 (n · x)A µ K µν A ν , is relevant to this accuracy. There are two one-loop diagrams, shown in figure 1, and the diagram symmetric to the one in figure 1b with the gluon attached to the quark field. The addition of such symmetric contributions is always implied.
The first diagram, figure 1a, corresponds to the attachment of the (n · x 1 ) or (n · x 2 ) factor to the external (anti)quark line,

·
The contribution of such a diagram to δS respectively. The factor 1/2 is due to the definition of H that involves a derivative a∂ a of the corresponding Z-factor, see eq. (2.28). For a generic ℓ-loop diagram D of this type one gets for the sum of contributions with (n · x k ) attachments to the quark and the antiquark line where H (1) (b) stands for the corresponding contribution to the evolution kernel. Taking into account the symmetric contribution with the gluon attached to another quark line, and adding the contribution of the diagram in figure 1a, eq. (4.10), we obtain for the sum of these terms with the complete one-loop evolution kernel (2.48), which is exactly the anticipated first term in eq. (3.45). The third term involves the insertion of the gluon and quark coordinates in the gauge linkn · (z u 21 n − z 2 n) = uz 12 (nn) .

JHEP03(2016)142
The same diagram without this insertion (i.e. without the thick arrow) gives rise to the contribution to the evolution kernel [21]: The characteristic structure ∼ [f (z 1 , z 2 ) − f (z α 12 , z 2 )] corresponding to the "plus" distribution in momentum space can be traced to the integration over the gluon position on the light-cone such that the above answer arises from the representation where u is the gauge link variable as in (4.12), as an intermediate step. The diagram "with an arrow" is given, therefore, by the same expression with an insertion of uz 12 (nn) and adding the factor 1/2 due to a different normalization The symmetric diagram with the gluon attached to the quark instead of the antiquark gives the same contribution up to a replacement z 1 ↔ z 2 so that in the prefactor z 12 → −z 12 . The symmetric contribution is, therefore, effectively subtracted and one obtains in the sum δS (1) reproducing the result quoted in (2.51) [11]. It is easy to show that this kind of relation between the diagrams with and without the arrow on the gauge link is general and true in all orders, the reason being that integration over the position of the gluon field in the light-ray operator does not interfere with the separation of singular parts. For a generic ℓ-loop contribution of this type one obtains with the same function h (ℓ) (D) (α). Thus diagrams with an arrow on the gauge link do not require a separate calculation as well.

Two loop calculation
To two-loop accuracy we have to have to take into account the Feynman diagrams of 16 different topologies shown in figure 2. Each of these diagrams gives rise to several contributions to ∆S + corresponding to a replacement of either one of the gluon propagators by the effective propagator (4.4), or of the three-gluon (and quark-antiquark gluon) vertex by the corresponding effective vertex (4.3). Below we tacitly imply using Feynman gauge, ξ = 1.
The "QED type" diagrams in figure 2 (a-c), (e), (f), (h-j), (l-o) can be calculated using the same strategy as the one-loop diagrams considered above. Let us illustrate this procedure on the example of the diagram figure 2(j). Inserting the effective gluon JHEP03(2016)142 · ¾ · · · · · · Figure 3. Rearrangement of three-gluon vertex insertions combined with effective propagators.
propagators one obtains four contributions: · · · · · ¾´ ½µ´ ¾µ · 3 ¿µ´ µ · The first and the last one, (j1) and (j4), combine to where H (j) is the contribution of the diagram in figure 2(j) to the evolution kernel. The remaining diagrams (j2) and (j3) with a modified gluon propagator (4.7) have to be calculated explicitly. The result can be found in appendix C.
The diagrams which contain the three-gluon vertex can be handled in the following way. It is convenient to consider the sum of contributions with the effective vertex and effective gluon propagators and rearrange it as shown in figure 3. Here the gray blobs in the diagrams in the first row stand for insertions generated by (4.3). Their sum can be rewritten as shown in the second row where the white box with an arrow denotes a new vertex gf abc g µρnν − g µνnρ . (4.21) Note that this vertex is symmetric under the interchange of the lower pair of gluons, (ν, b) ↔ (ρ, c), but the line (µ, a) is distinguished, hence an arrow in the notation. This special direction has to be chosen in such a way that the contributions with the insertion of (n·x k ) factors (gray boxes) in the external lines combine to produce a term ∼ H (D) (z 1 +z 2 ). For example, the contribution of the diagram in figure 2(k) can be split in the following five terms: The first two contributions give rise to the term 1 2 H (k) is the contribution of the diagram in figure 2(k) to the evolution kernel, and the remaining three have to be calculated explicitly, see appendix C.
Finally there are four diagrams with self-energy insertions, figure 2(a),(b),(h),(i). It is easy to see that · and also ¾ · ¾ · where the dark oval corresponds to the sum of the contributions of quark, gluon and ghost loops. Using these replacement rules one obtains immediately, e.g., for the diagrams in figure 2(h) and figure 2(i), respectively, where H (h,i) are the corresponding contributions to the evolution kernel.

JHEP03(2016)142
Finally, one has to consider insertions of (4.3) in the self-energy blob itself. For the gluon loop one obtains · · · · (4.23) It can be checked, however, that this contribution is cancelled identically by the similar diagrams with the ghost-gluon vertex insertions and the insertions in ghost propagators so that the insertions inside self-energy diagrams can be omitted altogether. The rest of the calculation is relatively straightforward. The complete results for the contribution of each Feynman diagram in figure 2 to the conformal anomaly and the evolution kernel are presented in appendix C.

Final results
Let us start with the evolution kernel The one-loop result reads [21] H (1) where f (z 1 , z 2 ) is a test function, and the two-loop kernel [11] can be written in the form Here P 12 f (z 1 , z 2 ) = f (z 2 , z 1 ) is the permutation operator and 4 where τ = αβ/ᾱβ. The function ϕ(α, β) is defined as The corresponding one-and two-loop anomalous dimensions are well known and can be found, e.g., in ref. [1]. Next, the generator of special conformal transformations reads where (3.47)

∆S
(1) The one-loop "conformal anomaly" contribution, ∆ + , is very simple [7,11] ∆ (1) The expression for the two-loop anomaly, ∆ + , represents our main result. It can be written as The kernels κ(t), ω(α, β), ω P (α, β) receive contributions of three different color structures Alternatively one can write the results separating the contributions of planar diagrams and the non-planar 1/N c suppressed corrections where, obviously, κ P = κ F F +2κ F A etc., and we took into account that the terms involving quark permutations on the light cone do not receive planar contributions, ω P F A = − 1 2 ω P F F . Note that the term ∼ β 0 in κ(t) arises by choice, rewriting the contribution proportional to the number of quark flavors N f in terms of β 0 . This rewriting is not mandatory and is only motivated in the present context by resulting in somewhat simpler expressions. In contrast, the terms (z 1 + z 2 )β ℓ−1 in the expression for ∆S (ℓ) + involve the "genuine" QCD β-function.
Explicit expressions for the "two-particle" kernels ω, ω P are: and for the planar combination For the "one-particle" kernels κ(t) we obtain The last expression can also be rewritten as The result for κ F F (t) can easily be obtained by subtracting κ P (t) − 2κ F A (t). Note that we prefer to write the corresponding contribution to ∆ + (second line in (5.10)) as a nested integral in auxiliary u, t variables. One of these integrations can be taken trivially (cf. (5.9)) resulting in somewhat more complicated expressions involving Li 3 -functions. Finally, the commutator term 1 4 H (2) , z 1 + z 2 in eq. (5.8) can be written as This term can be added to the result for ∆ (2) + (5.10) but keeping it separate seems to be more convenient for applications.

Conclusions
QCD evolution equations in minimal subtraction schemes have a hidden symmetry: one can construct three operators that commute with the evolution kernel and form an SL(2) algebra, i.e. they satisfy (exactly) the SL(2) commutation relations. In this paper we find explicit expressions for these operators to two-loop accuracy. On this way we make a digression to the 4 − 2ǫ dimensional world, where conformal symmetry of QCD is restored on quantum level at the specially chosen (critical) value of the coupling, and at the same time the theory is regularized allowing one to use the standard renormalization procedure for the relevant Feynman diagrams. We want to emphasize that the procedure is valid to all orders in perturbation theory and the result obtained in this way is complete, i.e. it includes automatically all terms that can be identified as due to a nonvanishing QCD βfunction (in the physical theory in four dimensions). To avoid misunderstanding, we stress that QCD in d = 4 dimensions is certainly not a conformal theory. The symmetry that we uncover is the symmetry of RG equations in QCD in a specially chosen (dimensional) regularization scheme.
It is well known that (QCD) has a non-trivial fixed point in strictly four space-time dimensions for a range of values of the number of quark flavors 9 ≤ N f ≤ 16 known as the conformal window, and in the last years there has been increasing interest in the study of the phase structure of such theories (Banks-Zaks fixed point [13]) on the lattice, see e.g. [26] for a recent review. Our result for the generators of conformal transformations, where one has to use the appropriate values of N f and N c , is valid for QCD within the conformal window as well, to the O(a 2 * ) accuracy in the critical coupling. The main motivation for this study is to obtain three-loop evolution equations for the generalized hadron parton distributions and light-cone meson distribution amplitudes that are relevant for the large-scale experimental studies of hard exclusive reactions in the coming decade. The present work presents the first step in this direction. The remaining calculation can be done in several ways. One possibility is to solve the system of linear differential equations (2.46c), as demonstrated in ref. [11] to the two-loop accuracy. Alternatively, one can exploit the well-known observation [27] that the evolution kernel must JHEP03(2016)142 be a function of the quadratic Casimir operator of the collinear conformal group. This function can be found from the spectrum of anomalous dimensions. Yet another possibility is to bypass the explicit construction of the kernel and try to find directly the solutions (conformal operators) by constructing a unitary transformation U that brings the generator S + to its canonical form, U S + U † = S (0) + . Utility of each of these methods requires a separate study that goes beyond the scope of this work.
Last but not least, the explicit perturbative construction of the generators of conformal transformations can be interesting in the context of the AdS/CFT correspondence for the maximally supersymmetric N = 4 Yang-Mills theory. The general structure of this expansion should be similar to what is obtained in this paper for QCD, but the answer is expected to be simpler. It would be interesting to do this calculation and compare the result with the algebraic approach in ref. [28].
Thanks to this identity the gauge fixing term in the action can be represented as a sum of the EOM and BRST exact operators Importantly, the matrix Z KK ′ has an upper triangular form. Thus

B Renormalization group analysis
The last two terms do not contribute, however, to the correlation function of two gauge invariant operators at different space-time points x = y, Indeed, class C, EOM operators, give rise to contact terms ∼ δ (d) (x−y), whereas operators of class B do not contribute to the correlation function due to the BRST invariance of the QCD action.
In this work we consider flavor-nonsinglet leading twist operators made of the quark and antiquark field and covariant derivatives. For these operators any counterterms (of any class) contain the same pair of quark fields,q, q. A simple dimensional analysis shows, however, that possible operators of class B and C cannot, in this case, be traceless in all Lorentz indices -so that they cannot appear as counterterms in leading-twist operators. Thus only gauge-invariant operators can contribute, The scale and conformal variation of the action (3.12) contains the "symmetry breaking" operator N (3.13). Making use of eq. (A.5) we can write it in the form Note that the last two terms drop out from the correlation functions of gauge invariant operators and N . Our goal is to express the operator N in terms of renormalized operators. From the general operator mixing pattern discussed above we expect the following structure: Next, consider for a moment the Ward Identities of the type (3.11) for the products of (renormalized) fields Φ k . Since the terms involving the variations of the fields are finite (they reduce to a differential operator applied to a finite Green function, cf. (3.15)), the term with the variation of the action must be finite as well; hence the operator N is also finite (up to possible terms containing two total derivatives) and, therefore, the product ǫ Z K is finite for all factors Z K appearing in (B.6). Taking into account that Z F , Z B = 1 + O(1/ǫ), whereas all other factors only contain poles, where the coefficients r K do not depend on ǫ; they are functions of the coupling constant and the gauge fixing parameter. These coefficients for the operators that do not involve total derivatives (alias whose matrix elements do not vanish for zero momenta) can be fixed from the study of the differential vertex operator insertions, see below. Note that Ωq−Ω q = ∂ µq γ µ q and Ωc−Ω c = ∂ µ [cD µ c − ∂ µc c] = ∂ µ Ω µ do reduce to total derivatives, so that we can determine the sum of the coefficients, r q + rq and r c + rc, but not their difference. Invoking charge symmetry arguments one can argue that r q − rq = 0, whereas r c − rc and the coefficient r B µ cannot be determined in this approach.
As well known, derivatives of Green functions of fundamental fields with respect to the couplings give rise to zero momentum insertions of the differential vertex operators and similarly for g∂ g → ξ∂ ξ . Since the expression on the l.h.s. is finite, one concludes that the correlation function with an insertion of g∂ g L R or ξ∂ ξ L R at zero momentum is also finite. Thus both g∂ g L R and ξ∂ ξ L R (up to total derivatives) can be written as a sum of renormalized operators with finite coefficients. For further analysis it is convenient to redefine, temporary, the bare gluon field A 0 → G 0 = g 0 A 0 so that Taking into account that ∂ g Φ 0 = Φ 0 ∂ g ln Z Φ , in particular for the redefined gluon field ∂ g G 0 = G 0 ∂ g ln(gZ g Z A ), and also together with ∂ g g 0 = −ǫg 0 /β(g) , ∂ g ξ 0 = 2ξ 0 ∂ g ln Z A , ∂ ξ ξ 0 = ξ 0 1 + 2∂ ξ ln Z A , (B.12) one obtains where L Y M R and L gf R are the gauge (Yang-Mills) and the gauge-fixing parts of the (renormalized) QCD Lagrangian.
The expressions on the r.h.s. of the two equations in (B.13) have the following structure: where the ellipses stand for a series in 1/ǫ. Since the operators on the l.h.s. are finite, the addition of these terms (ellipses) effectively amounts to a subtraction of divergences so that the sum is nothing but, by definition, a renormalized operator in MS scheme. Thus where D g = β(g)∂ g , and also We remind that these results are valid for zero-momentum insertions, or, equivalently, upon integration d d x over all space-time points.

JHEP03(2016)142
Finally, note that γ Φ = M ∂ M ln Z Φ = (β g ∂ g + β ξ ∂ ξ ) ln Z Φ and β ξ = −2ξγ A so that D g ln Z A = γ A (1 + 2ξ∂ ξ ln Z A ). Thus the last term in (B.16) can be rewritten as − 1 2ξ [(∂A) 2 ] − Φ Ω Φ ξ∂ ξ ln Z Φ and collecting all contributions we obtain where the ellipses stand for total derivative operators. From this expression one can read the results for the coefficients r K defined in eq. (B.7): r F = γ g , r B = r A = γ g + γ A , r q = rq = γ q , r c + rc = 2γ c + γ g + γ A . (B.19) To avoid misunderstanding note that in the derivation we did not use criticality so that the result in eq. (B.18) is valid for arbitrary coupling.
C Results for separate diagrams in Feynman gauge