Renormalization of twist-four operators in light-cone gauge

We compute one-loop renormalization group equations for non-singlet twist-four operators in QCD. The calculation heavily relies on the light-cone gauge formalism in the momentum fraction space that essentially rephrases the analysis of all two-to two and two-to-three transition kernels to purely algebraic manipulations both for non- and quasipartonic operators. This is the first brute force calculation of this sector available in the literature. Fourier transforming our findings to the coordinate space, we checked them against available results obtained within a conformal symmetry-based formalism that bypasses explicit diagrammatic calculations and confirmed agreement with the latter.


Introduction
The leading power approximation to QCD processes with large momentum transfer, such as the deep-inelastic and deeply virtual Compton scattering, admits an intuitive probabilistic description in the framework of the Feynman parton model [1]. According to the latter, physical cross sections are expressible in terms of (generalized) parton distribution functions. The QCD improved picture arises via systematic inclusions of quantum corrections to probe-parton scattering amplitudes as well as renormalization effects of leading twist Wilson operators that parametrize Feynman parton densities. More subtle effects arise from power-suppressed contributions to hadronic cross sections since they encode information on interference of hadronic wave functions with different number of partons. On the one hand, these are of interest in their own right since they provide access to intricate QCD dynamics [2]. On the other hand, they can be regarded as a QCD contaminating background to high precision measurement of New Physics, see, e.g., [3].
In either case, understanding these contributions quantitatively is indispensable at the precision frontier. Since data is typically taken at different values of the momentum transfer, at some point one has to incorporate effects of logarithmic scaling violation stemming from renormalization of higher twist operators. The task of their unravelling at twist-four level will be undertaken in the present study. Until very recently, only partial results for certain subsets of operators were available in the literature [4]. A special class of operators out of all higher twists is known as quasipartonic. They can be characterized either as composite fields built from on-shell fields of the Feynman parton model or understood as operators with their twist equal to their length, i.e., the number of fields that form them. For this class of operators a systematic approach to constructing high-twist evolution equations was developed about three decades ago by Bukhvostov, Frolov, Lipatov and Kuraev in Ref. [5]. At leading order in QCD coupling, the evolution kernel for these was found as a sum of pairwise interaction kernels between elementary fields comprising the operators in question. The particle number-preserving nature allows one to map it to a Hamiltonian quantum mechanical problem. This advantage was explored in a number of works at twist three level [6] 1 starting from [8]. Eventually, the problem was mapped into an exactly solvable lattice model [9,10]. However, while the quasipartonic operators form a subset closed under the renormalization group evolution [5,11], they do not exhaust the set of all operators contributing at a given twist. The remaining ones are dubbed non-quasipartonic and they contain at least one bad field component in the formalism of light-cone quantization. These operators are characterized by the property that their twist is greater than their length. Their evolution does not preserve the number of fields in quantum transitions and thus their study is more elaborate. In the twist three case alluded to above, this was not a pressing issue since the use of QCD equations of motion allows one to remove all non-quasipartonic operators from the basis. For even higher twists, this is not sufficient and particle number changing transitions involved in the analysis of non-quasipartonic operators have to be addressed explicitly.
The analysis of the renormalization problem for twist-four operators was completed recently in the coordinate space [12], i.e., in terms of light-ray composite operators. The formalism is based on the use of conformal symmetry preserved by leading order QCD evolution equations, Poincaré transformations in the transverse plane and a minimal input from Feynman graphs. Presently we will perform a brute-force computation of Feynman diagrams in the light-cone gauge and rely on the momentum-space technique which makes the underlying calculation rather straightforward. For an exception of a few subtleties with the use of QCD equation of motion to recover the particle-number increasing transitions, it reduces to a few algebraic, though rather tedious, steps.
The choice of the operator basis at higher twists is not unique due to multiple relations among a redundant set of operators via QCD equations of motion. Thus it is driven by requirements of simpler transformation properties under residual (conformal) symmetry as well as simplicity of underlying calculations. In the current work we will adopt the basis of twist-four operators suggested in Ref. [12]. This will allow us to verify our results obtained by an independent calculation based on a different technique. Since we will focus on the twist-four sector, we have three types of building blocks at our disposal as two-particle elements of operators in question: good-good, good-bad and bad-bad light-cone field components. According to traditional classification, they possess twists two, (at least) three and (at least) four, respectively. We will address only the first two types, since the last one can be eliminated in hadronic matrix elements in favor of the other ones containing more fields via QCD equations of motion, as discussed below. Our consideration will be limited to QCD nonsinglet sector, though partial results for two-to-two transitions will be reported for the singlet sector as well.
Our subsequent presentation is organized as follows. In the next section, we spell out the operator basis used in the current calculation and provide a dictionary between the twistor notations adopted in Ref. [12] and the light-cone conventions used in the present analysis. Then, we discuss the general structure of twist-four evolution equations and provide a Fourier transform bridge between the light-ray and momentum fraction space representations. In Sects. 4.1, 4.2 and 4.3, we present evolution kernels for two-to-two quasipartonic, non-quasipartonic and twoto-three transitions, respectively. As a result of this analysis we find a simplified form of light-ray evolution kernels for certain evolution kernels which are reported in the Appendices. The latter also contain technical details on the calculation of Feynman diagrams defining operator mixing as well as singlet two-to-two transitions.

Operator basis
The light-cone dominated processes are parametrized by matrix elements of composite operators built up by fields localized on a light-cone ray defined by the vector n µ = (1, 0, 0, 1)/ √ 2 that is reciprocal to the large light-cone component of the momentum transfer. Thus they have the following generic form where the X-field cumulatively stands for certain components of quark and gluon fields as explained below. The positions z − k =n · z k of the fields on the light-cone are defined with the help of a tangent null vectorn µ = (1, 0, 0, −1)/ √ 2 to the light-cone normalized such that n ·n = 1. The gauge invariance of O is achieved by means of an appropriate contraction of the color indices I k (either in the (anti-)fundamental I k = i k or adjoint representation I k = a k of the color group) into an SU(N) singlet with a tensor C I 1 I 2 ...I N and field coordinates parallel transported to an arbitrary position z − 0 with the help of the Wilson lines Here is the light-cone projection of the gauge field.

Good and bad light-cone fields
It is well-known that the light-cone gauge has a number of advantages. First, we observe that the gauge links are gone in Eq. (1) and, as a consequence, this results in reducing of the number of diagrams contributing to loop amplitudes. Second, the Feynman parton model arises naturally from the light-cone gauge QCD. Namely, one decomposes the quark Ψ and gluon fields A µ , in terms of the good X + = {Ψ + , A ⊥ ,Ā ⊥ } and the bad X − = {Ψ − , A − } components, respectively. Note that for the vector A µ we defined its minus as follows projection and, in addition, we decomposed the transverse gauge field in terms of its anti-and holomorphic componentsĀ with the help of the vector e µ ⊥ = (0, −1, −i, 0)/ √ 2 (and its complex conjugateē = e * ). These possess helicity h = ±1, respectively, being eigenvalues of the helicity operator [13] that is built from the spin tensor Σ µν entering the Lorentz generators iM µν . The bad components being non-dynamical in the light-cone time z + can be integrated out in the path integral and, thus, only the on-shell propagating modes Ψ + , A ⊥ andĀ ⊥ are left. We will not perform this step however and keep all non-propagating degrees of freedom in the QCD Lagrangian since the classification of operators will be easier in this case and moreover one does not loose Lorentz covariance. Finally, it is straightforward to construct an operator basis making use of the above building blocks, namely, the field X in Eq. (1) will have the following components (as well as their Hermitian conjugates X † ) with D ⊥ = e ⊥ · D being the holomorphic covariant derivative D µ = ∂ µ − igA µ .

Twistor representation
To make a connection to the basis of Ref. [14], let us recall the twistor formalism used there. We pass to the spinor representation for Lorentz vectors by contracting them with the twodimensional block σ µ of four-dimensional Dirac matrices in the chiral representation γ µ = antidiag(σ µ , σ µ ), e.g., where σ µ = (1, σ) and σ is the three-vector of Pauli matrices, whileσ µ = (1, −σ). The light-cone vectors n andn can be factorized into two twistors λ α and µ α n αα = λ αλα ,n αα = µ αμα , where λ * α = λα and µ * α = µα. For the light-cone vectors introduced in the previous section, we can choose the two-dimensional spinors as λ α = (0, 4 √ 2) and µ α = ( 4 √ 2, 0). These twistors will allow us to construct good and bad fields for specific helicities. Namely, using the decomposition of the Dirac quark field in chiral representation we can introduce their good and bad components as follows Identical relations hold for χ upon the obvious replacement ψ → χ. Here we introduced the bra and ket notations for undotted and dotted SL(2) indices, |λ = λ α , λ| = λ α and |λ] = λα, [λ| =λα that allow us to uniformly contract undotted indices from upper-left to lower-right and dotted ones from lower-left to upper right, i.e., λψ = λ α ψ α and [μψ] =μαψα. In a similar fashion, the gluon field strength can be decomposed as in terms of its chiral f αβ = i 4 σ µν αβ F µν and anti-chiralfαβ = i 4σ µναβ F µν components with the help of the self-dual σ µν = i 2 [σ µσν −σ νσµ ] and anti-self-dual tensorsσ µν = i 2 [σ µ σ ν −σ ν σ µ ]. The plus and minus fields are found by projections etc. Finally, as any four-vector, the covariant derivatives are decomposed in twistor components as follows

Bridging light-cone and twistor projections
The notations introduced in this and preceding sections allow us to establish a dictionary between the light-cone and twistor components. They are summarized by the following set of relations: for fermions, where γ 5 = diag (1, −1), and for gluons and, finally, covariant derivatives, Making use if these conversion formulas, we can adopt the basis introduce in Ref. [14], on the one hand, and use the momentum-space technique of Ref. [8] that makes the calculation more concise while using conventional four-component notations for Lorentz vectors and Dirac matrices.

SL(2) invariance and basis primary fields
Though massless QCD is not a conformal theory at quantum mechanical level since it induces a scale due to dimensional transmutation, the classical Lagrangian of the theory does enjoy SO (4,2) invariance. The one-loop evolution equations that we are set to analyze in this work inherit the latter since the symmetry breaking graphs do not make their appearance till twoloop order. Since the light-cone operators (1) involve fields localized on a light ray, the full conformal algebra reduces to its collinear conformal SL(2) subalgebra that acts only on the minus projections z − k ≡ z k of the Minkowski space-time coordinates z µ k . The differential representation of generators acting on the space spanned by the composite operators (1) reads The irreducible representations are characterized by the value of the conformal spin j n = (d n + s n )/2 determined by the canonical dimension d n and light-cone spin s n of the constituent fields X n . These generators commute with the generator of helicity introduced in Eq. (8) as well as twist E = N n=1 (d n − s n )/2, see, e.g., [15,16]. The field projections introduced in the previous section transform covariantly under SL (2) transformations and can be organized into "multiplets" of the same twist. Namely, the good Φ + and bad Φ − chiral fields as well as their conjugate anti-chiral analoguesΦ ± = Φ * ± , possess twist E = 1 and E = 2, respectively.
Since covariant derivatives D ++ , D ±∓ and D −− carry twist zero, one and two, respectively, they can be used to generate additional high-twist "single-particle" fields by acting on Φ ± . Obviously, we can ignore D ++ since they just induce an infinitesimal shift along the light cone. The D −− derivatives acting on Φ + will produce a twist-three constituent, which when accompanied by another good field component, will form a twist-four operator. However, this operator can be safely dropped from the basis thanks to QCD equations of motion. Next, the transverse derivatives D ±∓ can act either on good or bad fields. However, we can consider only their action on the former since we can always move derivatives from bad to good fields in a twist-four operator of the type Φ + D ±∓ Φ − . Moreover, since it is desirable to deal with conformal primary fields as individual building blocks, as they obey simple SL(2) transformations (22) and thus yield evolution equations with manifest conformal symmetry, one can reduce in half, as advocated in Ref. [14], the basis of good fields with transverse covariant derivatives acting on them. This is achieved by eliminating the ones with non-covariant transformation properties. The net result is that one introduces instead conformal primaries D −− Φ + which can be safely neglected as alluded to above. This procedure allows us to trades D −+ Φ + posing bad SL(2) transformation properties in favor of the primary D +− Φ + . Finally, two transverse derivatives acting on Φ + can be again be reduced to the irrelevant primary D −− Φ + . This concludes the recapitulation of the reasoning behind the choice of the twist-one X + and twist-two X − primaries which build up the operator basis at twist four. The latter is thus spanned by quasipartonic and nonquasipartonic operators (that read schematically) respectively.

Evolution equations
The twist-four light-ray operators (25) mix under renormalization. Their mixing matrix admits perturbative expansion in strong coupling α s = g 2 /(4π). The goal of our analysis is to calculate the leading term of the series, namely, Notice that the mixing matrix takes a triangular form (to all orders in coupling) since the quasipartonic operators form an autonomous set under renormalization group evolution. Here the transition kernels are some integral operators that shift fields on the light-cone towards each other. Their form is highly contained by conformal symmetry and was the subject of recent analysis [12]. The distinguished feature of nonquasipartonic operators is that they can change the number of fields upon evolution. Thus, while H (N →N ) for N = 3, 4 is merely given by the sum of pairwise transition kernels, the kernel H (3→4) involves both two-to-two and two-to-three transitions where σ is a positive/negative integer reflecting the mismatch in the field dimension in a given operator transition. The manifest SL(2) covariance on conformal primaries building up the composite operators and the fact that one loop transitions do not receive contributions from counter terms that break conformal invariance implies the commutativity of the kernels with the generators of the algebra and thus impose sever constraints on the form of the kernels.

Renormalization in momentum space
Though the conformal symmetry is more explicit in the coordinate space, the acutual calculations of one-loop graphs determining the mixing matrix are far more straightforward and elementary in the reciprocal momentum space. As we pointed out in the introduction, a technique to perform this analysis is available for quasipartonic operators from Ref. [5]. Presently we will get it generalized to nonquasipartnic operators as well. The formalism relies heavily on the light-cone gauge, where the gluon propagator takes the form As we can see, the integration over the loop momentum k decomposed in Sudakov variables will potentially produce divergences in the longitudinal variable k + due to an extra pole in the propagator, in addition to the conventional ultraviolet singularities regularized by a cut-off µ in the transverse momentum. The former arise due to incomplete gauge fixing by the condition (4) that allows one for light-cone time-independent residual transformations. To complete gauge fixing, one has to impose a condition on how to go around the 1/k + singularity. This will not be a pressing issue for the current work since we will focus on kinematics away from the phase space boundaries where these have to be treated properly. Let us point out however, that the advanced/retarded and symmetric boundary conditions on the gauge potential on the light-cone infinity impose ±i0 and principal value prescription [17]. While the condition consistent with the equal-time quantization yields the Mandelstam-Leibbrandt prescription [18]. Thus, we convert the light-ray operators to the momentum fraction space by means of the Fourier transform Then, the evolution kernels arise from the N to M-particle transition amplitude, where G(k 1 , . . . , k n |p 1 , . . . , p M ) is a sum of corresponding Feynman graphs. Extraction of the leading logarithmic divergence from the momentum integrals results in the sought momentumspace evolution equation where integral kernel in the momentum representation is with a notation introduced for the measure The N − 1 momentum integrals in (38) are eliminated by means of four-momentum conserving delta functions stemming from Feynman rules leaving us with a single four-dimensional loopmomentum integration measure (35). The extraction of 1/k 2 ⊥ contribution in the integrand can be easily achieved by rescaling the k − component of the loop momentum by introducing a new variable where k ⊥ is the transverse loop momentum. We remove the k + integral making use of the momentum fraction delta functions in Eq. (38), while the rescaled k − -integrals are evaluated in terms of the generalized step-functions [5,20] ϑ k α 1 ,...,αn (x 1 , . . . , These can be reduced to the simplest one ϑ 0 ) making use of a set of known relations [5,20]. So the advantage of this formalism is that there are no actual integrals to perform and the procedure of computing the evolution kernels is reduced to straightforwards but tedious algebraic manipulations with Dirac and Lorentz structures.
So all we need for the calculation is Fourier transforms of the conformal primary fields defining composite operators. When cast in four-dimensional light-cone notations, they read The results for the antichiral fieldsΦ are obtained from above by complex conjugation.

From coordinate to momentum space
It is straightforward to relate evolution kernels in the coordinate and momentum space by a Fourier transformation. For two-to-two transitions, we immediately find which are subject to the momentum-fraction conservation condition x 1 + x 2 = y 1 + y 2 . Analogously, for two-to-three transitions, we find where we assume that x 1 + x 2 = y 1 + y 2 + y 3 . Having results in one representation, one can easily obtain the other making use of the following two elementary operations where η = x 1 + x 2 = y 1 + y 2 implies momentum conservation. Finally, the coordinate kernels possess integrable end-point singularties, that have to be regularized in the course of the Fourier transform. We found the cut-off regularization to be the simplest one to handle the emerging divergencies At the end of the calculation, all singular terms cancel in the limit ε → 0 rendering the total result regular. We provide an example of explicit transformation in Appendix C.

Conformal symmetry in momentum space
Since conformal invariance played a crucial role for the coordinate-space calculations [12], let us analyze its consequences in the momentum space. To this end, following the same reasoning leading to the expression of Eqs 48 and taking care of the ordering of z and ∂ z , one observes the following identifications between the light-ray coordinates and the momentum fractions, where x n is the reciprocal momentum to the coordinate z n . Thus the conformal generators shown in Eq. (22) can be rewritten in the momentum space as Imposing the condition of commutativity we find that the evolution kernels obey the following differential equations (away from kinematical boundaries) x n |y 1 , . . . , y M ) = 0 , x n |y 1 , . . . , y M ) = 0 , for S + and S 0 , respectively. In arriving at the expressions of Eqs. (58) and (59), we have made use of the momentum-conserving Dirac delta function factorized from (40). Since essentially there are only M + N − 1 independent variables in the game, we are required to rewrite one of the variables as a linear combination of the rest before differentiation. Above we chose to eliminate x N . Finally, the S − simply yields the momentum fraction conservation condition which is trivially obeyed due to an overall delta function (41) that accompanies the transition kernel, Since we are interested in two-to-two and two-to-three transitions in this work, the expressions in Eqs. 58 and 59 simplify to where we use j n and j n ′ to refer to the conformal spins of the incoming and out going particles, respectively. Similarly, for three-particle transitions, we get

One-loop kernels
In this section we will report on our findings of all nonsinglet transition kernels. The latter will be quoted away from kinematical boundaries, i.e., when some of the momentum fractions (or their sums) could coincide. This will be sufficient to compare our results with the Fourier transform of the light-ray evolution kernel derived in Ref. [12] by dropping all contact, i.e., deltafunction, terms emerging from the latter. Of course, we can keep track of the latter as well and reproduce them from the momentum-fraction formalism by properly incorporating QCD field renormalization (as well as certain contact terms stemming from vertex graphs) into the game.
Since the light-cone gauge explicitly breaks Lorentz symmetry, the good and bad components receive different renormalization constants as can be immediately seen from the quark and gluon propagators [8,19]  computed to one-loop accuracy. Here the Z-factors become momentum-fraction dependent (contrary to covariant gauges) due to assumed principal value prescription for the 1/k + -pole in the gluon density matrix, where and for quark and gluon, respectively. Their contribution to the renormalization of the operator blocks reads where Γ µ 1 ...µn is the Dirac-Lorentz tensor defining the composite operator in question. The collinearly divergent integrals entering Σ's and Π's regulate the end-point singularities in the momentum-fraction kernels promoting them to conventional plus-distributions that become integrable over the entire range of momentum fractions [21].

Two-to-two transitions: quasipartonic operators
To make our expressions more compact, we introduce a set color structures with open indices that show up in our expressions where f abc is the SU(N) structure constants while t a andt a are the SU(N) and SU(N) generators in the fundamental representation and its conjugate, respectively.
In this quark-quark sector the fields have open fundamental color indices i 1 and i 2 . The operator renormalization kernel acts on them as follows and its explicit expression arises from the graph shown in Fig. 1. It is given by For the nonsinglet sector, the Feynman diagram responsible for the evolution is determined by the very same one-gluon exchange in Fig. 1 where For the quark-antiquark operators of the same flavor Φ }, there are two extra transitions corresponding to annihilation channels. Although we do not focus on the flavor singlet quark operators and the operators built up solely by gluon fields, we do provide corresponding results for the 2 → 2 evolution kernels in Appendix B.
For the quark-gluon operator blocks, the renormalization opens up more than one color channel, with corresponding transition kernels calculated from the graphs shown in Fig. 2 being x 2 ) Similar results are obtained by replacing the quark and an antiquark field, with The Feynman graphs involved in the analysis differ from Fig. 2 only by the orientation of one of the quark lines. All of these expressions agree with well-known particle transitions for quasipartonic operators [5,8,6,10,12].

Two-to-two transitions: non-quasipartonic operators
Now, we turn to the analysis of non-quasipartonic operators. According to the adopted basis (25), the new two-particle blocks that we have to address contain a bad field component accompanied by a good one, namely

Quark-quark transitions
To start with, we consider the quark-quark transitions first. To this end we introduce nonquasipartonic two-particle operator built up from primary fields and arranged as doublets since they mix under renormalization group evolution, The Feynman diagram responsible for the mixing addressed in this section is the same one as in Fig. 1. We elaborate on an example of a specific transition in great detail in Appendix A to demonstrate the inner workings of the formalism. As a result, we find where the two-by-two mixing matrix possesses the elements For O ij − operator sets, we similarly get where This concludes our discussion on non-singlet operators. In Appendix B.2.1, we also provide transitions into gluonic operators involved in this class when it is generalized to the singlet channel as well.

Quark-antiquark transitions
Next, we introduce doublets of quark-antiquark fields where we assume that the two-particle blocks possess different flavor such that they do not undergo annihilation transitions into gluon fields. Then the evolution equation can be written as with the elements of the evolution matrix given by

Quark-gluon transitions
For the operators involving one quark and one gluon field, we introduce the following two-vectors Feynman diagrams corresponding to two-to-two transition of non-quasipartonic quarkgluon blocks in Sect. (4.2.3).
Then, calculating Feynman diagrams responsible for their one-loop renormalization demonstrated in Fig. 3, we deduce that as in the quasipartonic case, there are two color-flow channels that induce the transitions where the elements of the kernels K ij and K ij admit the form Similarly for the operators in the O − group, we get with Having studied the operators generated by primary fields with the same chiralities, we now turn our attention to the cases where the operators are built up by fields of opposite chiralities, namely, Their one-loop evolution equation is driven by and the transition kernels involved read

Antiquark-gluon transitions
For antiquark-gluon blocks, we introduce the doublets whose transitions are determined by computing Feynman diagrams in Fig. 4 and yield x 2 1 y 2 2 6y 2 1 + 6y 1 y 2 + y 2 2 − 4y 1 y 2 2 x 2 + y 1 x 2 2 (y 1 + 2y 2 ) This concludes our analysis of non-singlet two-to-two transitions of nonquasipartonic operators. They agree with corresponding findings in [8], the last paper of Ref. [6] and [12] after the Fourier transformations. For a partial result in the singlet channel, we refer the reader to Appendix B.

Two-to-three transitions: non-quasipartonic operators
Finally, we come to the analysis of particle number-changing transitions. This is the most elaborate sector of twist-four operators. Apart from proliferation of Feynman graphs, there are also subtle effects related to transitions induced by the use of QCD equations of motion. We provide for the latter a diagrammatic representation that puts it on the same footing as the rest of the calculation and thus reduces the procedure to tedious algebraic manipulations. An example exhibiting the formalism is worked out in Appendix D.  4.3.1 ψ − ψ + and 1 2 D −+ψ+ψ+ Let us start our study with the quark-quark bad-good and good-good operators with a transverse derivative, respectively, mixing with the following three-particle operator constructed from good field components In both cases, there are three nontrivial color-flow channels such that the evolution equation takes the form First, for the O ij (z 1 , z 2 ) = ψ i − (z 1 )ψ j + (z 2 ) case, the evolution kernels computed from the diagrams in Fig. 5 read Making use of the differential operators introduced in section 3.3, it is straightforward to verify that these kernels are all conformally invariant.

f
undergo a transition to the following three-field operator according to where the decomposition runs over the color structures given in Eq. (158). The evolution kernels for the two cases O ai (x 1 , 195) and respectively. The first set comes from Fig. 10, while the second one from both Figs. 10 and 11. This time we have found an exact agreement with the findings of Ref. [12] without further implementation of symmetry properties. This can be easily explained by the fact that the operators in this case lack the "field exchanging" symmetries. As a result, no redundancies in the evolution kernels are allowed to be left over. By taking the heavy quark limit, we also reproduced the results reported in Ref. [22].
transition to It is described by with the color structures given in Eq. (158).

Outlook and Conclusion
In this paper, we generalized the formalism suggested by Bukhvostov-Frolov-Lipatov-Kuraev for renormalization of quasipartonic operators to include nonquasipartonic operators as well. The advantage of the method is that at one-loop order, the procedure is purely algebraic requiring straightforward though quite tedious manipulations with Dirac and Lorentz structure of Feynman graphs. The focus of the present study was the evolution equations for non-singlet twist-four operators. Their basis consists of four-particle quasipartonic and three-particle good-good-bad light-cone operators. While the former were studied at length in existing literature, the latter were addressed here starting from Feynman graphs, providing an explicit brute force calculation of these evolution kernels. The main ingredients for these transitions are good-bad two-to-two and two-to-three components. A crucial role in both cases is placed by proper use of QCD equations of motion which induce extra contribution that are required for proper closure of evolutions equations. Since the basis of twist-four operators is built from conformal primary fields, the resulting evolution kernels should obey a very stringent consistency constraint of being conformally invariant. This was explicitly confirmed by our analysis. We provided a Fourier transform from the momentum to coordinate space and back and checked our findings against the only available earlier results for nonquasipartonic operators that were derived for light-ray operators making use of the conformal symmetry and dynamical part of the Poincaré algebra. We found agreement in all cases and also provided a simplified form of light-ray kernels in certain channels that made use of the exchange symmetry of the operators involved.

Acknowledgments
We would like to thank A. Manashov for very instructive discussions and clarification with regards to results of Ref. [12] and N. Offen for helpful correspondence about heavy light-ray operators and, last but not least, M. Ramsey-Musolf and M. Glatzmaier for their interest in the project at its early stages. This work was supported by the U.S. National Science Foundation under the grant PHY-1068286.

A Sample calculations in light cone gauge
In this appendix we provide an explicit calculation of the transitions kernel for good-bad two-totwo quark transitions χ + ⊗ ψ − → χ + ⊗ ψ − + χ − ⊗ ψ + , shown in Fig. 12. The operator in question can be written at one loop in the form where in the gluon propagator in the light-cone gauge was introduced in Eq. (34), while for reader's convenience we provide below expressions for the quark propagator and and the vertex function, respectively, Denoting the string introduced in curly brackets as N /D, we can work out the denominator D stemming from the propagators as D = (p 1 + p 2 − k 1 ) 2 (p 1 − k 1 ) 2 k 2 1 . Choosing the loop momentum as k = k 1 , we expand D in inverse powers of the transverse momentum k ⊥ and find immediately for the leading and first subleasing contributions We will parametrize the contributions of the first, second and third terms in the square brackets as A, B and C contributions, respectively, i.e., A − B − C.
To clarity the manipulations involved in the analysis, the numerator will be calculated term by term. To start with, notice that p − 1 and p − 2 can be automatically neglected in the calculation as they vanish for Fourier transform of light-ray operators that we consider in this work. Let us start with the g µν piece and denote its contraction with the strong of Dirac matrices in Eq. (228) as I 1 . Then after Sudakov decomposition of all momenta and little Dirac algebra, we find after rescaling the k − momentum component according to Eq. (42) where we have neglected all terms that do not produce any divergences, i.e., terms scaling as k n ⊥ with n < 2. Next, we turn to the second (k − p 1 ) µ n ν and third (k − p 1 ) ν n µ terms. For their contraction with the scare bracket, we find in a completely analogous manner Now, combining the above in the integrand, we trace only terms with 1/k 2 ⊥ behavior since these are the only contributions yielding logarithmic divergence. Integrating over the longitudinal k + component with the help of the Dirac delta function in Eq. (226) where in the last step we restored the omitted causal i0 prescription in the longitudinal denominators use the defining integral Eq. (43) for the generalized step functions. Similarly, we find To proceed further with other contributions, we compute the following integral first 2π 0 dϕ 2π Here we used the fact that the integrand does not have any vectors but k ⊥ so that one can immediately calculate the average in the two-dimensional transverse plane k α ⊥ k β ⊥ → k 2 ⊥ δ αβ /2. Thus we obtain Putting all the pieces together, we get Finally, using equations of motion for the (anti)quark fields, with neglected gluon field since we after the two-to-two transitions only, (p + 2 γ − + / p 2⊥ )ψ(p 2 ) = 0 andψ(p 1 )(p + 1 γ − + / p 1⊥ ) = 0, we can trade transverse momenta accompanying the good component of the quark to the bad quark fields. This way we arrive at Eqs. (89)-(92).

B Flavor singlet 2 → 2 transitions
We complement the non-singlet analysis performed in the body of the paper with particle results involving the singlet sector. In all cases we found agreement with corresponding expressions reported in Ref. [12].

B.1 Quasi-partonic operators
To start with, we present the quasipartonic quark-antiquark to gluon-gluon kernels and gluongluon transitions as well. Figure 13: Two-to-two quark-gluon transitions in Eq. (239). The color structures C c are defined in Eq. (73).
In the singlet sector, the quark-antiquark evolution will produce extra annihilation-type contributions shown by the first two graph in Fig. 13 [ in addition to already computed transitions, denoted above by ellipses, and given in Eqs. (76) and (77). The index f runs over all quark flavors. The last line represents transitions into gluons, exhibited by the last two graphs in Fig. (13). The color structures are displayed in Eq. (73). The transition kernels then read In the following two subsection, we will list our results of the evolution kernels for the pure gluonic transitions.
++ }(x 1 , x 2 ) For gluon blocks of the same chirality, the nonvanshing Feynman graphs that induce the transition are given in Fig. 14 and produce Here we again observe the "exchange symmetry" elaborated in details in Sects. 4.3.2 and 4.3.3. In the present case, it implies the simultaneous interchange of a ↔ b and z 1 ↔ z 2 .
Finally, the opposite-chirality gluon sector evolves as according to nontrivial Feynman diagrams in Fig. 15 with All other quasipartonic singlet transitions can be found in the literature [5,8,10,12]

B.2 Non-Quasipartonic operators
In this Appendix, we complement non-quasipartonic operators with purely gluonic transitions, thus extending the consideration of Sect. 4.2.

D Equation-of-motion graphs
In the preceding Appendix A illustrating the good-bad two-to-two transitions, we already had to rely on the use of equations of motion to produce correct evolution kernels. We ignored however the effects of the gluon field. In the analysis of the two-to-three transitions, we have to restore this neglected contributions. Here we demonstrate how this can be achieved. We will choose a nontrivial diagram from Sect. 4.3.3 to illustrate our point. For simplicity we only consider quark fields. For gluons, similar logic is applicable in a straightforward fashion though the algebra gets a bit more involved. We start with the equation for a massless quark field, / Dψ(z) = / ∂ − ig / A(z) ψ(z) = 0 (273) that translates in momentum space to / pψ(p) = −g d 4 p ′ / A(p ′ )ψ(p − p ′ ) .
As before, the fields' momenta do not possess the "−" components. A way to extract the gluon splitting off the quark line due to the equation of motion is to collect terms proportional to / p. In practice, however, this proves to be extremely difficult. The trick instead is to explicitly spill a gluon off the quark line with ordinary Feynman rules, ψ(p) → ψ(p 1 )A(p 2 ) as shown in the graph The emerging quark propagator P(p 1 + p 2 ) has a pole as it goes on-shell due to the collinearity of emitted gluon off a collinear quark. However, we can introduce nonvanishing transverse momenta for outgoing quarks and gluons and then collect the terms of the form (p 1 + p 2 ) 2 = (p 1 + p 2 ) 2 ⊥ to cancel the singular denominator in P(p 1 + p 2 ).
As an explicit example, let us consider the diagram + (x 2 1 + (x 2 − y 1 )y 1 + x 1 (x 2 + y 1 )ϑ 1 112 (x 1 , x 1 − y 1 , −x 2 )) + (x 1 + y 1 )ϑ 1 111 (x 1 , where we have restored the color structures and Γ = γ + (1 − iγ 1 γ 2 ) is the matrix structure corresponding to the projection onto operatorf a ++ (y 1 )ψ + (y 2 )f d ++ (y 3 ). Here only terms proportional to q 2 ⊥ are kept since they are the ones corresponding to the equation of motion operators that we are dealing with. We should point out that our management of the equation-of-motion diagrams coincide with the procedures described in Ref. [23].
E Light-ray kernels in section 4.3.2 and 4.3.3 Here we convert our diagrammatic results given in Sects. 4.3.2 and 4.3.3 into the coordinate space. This is done by inverse-Fourier transform the corresponding kernels in momentum space making use of formulas in Sect. 3.2 and applied in Appendix C.