Constraints for twist-two alien operators in QCD
G. Falcioni F. Herzog S. Moch S. Van Thurenhout













T> c
M>c <23.0cm  18.5cm
 -0.9cm  -0.9cm
 -1.5cm   4colour#1#1g_s(#1)α_s^ #1g_s^2equationsection
ZU-TH 43/24  August 2024 

DESY-24-108 



 Constraints for twist-two alien operators in QCD




G. Falcioni^  a,b,
F. Herzog^  c,
S. Moch^  d and
S. Van Thurenhout^  e


^a Dipartimento di Fisica, Università di Torino, Via Pietro Giuria 1, 10125 Torino, Italy

^b Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland

^c Higgs Centre for Theoretical Physics, School of Physics and Astronomy, 


The University of Edinburgh, Edinburgh EH9 3FD, Scotland, UK

^dII. Institute for Theoretical Physics, Hamburg University


Luruper Chaussee 149, D-22761 Hamburg, Germany

^e HUN-REN Wigner Research Centre for Physics, Konkoly-Thege Miklós u. 29-33, 1121
Budapest, Hungary



 Abstract




Parton evolution equations in QCD are controlled by the anomalous dimensions of gauge-invariant twist-two spin-N quark and gluon operators.
Under renormalization, these mix with gauge-variant operators of the same quantum numbers, referred to as alien operators.
Our work addresses the systematic study of these alien operators at arbitrary spin N, using generalized BRST symmetry relations to derive their couplings and Feynman rules at all values of N.
We observe how the all-N structure of the generalized (anti-)BRST constraints relates the couplings of alien operators with n+1 gluons to those with n gluons. 
Realizing a bootstrap, we present all one-loop results necessary for performing the operator renormalization up to four loops in QCD.



§ INTRODUCTION


The study of twist-two operators of spin-N for quarks and gluons in quantum chromodynamics (QCD) and their renormalization dates to the origins of QCD as the gauge theory of the strong interaction <cit.>. 
The renormalization of off-shell operator matrix elements (OMEs) in QCD, i.e. Green's functions with off-shell external momenta and insertions of these quark and gluon operators, gives access to their anomalous dimensions.
These coincide with the Mellin transforms of the standard QCD splitting functions, that govern the scale evolution of the parton distribution functions.
It is well-known that the twist-two operators of spin-N mix under renormalization
with a set of gauge-variant operators of the same quantum numbers, which
involve equation-of-motion (EOM) and ghost operators.
The latter, often referred to in summary as alien operators, can be constructed systematically, by employing a generalized gauge symmetry of the QCD Lagrangian in covariant gauge with the addition of the physical quark and gluon operators <cit.>.
The generalized gauge symmetry can be promoted to a generalized BRST (gBRST) symmetry <cit.>. 
This provides an algebraic approach for the derivation of a complete set of operators to be considered in the renormalization of the off-shell OMEs at a given loop order in perturbative QCD in an expansion in the strong coupling g_s, α_s=g_s^2/(4π). The complete set of operators required up to four loops has been listed in <cit.>.

Each alien operator features a coupling constant that can be interpreted as the renormalization constant that generates mixing of the gauge-invariant operators into each alien. In order to renormalize the physical OMEs, these coupling constants must be computed order-by-order in perturbation theory. The required couplings to renormalize the two-loop OMEs were computed in <cit.> in closed form for all values of N. A method to determine the alien counterterms, i.e. the Feynman rules obtained by summing all the alien operators with their associated couplings, was presented in <cit.> together with results up to the three-loop level for a covariant gauge and all values of N. From this, the n_f^2 contributions to the pure-singlet splitting functions were obtained at four loops <cit.>. Beyond three loops, ref. <cit.> determined a set of all-order constraints on the couplings, induced by gBRST and generalized anti-BRST symmetries <cit.>. In <cit.>, these constraints were solved at arbitrary loop order for fixed N≤20, leaving the systematic study of the alien operators at arbitrary spin N as an open problem. In this paper, we follow a different strategy. Namely, we will solve the constraints on the alien couplings to leading order in g_s but for all values of N. The main results of our study are:


  *  The all-N structure of the couplings is fixed in terms of a small set of constants. The latter can be determined by explicitly computing the couplings for some fixed values of N.

  *  The structure of the couplings of alien operators with n+1 gluons is related to the ones with n gluons, allowing for a bootstrap in the determination of complicated higher-order couplings in terms of simpler lower-order ones. 

The outline of the article is as follows. 
In Sec. <ref> we set the stage, review the generalized gauge symmetry and provide a brief summary of the set of relevant alien operators. In Sec. <ref> we study the identities that exist among the couplings of the alien operators and show how they can be used to restrict the all-N structure of the couplings. The results of this analysis are then used in Sec. <ref> to derive the Feynman rules of the alien operators, suitable for the renormalization of OMEs at all N up to four loops in QCD. Finally, in Sec. <ref>, we summarize our findings and provide an outlook on further developments.



§ SETTING THE STAGE

In this section, we review the construction of the alien operators and summarize our conventions. The complete gauge-fixed QCD action is written as

    S=
        ∫d^Dx (ℒ_0+ℒ_GF+G) .

Here ℒ_0 represents the classical part of the QCD Lagrangian

    ℒ_0 = -1/4 F^μν_a F_μν^a + ∑_f=1^n_fψ^f(iD-m_f)ψ^f ,

with the field strength defined as

    F_μν^a = ∂_μ A_ν^a - ∂_ν A_μ^a + g_s f^abc A_μ^bA_ν^c .
f^abc are the standard QCD structure constants. The covariant derivative in Eq. (<ref>) is D= γ^μ(∂_μ-ig_s T^aA^a_μ) with T^a the generator of the gauge group in the fundamental representation. The gauge-fixing and ghost terms are

    ℒ_GF+G=-1/2ξ(∂^μ A^a_μ)^2-c^a ∂^μ D_μ^ab c^b

with ξ the covariant gauge parameter and c^a and c^a the anti-ghost and ghost fields, respectively. The covariant derivative in the adjoint representation is D_μ^ac=∂_μδ^ac+g_s f^abcA_μ^b. The QCD Lagrangian can be extended to also include spin-N gauge-invariant operators of twist two, which we define as

    _ g^(N)(x)    =   1/2Tr[F_ν(x) D^N-2 F^ν(x)] ,
    _ q^(N)(x)    =   1/2Tr[ψ(x)Δ D^N-1ψ(x)]  .

Here Δ_μ is a lightlike vector and we introduced the notation

    F^μ;a=Δ_ν F^μν;a,      A^a=Δ_μ A^μ;a,     D=Δ_μ D^μ,    ∂=Δ_μ∂^μ .

Under renormalization the operators in Eq. (<ref>) mix with operators proportional to the (classical) EOM and with BRST-exact operators <cit.>. Following <cit.>, we begin by presenting the EOM aliens in the form

    _EOM^(N)=(D· F^a  +g_sψT^aΔψ) 𝒢^a(A^a,∂ A^a,∂^2  A^a,... )

with D· F^a = D_ν F^ν;a and 𝒢^a a generic local function of the gauge field and its derivatives. It is convenient to expand 𝒢^a in a series of contributions with an increasing number of gauge fields. This leads to

    _^(N) = _^(N),I+_^(N),II+_^(N),III+_^(N),IV+ …

with


    _EOM^(N),I   = 
        η(N) (D· F^a + g_s ψ Δ T^a 
        ψ) (∂^ N-2A^a ), 
    _EOM^(N),II   = 
     (D· F^a + ψ Δ 
       T^a ψ) ∑_i+j
    =N-3C_ij^abc(∂^iA^b)(∂^jA^c),
    _EOM^(N),III   = 
       (D· F^a + ψΔ T^aψ) ∑_i+j+k
    =N-4C_ijk^abcd(∂^iA^b)(∂^jA^c)(∂^kA^d),
    _EOM^(N),IV   = 
      g_s^3  (D· F^a + ψΔ T^aψ)∑_i+j+k+l
    =N-5C_ijkl^abcde(∂^iA^b)(∂^jA^c)(∂^kA^d)(∂^lA^e).

The coefficients C^a_1… a_n_i_1… i_n-1 appearing in Eqs. (<ref>)-(<ref>) can be written in terms of a set of independent colour tensors, each of them multiplying an associated coupling constant, as follows

    C_ij^abc = f^abcκ_ij, 
       C_ijk^abcd=(f  f)^abcdκ_ijk^(1)+d_4^abcdκ_ijk^(2)+d_4ff^abcdκ_ijk^(3), 
       C_ijkl^abcde = (f f f)^abcdeκ_ijkl^(1)+d_4f^abcdeκ_ijkl^(2)

with

    (f  f)^abcd = f^abef^cde, 
       (f f f)^abcde = f^abmf^mcnf^nde, 
       d_4^abcd = 1/4![Tr(T_A^aT_A^bT_A^cT_A^d)+symmetric permutations], 
       d_4ff^abcd = d_4^abmnf^mcef^edn, 
       d_4ff^abcd = d_4ff^abcd-1/3C_A d_4^abcd, 
       d_4f^abcde = d_4^abcmf^mde.

Here (T_A)^b_ac = i f^abc are the generators of the adjoint representation of the colour group. We now extend the classical Lagrangian ℒ_0 in Eq. (<ref>) to include the gauge-invariant operators of twist two as well as the EOM aliens

    ℒ_GGI = ℒ_0+w_ i _ i^(N)+_EOM^(N),

where w_ i is a coupling for the operator _ i with i=g,q, playing the same role as the coefficients η(N), κ_ij, … defined in Eqs. (<ref>)-(<ref>). The Lagrangian ℒ_GGI is invariant under the generalized gauge transformation <cit.>A^a_μ→ A^a_μ + δ_ω A^a_μ + δ_ω^Δ A^a_μ, where 

    δ_ω A^a_μ   = D^ab_μω^b(x),
    δ_ω^Δ A^a_μ   = -Δ_μ[η(N) ∂^N-1ω^a + g_s∑_i+j
    =N-3C^aa_1a_2_ij (∂^iA^a_1) (∂^j+1ω^a_2)
       +g_s^2∑_i+j+k
    =N-4C^aa_1a_2a_3_ijk (∂^iA^a_1) (∂^jA^a_2) (∂^k+1ω^a_3)
       +g_s^3 ∑_i+j+k+l
    =N-5C^aa_1a_2a_3a_4_ijkl (∂^iA^a_1) (∂^jA^a_2) (∂^kA^a_3) (∂^l+1ω^a_4) + O(g_s^4)]

and 

    C_ij^abc = f^abcη_ij,  
       C_ijk^abcd=(f   f)^abcdη_ijk^(1)+d_4^abcdη_ijk^(2)+d_4ff^abcdη_ijk^(3), 
       C_ijkl^abcde = (f f f)^abcdeη_ijkl^(1)+d_4f^abcdeη_ijkl^(2a)+d_4f^aebcdη_ijkl^(2b).

The generalized gauge symmetry implies that the couplings η^(k)_n_1… n_j are related to κ^(k)_n_1… n_j in Eqs. (<ref>)-(<ref>)

    η_ij   =2κ_ij+η(N)i+j+1i,
    η^(1)_ijk   =2κ_i(j+k+1)j+k+1j+2[κ_ijk^(1)+κ_kji^(1)],
    η^(2)_ijk   =3κ^(2)_ijk,
    η^(3)_ijk   =2[κ^(3)_ijk-κ^(3)_kji],
    η^(1)_ijkl   =2[κ_ij(l+k+1)^(1)+κ_(l+k+1)ji^(1)]l+k+1k+2[κ_ijkl^(1)+κ_ilkj^(1)+κ_likj^(1)+κ_lkij^(1)],
    η^(2a)_ijkl   =3κ_ij(k+l+1)^(2)k+l+1k+2κ_ijkl^(2),
    η^(2b)_ijkl   =2κ^(2)_lijk.

The new gauge transformations in Eq. (<ref>) are promoted to a nilpotent generalized BRST (gBRST) operator, by replacing the transformation parameter ω^a with the ghost field c^a<cit.>. In turn the ghost alien operator is generated by the action of such gBRST operator on a suitable ancestor operator <cit.>, giving

    _c^(N) = _c^(N),I+_c^(N),II+_c^(N),III+_c^(N),IV+ …

with

    _c^(N),I = -η(N) (∂c^a)(∂^N-1c^a), 
       _c^(N),II = -g_s ∑_i+j
    =N-3C_ij^abc(∂c^a)(∂^iA^b)(∂^j+1c^c), 
       _c^(N),III = -g_s^2∑_i+j+k
    =N-4C_ijk^astu(∂c^a)(∂^iA^s)(∂^jA^t)(∂^k+1c^u), 
       _c^(N),IV = -g_s^3∑_i+j+k+l
    =N-5C_ijkl^abcde(∂c^a)(∂^iA^b)(∂^jA^c)(∂^kA^d)(∂^l+1c^e) .


  
Renormalization The complete Lagrangian, including the twist-two physical and alien operators, can be written as


    ℒ   = ℒ_0 + ℒ_GF+G + w_ i _ i + _EOM^(N) + _c^(N) =ℒ_0(A^a_μ,g_s) + ℒ_GF+G(A^a_μ,c^a,c̅^a,g_s,ξ) + ∑_ k 𝒞_ k _ k,

where 𝒞_ k labels all the distinct couplings of the operators, e.g. 𝒞_ k=(w_ i,η(N), κ_0 1, κ_1 2…). The ultraviolet (UV) singularities associated with the QCD Lagrangian are absorbed by introducing the bare fields and parameters

    A^a;bare_μ(x) = √(Z_3) A^a_μ(x),       c^a;bare(x)=√(Z_c) c^a(x),       c̅^a;bare(x)=√(Z_c) c̅^a(x),
          g_s^bare=μ^ϵ Z_g g_s,       ξ^bare=√(Z_3) ξ.

We renormalize the singularities originating from the insertion of the composite operators using


    _ i^ren(x)    = Z_ ij _ j^bare(x),

where _ j^bare indicates the operators in Eqs. (<ref>), (<ref>) and (<ref>) written in terms of the bare fields. Note that throughout this work we use D=4-2 dimensional regularization, combined with the  renormalization scheme. Z_ ij is the renormalization matrix of the operators, which makes the OMEs featuring an insertion of _ i^ren finite. The renormalized Lagrangian becomes


    ℒ = ℒ_0(A^a;bare_μ,g_s^bare) + ℒ_GF+G(A^a;bare_μ,c^a;bare,c̅^a;bare,g_s^bare,ξ^bare) + ∑_ k𝒞_ k^bare _ k^bare,
       𝒞_ i^bare = ∑_ k𝒞_ k Z_ k i,

where 𝒞_ k is the (finite) renormalized coupling of the operator _ k. The UV-finite OMEs featuring a single insertion of _ g^ren are computed by setting the renormalized couplings 𝒞_ i=δ_ i g in Eq. (<ref>), which gives


    𝒞_ i^bare=Z_ g i.

Similarly, the renormalized OMEs with an insertion of _ q are obtained with 𝒞_ i^bare=Z_ q i. Therefore, the couplings of the bare operators η^bare(N), … are interpreted as the renormalization constants that mix the physical operators into the aliens. These quantities can be extracted from the direct calculation of the singularities of the OMEs with an insertion of _ g^bare (_ q^bare). For instance, the coupling η^bare(N), which is associated to an operator with a two-point vertex, was determined in <cit.> from the renormalization of the OMEs of _ g with two external ghosts and it was found to be [Note that the expression for η in <cit.> has an additional factor of 2. This is a consequence of the chosen conventions for dimensional regularization. In particular, we use D=4-2 while <cit.> employs D=4+.]

    η^bare(N)=Z_g c = -a_s/ϵC_A/N(N-1)+O(a_s^2),

where C_A is the quadratic Casimir in the adjoint representation and 
a_s=α_s/(4π)=g_s^2/(4π)^2. 
The value of η was determined at two loops in <cit.> and at three loops in <cit.>. Throughout this paper we will mainly be interested in the one-loop alien couplings. As such, it will be convenient to select just the N-dependent part of the one-loop result of η^bare(N), which in the following we simply denote by η(N), i.e.

    η(N) = -1/N(N-1).

The couplings of the operators featuring multiple fields, e.g., the couplings κ_ij multiply at least three fields, are determined by renormalizing OMEs with the corresponding external fields. Recently, a method to compute the counterterms of the OMEs with insertions of the gauge-invariant operators as a function of the spin N was put forward in ref. <cit.>. The result of that paper can be used to extract the coefficients κ^bare_i j up to O(a_s^2) and those of κ^(p);bare_ijk, for p=1,2, at O(a_s), finding agreement with the low-N values reported in <cit.>. In addition, the calculation of the five-point counterterms at O(a_s), which can be used to determine κ^(p)_ijkl, for p=1,2, has been announced recently <cit.>. 

In this paper, we would like to determine the renormalization constants Z_ gi by solving the constraints on the couplings 𝒞_ i^bare, which are imposed by the symmetries of Eq. (<ref>). The latter is the Lagrangian in Eq. (<ref>) evaluated with bare fields and couplings constants. Therefore the two Lagrangians share the same symmetry properties, with the obvious substitutions. 
For simplicity, in the rest of this paper we drop the superscript `bare', wherever it does not create any ambiguity. 



  
Independent operators and couplings
The symmetry constraints on the couplings in Eq. (<ref>) have been derived in ref. <cit.>. Without repeating the derivation of that paper, we distinguish three types of relations, which follow from the way we have constructed the operators at the beginning of this section.

First of all, the couplings introduced in the EOM operators, see Eqs. (<ref>)-(<ref>) and (<ref>)-(<ref>), are chosen to inherit the properties of the colour structures they multiply. For example, because of the anti-symmetry of the structure constants, we take

    κ_ij=-κ_ji.

This implies, e.g., that at spin N=4, where i,j=0,1, there is only one independent coupling, e.g., κ_0 1.

The second type of constraints regards the couplings that enter the ghost operators, Eqs. (<ref>)-(<ref>). Because these operators were constructed directly from the EOM ones using gBRST, the η couplings are connected to the κ ones. The relevant identities have been listed in Eqs. (<ref>)-(<ref>). 

Finally, we impose the invariance of Eq. (<ref>) under the generalized transformations of anti-BRST type <cit.>, which stem from Eq. (<ref>), by replacing the transformation parameter ω^a with the anti-ghost field c̅^a. This implies the following condition on the ghost operator _c^(N) defined in Eq. (<ref>)

    O_ c^(N)(A^a_μ,c^a,c̅^a) = O_ c^(N)(A^a_μ,c̅^a,c^a),

which translates into a set of constraints on the couplings in Eqs. (<ref>)-(<ref>) and, in turn, on those of the EOM operators. Taking the example of N=4, the anti-BRST relation imposes κ_01=2η(4), thus reducing the number of independent couplings even further <cit.>.

It is highly non-trivial to find all-N solutions for all the constraints. In refs. <cit.>, they were solved only for fixed values of N, in order to fix bases of independent alien operators up to N=20. In the following sections, we solve the relations with exact N dependence. This is done by setting up an ansatz for the function space that enters to leading order in a_s. The construction of this ansatz is primarily based on constraints from (anti-)gBRST. We will see below that the latter allow one to bootstrap the functional form of higher-order couplings from that of the lower-order ones. The determination of the unknown parameters in the ansatz is then performed by using the full set of colour, gBRST and anti-gBRST relations. As will become clear below, this allows one to fix most, but not all, free parameters. The few that remain then need to be determined from the explicit renormalization of a limited number of fixed-N operator matrix elements. This is particularly important for finding any overall N-dependent function.



§ IDENTITIES AMONG THE ALIEN COUPLINGS


In this section we will discuss in detail the identities between the couplings coming from the (anti)-gBRST relations. In particular, we will show that they allow one to restrict the function space of the couplings and hence constrain their generic N-dependence.



 §.§ Class II couplings

The class II operators are defined in terms of two couplings, κ_ij and η_ij, which obey the following relations

    κ_ij+κ_ji=0,    [anti-symmetry of f]
       η_ij=2κ_ij+η(N)i+j+1i,    [gBRST]
       η_ij+∑_s=0^i(-1)^s+js+jjη_(i-s)(j+s) = 0.    [anti-gBRST]

Note that one can generate an equation for the ghost coupling alone by combining the anti-symmetry of κ_ij, Eq. (<ref>), with the gBRST relation, Eq. (<ref>),

    η_ij+η_ji = η(N)[i+j+1i+i+j+1j].

The one-loop value of this coupling was first computed in <cit.>
and later corrected in  <cit.>.
In our conventions it reads [Note that there are typos in the corresponding expression in <cit.>. In particular, the right-hand side of Eq. (4.38) in <cit.> should be replaced by the right-hand side of Eq. (<ref>) here.
]

    η_ij  =-η(N)/4[(-1)^j-3N-2i+1-N-2i]

which implies

    κ_ij = -η(N)/8[(-1)^j+3N-2i-3N-2i+1].

The power of the relations described above is that they can be used to gain valuable information about the structure of the couplings at arbitrary N. For example, one can use Eq. (<ref>) to write down an ansatz for η_ij of the form

    η_ij = η(N) [ c_1i+j+1i+c_2i+j+1j].

Here c_1 and c_2 are constants to be determined.  We assume here that the dependence on η(N) is factorized at leading order, 
as suggested by Eq. (<ref>) and observed in Eq. (<ref>). 

This ansatz can then be substituted in the anti-gBRST consistency relation, Eq. (<ref>), yielding

    η_ij+∑_s=0^i(-1)^s+js+jjη_(i-s)(j+s) =η(N) [ (-1)^jc_1-c_2i+j+1j]

for even values of N. Hence, only the trivial solution c_1=c_2=0 obeys the anti-gBRST relation. However, the right-hand side of Eq. (<ref>) suggests the inclusion of a term proportional to (-1)^j to the ansatz,

    η_ij =η(N) [ c_1i+j+1i+c_2i+j+1j+c_3(-1)^j].

The anti-gBRST relation now becomes

    η_ij+∑_s=0^i(-1)^s+js+jjη_(i-s)(j+s) = η(N) (c_1+c_3)[(-1)^j+i+j+1i]

such that c_3=-c_1 is a consistent solution. If we now impose also Eq. (<ref>) then we obtain the relation c_1+c_2=-4, leaving just one free parameter unconstrained. 
Hence

    η_ij =η(N) [ c_1[i+j+1i-(-1)^j]-(4+c_1)i+j+1j].


It should be noted that, if an ansatz is generated using (anti-)gBRST relations, one is in principle free to add non-zero functions that live in the kernel of these relations. For example, if one adds a term of the form

    -f(N)/4((-1)^j + N - 2i + 1 - N - 2i)

to Eq. (<ref>), the corresponding expression for η_ij still obeys the constraints. Here f(N) represents an arbitrary function of N, with the actual solution being recovered by setting f(N)=0 for even values of N. In particular, substituting Eq. (<ref>) in the constraint coming from anti-symmetry and gBRST, cf. Eq. (<ref>), one finds

    [(-1)^i+(-1)^j ]f(N)=0.

The left-hand side of this expression always vanishes for all physical (even) values of N, independent of the functional form of f(N). In general, the exclusion of this type of function can only be confirmed by comparison with fixed-N computations.

An important consequence is that now we have recovered the full function space of the actual solution, Eq. (<ref>), using only the symmetry relations of the couplings. More generally, note that Eq. (<ref>) is an example of a conjugation relation, in the sense that a second application of the sum leads to

    ∑_t=0^i(-1)^t+jt+jjη_(i-t)(j+t) = -∑_t=0^i(-1)^t+jt+jj∑_s=0^i-t(-1)^s+j+ts+j+tj+tη_(i-t-s)(j+t+s)

and hence

    η_ij = ∑_t=0^it+jj∑_s=0^i-t(-1)^ss+j+tj+tη_(i-t-s)(j+t+s).

The latter identity is actually always true for any discrete two-variable function η_ij. This type of conjugation relation has already been encountered in the computation of the anomalous dimensions of twist-two operators in non-forward kinematics, see e.g. <cit.>, and holds great predictive power. In particular, it provides valuable information about the function space of the object at hand. To take full advantage of such relations, one needs to be able to evaluate them analytically. This is possible by using principles of symbolic summation, in particular by application of the creative telescoping algorithm <cit.>. The latter is a generalization of classical telescoping and attempts to evaluate the sum of interest by rewriting it as a recursion relation using Gosper's algorithm <cit.>. The closed-form expression of the sum then corresponds to the linear combination of the solutions of the recursion that has the same initial values as the sum. This methodology is neatly implemented in the  Mathematica package  Sigma<cit.>. For the class III and IV couplings to be described below we will also encounter identities involving multiple sums, for which the package  EvaluateMultiSums<cit.> can be used.





 §.§ Class III couplings


  §.§.§ κ_ijk^(1) and η_ijk^(1)

The couplings η_ijk^(1) and κ_ijk^(1) can be thought of as direct generalizations of η_ij and κ_ij in the class II operators. They obey the following relations

    κ_ijk^(1)+κ_ikj^(1)=0,    [anti-symmetry of f]
       κ_ijk^(1)+κ_jki^(1)+κ_kij^(1) = 0,    [Jacobi identity]
       η_ijk^(1)=2κ_i(j+k+1)j+k+1j+2[κ_ijk^(1)+κ_kji^(1)],    [gBRST]
       η_ijk^(1)=∑_m=0^i∑_n=0^j(m+n+k)!/m!n!k!(-1)^m+n+kη_(j-n)(i-m)(k+m+n)^(1).    [anti-gBRST]

Note that now the indices are constrained as i+j+k=N-4. As before, one can combine the relations of the EOM coupling with the gBRST relation to connect η_ijk^(1) to κ_ijk^(1). In particular we find

    η_ijk^(1)+η_ikj^(1) = 2κ_i(j+k+1)j+k+2j+1+2[κ_kji^(1)+κ_jki^(1)]

when combining the anti-symmetry property of κ_ijk^(1) with the gBRST identity. Similarly the combination of the Jacobi identity with gBRST leads to

    η_ijk^(1)+η_kij^(1)+η_jki^(1) = 2κ_i(j+k+1)j+k+1j+2κ_k(i+j+1)i+j+1i+2κ_j(i+k+1)i+k+1k.

The latter identity relates the class III coupling η_ijk^(1), which is O(g_s^2), 
to the class II coupling κ_ij of O(g_s), i.e. at one order lower in perturbation theory. 
As such, we can use it to determine the function space of η_ijk^(1). Taking into account all independent permutations of i, j and k, we find that this function space is 18-dimensional

    {   (-1)^i+ji+j+1i,N-2k+1i+j+1i,N-2ki+j+1i,(-1)^j+kj+k+1j,
       N-2i+1j+k+1j,N-2ij+k+1j,(-1)^i+ki+k+1k,N-2j+1i+k+1k,
       N-2ji+k+1k +  independent permutations of i, j and k}.

Furthermore, due to the close relationship between η_ijk^(1) and κ_ijk^(1), we assume that the functional form of the latter is constructed from the same functions.  
Hence in total we have 36 free parameters. Using the relations described above, cf. Eqs. (<ref>)-(<ref>), we are able to fix 34 of these. The final two free parameters are then determined using the one-loop results κ_110^(1)=0 and κ_121^(1)=13 C_A/336, which follow from the explicit operator renormalization for N=6 and N=8 respectively. Our final result for κ_ijk^(1) then becomes

    κ_ijk^(1) = η(N)/48{2(-1)^i + ji+j+1 i+ (-1)^i + ki+k+1 k + 3(-1)^j + k+1j+k+1j +i+k+1 i[ 2(-1)^i + k+1 + 5N-1j+1] + j+k+1k[3(-1)^j + k - 
       10N-2i + 4N-2i+1] + 
       i+j+1j[(-1)^i+j+1 + 5N-2k - 9N-2k+1] }.

To verify this expression, agreement with explicitly computed fixed-N values has been established up to N=20. The necessary direct computations at fixed values of N of Feynman diagrams for the OMEs 
with (physical or alien) spin-N twist-two operators ^(N) inserted 
in Green's functions with off-shell quarks, gluons or ghosts are performed with the setup used and described in <cit.> for the computation of moments of four-loop QCD splitting functions. 
In particular, the  Forcer package <cit.>, written in  Form <cit.>, is used for the parametric reductions of the two-point functions up to four loops for fixed even integer values of N. Substituting our result for κ_ijk^(1) into the gBRST relation, Eq. (<ref>), allows one to also reconstruct the full N-dependence of η_ijk^(1)
    η_ijk^(1) = -η(N)/24{5(-1)^i+j+1i+j+1i+(-1)^i+ki+k+1k+2(-1)^j+k+1j+k+1j+i+k+1i[(-1)^i+k+4N-2j+1]+j+k+1k[5(-1)^j+k+1-3N-2i+N-2i+1]+i+j+1j[4(-1)^i+j-15N-2k-5N-2k+1]}.

We have verified that Eqs. (<ref>) and (<ref>) are in agreement with the results of ref. <cit.>, as explained in Sec. <ref> below.


  §.§.§ κ_ijk^(2) and η_ijk^(2)

The next alien couplings we consider are κ_ijk^(2) and η_ijk^(2), which obey the following relations

    κ_ijk^(2)=κ_jik^(2)=κ_ikj^(2)=κ_kji^(2)=κ_jki^(2)=κ_kij^(2),    [symmetry of d_4]
       η_ijk^(2) = 3κ_ijk^(2),    [gBRST] 
       η_ijk^(2) = ∑_m=0^i∑_n=0^j(-1)^m+n+k(m+n+k)!/m!n!k!η_(i-m)(j-n)(m+n+k)^(2).    [anti-gBRST]

As the anti-gBRST equation has a similar form as the one for η_ijk^(1), cf. Eq. (<ref>), we assume the function space for η_ijk^(2) and κ_ijk^(2) to be the same as above, cf. Eq. (<ref>). Imposing Eqs. (<ref>)-(<ref>) then allows one to fix all but one of the unknowns. Hence we find expressions η_ijk^(2) and κ_ijk^(2) with only one (overall) free parameter

    κ_ijk^(2) = c{(-1)^i+ji+j+2i+1+(-1)^i+ki+k+2i+1+j+k+2j+1[(-1)^j+k+N-1i+1]},
       η_ijk^(2) = 3κ_ijk^(2).

Note that the c parameter can a priori be some N-dependent function. A computation of the OMEs at a few fixed values of N with the procedure outlined in Sec. <ref> for the renormalization of the respective operators fixes c=1/N/(N-1) such that

    κ_ijk^(2) = 1/N(N-1){(-1)^i+ji+j+2i+1+(-1)^i+ki+k+2i+1+j+k+2j+1[(-1)^j+k+N-1i+1]},
       η_ijk^(2) = 3κ_ijk^(2).

Noting that

    1/N(N-1) = -η(N)

we see that also for these two couplings η(N) factorizes, which is not expected a priori from the constraints.



  §.§.§ κ_ijk^(3) and η_ijk^(3)

The last set of couplings in the class III alien operators, κ_ijk^(3) and η_ijk^(3), obey the following relations

    κ_ijk^(3)=κ_ikj^(3),    [symmetry] 
       κ_ijk^(3)+κ_kij^(3)+κ_jki^(3)=0,    [generalized Jacobi identity] 
       η_ijk^(3) = 2(κ_ijk^(3)-κ_kji^(3)),   [gBRST]
       η_ijk^(3) = ∑_m=0^i∑_n=0^j(-1)^m+n+k(m+n+k)!/m!n!k!η_(j-n)(i-m)(m+n+k)^(3).    [anti-gBRST]

As before, we suggest the same function space as for κ_ijk^(1) and η_ijk^(1), cf. Eq. (<ref>). The above relations then only leave two parameters unfixed, such that we have

    κ_ijk^(3) = c_1 (-1)^i+ji+j+1i+c_2 (-1)^i+ki+k+1k+j+k+1j[(c_1+c_2) (-1)^j+k+1 
    +(c_1+c_2) (-1)^j+k+1 +c_1 N-2i+1]+i+k+1i[c_1 (-1)^i+k+(2 c_1+c_2) 
    N-2j+c_2 N-2j+1]+i+j+1j[c_2 (-1)^i+j-(2 c_1+c_2) N-2k-(c_1+c_2)  N-2k+1]

with c_1, c_2 to be determined. We emphasize that, as before, these could be N-dependent functions [In this case we expect c_1∼ c_2∼η(N).]. The corresponding expression for η_ijk^(3) depends on the same parameters through the gBRST relation, cf. Eq. (<ref>).
Since the couplings κ_ijk^(3) and η_ijk^(3) do not appear through operator mixing in the renormalization of physical OMEs up to four loops, we leave the two free parameters c_1, c_2 in Eq. (<ref>) undetermined, for the time being. 
We will address this issue again when extending the computation of low-N non-singlet anomalous dimensions at five loops <cit.> to the flavor-singlet sector.




 §.§ Class IV couplings


  §.§.§ κ_ijkl^(1) and η_ijkl^(1)

We have the following set of relations

    κ_ijkl^(1)+κ_ijlk^(1) = 0,   [anti-symmetry] 
       κ_ijkl^(1)+κ_iklj^(1)+κ_iljk^(1) = 0,    [Jacobi] 
       κ_ijkl^(1)+κ_jilk^(1)+κ_lkji^(1)+κ_klij^(1) = 0,   [double Jacobi] 
       η_ijkl^(1) = 2[κ_ij(l+k+1)^(1)+κ_(l+k+1)ji^(1)]l+k+1k+2[κ_ijkl^(1)+κ_ilkj^(1)+κ_likj^(1)+κ_lkij^(1)],    [gBRST] 
       η_ijkl^(1) = -∑_s_1=0^i∑_s_2=0^j∑_s_3=0^k(s_1+s_2+s_3+l)!/s_1!s_2!s_3!l!(-1)^s_1+s_2+s_3+lη_(k-s_3)(j-s_2)(i-s_1)(s_1+s_2+s_3+l)^(1)     [anti-gBRST]

with now i+j+k+l=N-5.
Combining the double Jacobi identity, Eq. (<ref>), with the gBRST one, Eq. (<ref>), allows one to write η_ijkl^(1) in terms of κ_ijk^(1) appearing already in the class III operators at one order in perturbation theory lower,

    η_ijkl^(1) + η_jilk^(1) + η_lkji^(1) + η_klij^(1) = 2[κ_ij(k+l+1)^(1)+κ_(k+l+1)ji^(1)]k+l+1k + 2[κ_ji(k+l+1)^(1)+κ_(k+l+1)ij^(1)]k+l+1l  + 2[κ_lk(i+j+1)^(1)+κ_(i+j+1)kl^(1)]i+j+1j  + 2[κ_kl(i+j+1)^(1)+κ_(i+j+1)lk^(1)]i+j+1i.

As such, we can use the expression we have computed for κ_ijk^(1), cf. Eq. (<ref>), to determine the function space of η_ijkl^(1). Taking into account all the independent permutations of the indices i, k, j and l this space is now 264-dimensional. Assuming that the functional form of κ_ijkl^(1) is similar to the one of η_ijkl^(1) then implies that in total we now have 528 parameters to fix. However, after implementing all of the above relations, only 8 remain in the end. The latter can again be fixed from the explicit renormalization of a few fixed-N matrix elements. More specifically we extracted them by performing a small momentum expansion around the limit p_3,p_4,p_5 → 0 of the OME 

    ⟨ O_ g^(N); c̅(p_1) c(p_2) g(p_3) g(p_4) g(p_5))⟩ .

This expansion is achieved on a diagram-by-diagram basis using the expansion-by-subgraph method <cit.> 
to second order at N=10 and to third order at N=12. By expanding sequentially in the external gluon momenta p_3,p_4 and p_5  the integrals are reduced to simple one-scale propagator integrals. We have implemented the expansion-by-subgraph in  Maple<cit.> and then subsequently evaluated the expressions in  Form. This methodology was also used to cross-check the expressions for κ_ij^(1) and κ_ijk^(r=1,2) up to N=20. At one loop the poles of the OME, Eq. (<ref>), are generated purely by the ghost alien operator O_c^(N),IV allowing for a clean extraction of η_ijkl^(1) renormalization constants, from which the  κ_ijkl^(1) values can be obtained. In particular, in order to determine the remaining constants in the all-N ansatz for κ_ijkl^(1), we use

    κ^(1)_0210=-1/128C_A,  
        κ^(1)_0050=109/1440C_A,  
        κ^(1)_0104=-935/6912C_A,  
        κ^(1)_1006=-2537/16896C_A.

We then find
    κ_ijkl^(1)=-η(N)/384{[6(-1)^j+ki+l+1i-3(-1)^j+ki+l+1l+7(-1)^j+k+lj+k+l+2l+7(-1)^j+k+lj+k+l+2j+k+1-27N-1i+1j+k+l+2j+k+1+2i+j+k+2j+k+1[2(-1)^i+j+k+9N-2l]-2i+j+k+2i[4(-1)^i+j+k+15N-2l+1]]j+k+1j-[5(-1)^j+ki+l+1i-4(-1)^j+ki+l+1l-14(-1)^j+k+lj+k+l+2l+7(-1)^j+k+lj+k+l+2j+k+1+54N-2ij+k+l+2j+k+1+i+j+k+2j+k+1[-3(-1)^i+j+k+4N-1l+1]+i+j+k+2i[3(-1)^i+j+k+13N-1l+1]]j+k+1k+(-1)^i+j+k+1i+k+1ki+j+k+2j-2(-1)^i+j+ki+j+1ji+j+k+2k-6(-1)^j+li+k+1ij+l+1j+3(-1)^j+li+k+1kj+l+1j+5(-1)^j+li+k+1ij+l+1l-4(-1)^j+li+k+1kj+l+1l+30N-2k+1j+l+1ji+j+l+2i-5(-1)^i+j+lj+l+1li+j+l+2i+13N-1k+1j+l+1li+j+l+2i-4(-1)^i+j+li+l+1ii+j+l+2j+12N-2ki+l+1ii+j+l+2j+24N-2k+1i+l+1ii+j+l+2j-3(-1)^i+j+li+l+1li+j+l+2j+49N-1k+1i+l+1li+j+l+2j-30N-1k+1i+j+1ii+j+l+2l+2(-1)^i+j+li+j+1ji+j+l+2l-60N-2k+1i+j+1ji+j+l+2l+8(-1)^i+j+li+l+1ii+j+l+2i+l+1-6N-2k+1i+l+1ii+j+l+2i+l+1-3(-1)^i+j+li+l+1li+j+l+2i+l+1+71N-1k+1i+l+1li+j+l+2i+l+1-11(-1)^k+li+j+1ik+l+1k+7(-1)^k+li+j+1jk+l+1k+11(-1)^k+li+j+1ik+l+1l-7(-1)^k+li+j+1jk+l+1l+60N-2j+1k+l+1ki+k+l+2i-10(-1)^i+k+lk+l+1li+k+l+2i+26N-1j+1k+l+1li+k+l+2i+(-1)^i+k+l+1i+l+1ii+k+l+2k-4N-1j+1i+l+1ii+k+l+2k-2(-1)^i+k+li+l+1li+k+l+2k-44N-2ji+l+1li+k+l+2k-26N-2j+1i+l+1li+k+l+2k-15N-1j+1i+k+1ii+k+l+2l+(-1)^i+k+li+k+1ki+k+l+2l-30N-2j+1i+k+1ki+k+l+2l+5(-1)^i+k+li+l+1ii+k+l+2i+l+1-10N-1j+1i+l+1ii+k+l+2i+l+1+(-1)^i+k+li+l+1li+k+l+2i+l+1-18N-2j+1i+l+1li+k+l+2i+l+1-14(-1)^j+k+lk+l+1lj+k+l+2j-7(-1)^j+k+lj+l+1lj+k+l+2k}

and

    η_ijkl^(1)=-η(N)/96{[[(j+k+l+2l+5j+k+l+2j+k+1)(-1)^l+1+3i+l+1i+3i+l+1l](-1)^j+k-17(-1)^i+j+ki+j+k+2i+i+j+k+2j+k+1[13(-1)^i+j+k+54N-2l]]j+k+1j+[-3(-1)^j+ki+l+1i-3(-1)^j+ki+l+1l+17(-1)^j+k+lj+k+l+2l+7(-1)^j+k+lj+k+l+2j+k+1+6N-2ij+k+l+2j+k+1+i+j+k+2i[(-1)^i+j+k+6N-1l+1]+i+j+k+2j+k+1[(-1)^i+j+k+6N-1l+1]]j+k+1k-12(-1)^i+j+ki+j+1ji+j+k+2k+(-1)^j+l+1i+k+1ij+l+1j+(-1)^j+l+1i+k+1kj+l+1j+(-1)^j+l+1i+k+1ij+l+1l+(-1)^j+l+1i+k+1kj+l+1l-3(-1)^i+j+li+l+1ii+j+l+2j-18N-2ki+l+1ii+j+l+2j+18N-2k+1i+l+1ii+j+l+2j-3(-1)^i+j+li+l+1li+j+l+2j+3N-1k+1i+l+1li+j+l+2j+18N-1k+1i+j+1ii+j+l+2l-6(-1)^i+j+li+j+1ji+j+l+2l-30N-2k+1i+j+1ji+j+l+2l+3(-1)^i+j+li+l+1ii+j+l+2i+l+1+12N-2k+1i+l+1ii+j+l+2i+l+1+3(-1)^i+j+li+l+1li+j+l+2i+l+1+3N-1k+1i+l+1li+j+l+2i+l+1+7(-1)^k+li+j+1ik+l+1k-5(-1)^k+li+j+1jk+l+1k+17(-1)^k+li+j+1ik+l+1l-13(-1)^k+li+j+1jk+l+1l-18N-2j+1k+l+1ki+k+l+2i-2(-1)^i+k+lk+l+1li+k+l+2i-10N-1j+1k+l+1li+k+l+2i+(-1)^i+k+l+1i+l+1ii+k+l+2k-5N-1j+1i+l+1ii+k+l+2k+(-1)^i+k+l+1i+l+1li+k+l+2k+10N-2ji+l+1li+k+l+2k-26N-2j+1i+l+1li+k+l+2k+(-1)^i+k+l+1i+l+1ii+k+l+2i+l+1-5N-1j+1i+l+1ii+k+l+2i+l+1-3(-1)^i+k+li+l+1li+k+l+2i+l+1-24N-2j+1i+l+1li+k+l+2i+l+1+4(-1)^j+k+lk+l+1lj+k+l+2j}.

We have checked the correctness of these expressions by comparing with fixed-N computations up to N=14.



  §.§.§ κ_ijkl^(2), η_ijkl^(2a) and η_ijkl^(2b)

For this final set of couplings we have the following relations

    κ_ijkl^(2)+κ_ijlk^(2) = 0,   [anti-symmetry] 
       κ_ijkl^(2) = κ_jikl^(2),   [symmetry of d_4] 
       η_ijkl^(2a) = 3κ_ij(k+l+1)^(2)k+l+1k+2κ_ijkl^(2),    [gBRST (a)] 
       η_ijkl^(2b) = 2κ_lijk^(2),    [gBRST (b)] 
       η_ijkl^(2a) = -∑_s_1=0^i∑_s_2=0^j∑_s_3=0^k(s_1+s_2+s_3+l)!/s_1!s_2!s_3!l!(-1)^s_1+s_2+s_3+l ×
       × η_(i-s_1)(j-s_2)(k-s_3)(s_1+s_2+s_3+l)^(2a),     [anti-gBRST (a)]  
       η_ijkl^(2b) = η_ikjl^(2a) - η_ijkl^(2a) + ∑_s_1=0^i∑_s_2=0^j∑_s_3=0^k(s_1+s_2+s_3+l)!/s_1!s_2!s_3!l!(-1)^s_1+s_2+s_3+l ×
       × η_(i-s_1)(j-s_2)(k-s_3)(s_1+s_2+s_3+l)^(2b).     [anti-gBRST (b)]

Note that Eqs. (<ref>) and (<ref>) can be combined to express η_ijkl^(2a) in terms of the class III coupling κ_ijk^(2) as

    η_ijkl^(2a) + η_ijlk^(2a) = 3κ_ij(k+l+1)^(2)k+l+2k+1.

Using the expression we derived for κ_ijk^(2), cf. Eq. (<ref>), this becomes

    η_ijkl^(2a) + η_ijlk^(2a) = 3c{(-1)^i+ji+j+2i+1-(-1)^i+k+li+k+l+3i+1+j+k+l+3j+1[-(-1)^j+k+l+N-1i+1]}k+l+2k+1

with c to be determined. Likewise one can use Eqs. (<ref>) and (<ref>) to write

    η_ijkl^(2a) - η_jikl^(2a) = 0.

To obtain this last identity we used the symmetry property of κ_ijk^(2), cf. Eq. (<ref>).
The complete solution of the gBRST constraints in Eqs. (<ref>)-(<ref>) proceeds in complete analogy to the previous cases.
However, similar to the class III couplings in Sec. <ref>, also 
κ_ijkl^(2), η_ijkl^(2a) and η_ijkl^(2b) 
do not enter in the operator renormalization of physical OMEs up to four loops, 
hence we will not consider them further here.




§ FEYNMAN RULES OF ALIEN OPERATORS

In this section we derive the Feynman rules of the alien operators. These were computed up to two loops in <cit.>, and an extension to the three-loop level was recently presented in <cit.>. 
The Feynman rules for the gauge-invariant (physical) quark and gluon operators, up to the four-loop level, can be found e.g. in <cit.> and references therein. 
The generalization to arbitrary orders in perturbation theory is given in <cit.>[Note that  <cit.> also presents the corresponding rules for the operators with total derivatives, relevant for non-zero momentum flow through the operator vertex.]. 
We assume all momenta to be incoming and the total momentum flowing through the operator vertex to be zero, implying

    ∑_ip_i = 0.


 §.§ Ghost operators

The momenta of the ghost fields are taken to be p_1 and p_2, while p_3,p_4,… denote the momenta of any additional gluons. 
As a check, we will compare our Feynman rules against the known ghost vertices with up to two additional gluons, which were computed in <cit.>. 
Because of different conventions for the operator definitions, the rules for the ghost vertices in the latter have to be divided by i^N. We can write the perturbative expansion of the ghost operator, cf. Eq. (<ref>), as

with

    ε^ab = 1+(-1)^N/2i^Nη(N)δ^ab(Δ· p_1)^N,

    ε^ab,c_1_μ = 1+(-1)^N/2i^N-1f^a c_1 b∑_i+j
    =N-3η_ij(Δ· p_1)(Δ· p_3)^i(Δ· p_2)^j+1,

    ε^ab,c_1 c_2_μν(p_1,p_2,p_3,p_4) = 1+(-1)^N/2i^NΔ_μΔ_ν{(f f)^a c_1 c_2 b∑_i+j+k
    =N-4η_ijk^(1)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_2)^k+1  +d_4^a c_1 c_2 b∑_i+j+k
    =N-4η_ijk^(2)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_2)^k+1 +d_4ff^a c_1 c_2 b∑_i+j+k
    =N-4η_ijk^(3)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_2)^k+1} + [(p_3,μ,c_1)↔ (p_4,ν,c_2)],

    ε^ab,c_1 c_2 c_3_μνρ(p_1,p_2,p_3,p_4,p_5) = -1+(-1)^N/2i^N-1{(f f f)^a c_1 c_2 c_3 b∑_i+j+k+l
    =N-5η_ijkl^(1)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j×(Δ· p_5)^k(Δ· p_2)^l+1+d_4f^a c_1 c_2 c_3 b∑_i+j+k+l
    =N-5η_ijkl^(2a)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k(Δ· p_2)^l+1+d_4f^a b c_1 c_2 c_3∑_i+j+k+l
    =N-5η_ijkl^(2b)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k(Δ· p_2)^l+1} +   permutations

where the `+ permutations' in the O(g_s^3) rule in Eq. (<ref>) denotes the fact that all permutations of the gluonic quantities (momenta, Lorentz and colour indices) have to be added. Note that p_2 in ε^ab in Eq. (<ref>) was eliminated using momentum conservation, p_2=-p_1. This then agrees with Eq. (5.20) in <cit.> after dividing the latter by i^N, as discussed above. 
Similarly, after performing the summation, the O(g_s) rule exactly matches Eq. (5.21) in <cit.>. 
At O(g_s^2), our expression for ε^ab,c_1 c_2_μν(p_1,p_2,p_3,p_4) should be compared against Eq. (5.22) in <cit.>. 
With η_ijk^(1) given by Eq. (<ref>) and η_ijk^(2) by Eq. (<ref>), we find exact agreement with that expression [The term in our expression proportional to (f f)^a c_1 c_2 b should be compared to the f^a_1 a_3 af^a_2 a_4 a part of Eq. (5.22)  in <cit.> while our d_4^a c_1 c_2 b rule should be compared to the one proportional to d_4^a_1 a_2 a_3 a_4/C_A.]. 
Finally the Feynman rule for the d_4ff part of the ghost operator is computed using

    η_ijk^(3) = 2{(c_1-c_2)(-1)^i+k[i+k+1i-i+k+1k]+(-1)^j+k+1[(c_1+2c_2)j+k+1j+(2c_1+c_2)j+k+1k]+i+j+1i[(2c_1+c_2)(-1)^i+j-c_1N-2k+1]+i+j+1j[c_1(-1)^i+j+2c_2(-1)^i+j-2c_1N-2k-c_1N-2k+1-c_2N-1k+1]+i+k+1i[2c_1N-2j+c_2N-1j+1]-i+k+1k[2c_1N-2j+c_2N-1j+1]+j+k+1j[2c_1N-2i+c_1N-2i+1+c_2N-1i+1]+c_1N-2i+1j+k+1k}

which follows from Eqs. (<ref>) and (<ref>). 
As discussed in Sec. <ref>, the free parameters c_1,c_2 can be determined by a computation of fixed-N OMEs.




 §.§ Alien gluon operators

Next we derive the Feynman rules for the gluonic EOM operator, whose perturbative expansion can be written as

with

    𝒢_μν^c_1 c_2(p_1,p_2) = 1+(-1)^N/2i^Nη(N)δ^c_1 c_2(Δ· p_1)^N-2[2p_1^2Δ_μΔ_ν-(Δ· p_1)(Δ_μp_1ν+Δ_νp_1μ)],

    𝒢_μνρ^c_1 c_2 c_3(p_1,p_2,p_3) = -1+(-1)^N/2i^N-1f^c_1 c_2 c_3{η(N)(Δ· p_1)^N-2Δ_μ[p_3νΔ_ρ-g_νρ(Δ· p_3)+Δ_ρ(p_2+p_3)_ν]+Δ_νΔ_ρ[p_1^2Δ_μ-p_1μ(Δ· p_1)]∑_i+j
    =N-3κ_ij(Δ· p_2)^i (Δ· p_3)^j }+ permutations,

    𝒢_μνρσ^c_1 c_2 c_3 c_4(p_1,p_2,p_3,p_4) = 1+(-1)^N/2i^N-2f^c_1 c_2 xf^x c_3 c_4{[Δ_νΔ_ρΔ_σ(p_1+2p_2)_μ-g_μνΔ_ρΔ_σ(Δ· p_2)]∑_i+j
    =N-3κ_ij(Δ· p_3)^i(Δ· p_4)^j-g_νρΔ_μΔ_σ(Δ· p_1)^N-2+[p_1^2Δ_μ-p_1μ(Δ· p_1)]Δ_νΔ_ρΔ_σ∑_i+j+k
    =N-4κ^(1)_ijk(Δ· p_2)^i(Δ· p_3)^j(Δ· p_4)^j}+1+(-1)^N/2[p_1^2Δ_μ-p_1μ(Δ· p_1)]Δ_νΔ_ρΔ_σ{d_4^c_1 c_2 c_3 c_4∑_i+j+k
    =N-4κ^(2)_ijk(Δ· p_2)^i(Δ· p_3)^j(Δ· p_4)^j+d_4ff^c_1 c_2 c_3 c_4∑_i+j+k
    =N-4κ^(3)_ijk(Δ· p_2)^i(Δ· p_3)^j(Δ· p_4)^j}+ permutations,

    𝒢_μνρστ^c_1 c_2 c_3 c_4 c_5(p_1,p_2,p_3,p_4,p_5) = 1+(-1)^N/2i^N-1f^c_1 c_2 xf^x c_3 yf^y c_4 c_5{-g_μρΔ_νΔ_σΔ_τ∑_i+j
    =N-3κ_ij(Δ· p_4)^i(Δ· p_5)^j+Δ_ρΔ_σΔ_τ[(p_1+2p_2)_μΔ_ν-(Δ· p_2)g_μν]∑_i+j+k
    =N-4κ_ijk^(1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k+[p_1^2Δ_μ-p_1μ(Δ· p_1)]Δ_νΔ_ρΔ_σΔ_τ∑_i+j+k+l
    =N-5κ_ijkl^(1)(Δ· p_2)^i(Δ· p_3)^j(Δ· p_4)^k(Δ· p_5)^l}+1+(-1)^N/2i^N-1d_4f^c_1 c_2 c_3 c_4 c_5{Δ_μΔ_νΔ_ρ[(p_4+2p_5)_σΔ_τ-(Δ· p_5)g_στ]∑_i+j+k
    =N-4κ_ijk^(2)(Δ· p_1)^i(Δ· p_2)^j(Δ· p_3)^k+[p_1^2Δ_μ-p_1μ(Δ· p_1)]Δ_νΔ_ρΔ_σΔ_τ∑_i+j+k+l
    =N-5κ_ijkl^(2)(Δ· p_2)^i(Δ· p_3)^j(Δ· p_4)^k(Δ· p_5)^l}+ permutations,

where again all permutations of gluon momenta, Lorentz and colour indices have to be added, if indicated by `+ permutations'. Note that p_2 in 𝒢_μν^c_1 c_2(p_1,p_2) in Eq. (<ref>) was again eliminated using momentum conservation. This then agrees with Eq. (5.23) in <cit.> and Eq. (243) in <cit.> after dividing the latter rules by i^N to match to our conventions. 

For the O(g_s) EOM vertex three contributions need to be taken into account,

  *  the non-Abelian part of the field strength in the class I operator with D→∂,
    
  *  the O(g_s) part of the covariant derivative acting on the Abelian part of the field strength in the class I operator and
    
  *  the class II operator, cf. Eq. (<ref>), with D→∂ and keeping only the Abelian part of the field strength.

Our result matches the corresponding rules in the literature, cf. Eq. (5.24) in <cit.> and Eq. (244) in <cit.> respectively (again after dividing by the overall i^N). 

Next the four-gluon vertex gets four contributions,

  *  the O(g_s) part of the covariant derivative acting on the non-Abelian part of the field strength in the class I operator,
    
  *  the non-Abelian part of the field strength in the class II operator with D→∂,
    
  *  the O(g_s) part of the covariant derivative acting on the Abelian part of the field strength in the class II operator and
    
  *  the class III operator, cf. Eq. (<ref>), with D→∂ and keeping only the Abelian part of the field strength.

The second and third contributions depend on the lower-order coupling κ_ij, while the fourth one is written in terms of the couplings κ_ijk^(1), κ_ijk^(2) and κ_ijk^(3) given by Eqs. (<ref>), (<ref>) and (<ref>) respectively. The (f f) and d_4 parts of our rule agree with Eq. (5.25) in <cit.>[Note however that our result proportional to d_4 needs to be multiplied by a symmetry factor of 1/4! to match Eq. (5.25) in <cit.>.], while the d_4ff part is new. 

Finally, as a new result [The corresponding result within the framework of ref. <cit.> was recently announced in a conference talk <cit.>.], we consider the five-gluon vertex 
  𝒢_μνρστ^c_1 c_2 c_3 c_4 c_5(p_1,p_2,p_3,p_4,p_5) 
in Eq. (<ref>). 
Again we need to take into account higher-order contributions of the lower-point vertices. In particular, the (f f f) part of the five-gluon rule gets four contributions,

  *  the O(g_s) part of the covariant derivative acting on the non-Abelian part of the field strength in the class II operator,
    
  *  the non-Abelian part of the field strength in the class III operator with D→∂,
    
  *  the O(g_s) part of the covariant derivative acting on the Abelian part of the field strength in the class III operator,
    
  *  the class IV operator, cf. Eq. (<ref>), with D→∂ and keeping only the Abelian part of the field strength.

On the other hand the d_4f part only gets three contributions,

  *  the non-Abelian part of the field strength in the class III operator with D→∂,
    
  *  the O(g_s) part of the covariant derivative acting on the Abelian part of the field strength in the class III operator and
    
  *  the class IV operator, cf. Eq. (<ref>), with D→∂ and keeping only the Abelian part of the field strength.


 §.§ Alien quark operators

Finally in this section we provide the Feynman rules for the alien quark operators presented in Eqs. (<ref>)-(<ref>). As mentioned above, these operators are written in terms of the same couplings as those in the gluon EOM operators. Assuming the momenta of the external quark fields to be p_1 and p_2 we have the following perturbative expansion

with

    𝒬(p_1,p_2) = 0,

    𝒬_  μ^c_1(p_1,p_2,p_3) = -1+(-1)^N/2i^Nη(N)T^c_1Δ_μΔ(Δ· p_3)^N-2,

    𝒬_  μν^c_1 c_2(p_1,p_2,p_3,p_4) = [1+(-1)^N]i^N-1T^af^a c_1 c_2Δ_μΔ_νΔ∑_i+j
    =N-3κ_ij^(1)(Δ· p_3)^i(Δ· p_4)^j,

    𝒬_  μνρ^c_1 c_2 c_3(p_1,p_2,p_3,p_4,p_5) =  [1+(-1)^N]i^NT^aΔ_μΔ_νΔ_ρΔ{f^a c_1 xf^c_2 c_3 x∑_i+j+k
    =N-4κ_ijk^(1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k
        + d_4^a c_1 c_2 c_3∑_i+j+k
    =N-4κ_ijk^(2)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k + d_4ff^a c_1 c_2 c_3∑_i+j+k
    =N-4κ_ijk^(3)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k}+ [(p_3,μ,c_1)↔(p_4,ν,c_2)] + [(p_3,μ,c_1)→(p_5,ρ,c_3)→(p_4,ν,c_2)→(p_3,μ,c_1)]

The vertices with up to two additional gluons can be compared against the results presented in Eqs. (5.17)-(5.19) of <cit.>. Dividing the latter by i^N to match to our conventions, we find exact agreement. 
Note that Eqs. (<ref>)-(<ref>) contain an additional factor of two coming from the [(p_4,ν,c_2)↔(p_5,ρ,c_3)] permutation. This directly follows from the (anti-)symmetry properties of the κ-couplings, cf. Eqs. (<ref>), (<ref>) and (<ref>). 
Finally, because the κ-couplings enter the quark operator at one order in the strong coupling lower than in the gluon EOM one, we can push the perturbative expansion of the quark operator to one order higher. 
Consequently we also present the quark operator vertex at O(g_s^4) with four additional gluons

We find

    𝒬^c_1 c_2 c_3 c_4_μνρσ(p_1,p_2,p_3,p_4,p_5,p_6) = -1+(-1)^N/2i^N-1T^aΔ_μΔ_νΔ_ρΔ_σΔ{(f f f)^a c_1 c_2 c_3 c_4∑_i+j+k+l
    =N-5κ_ijkl^(1)×(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k(Δ· p_2)^l+1+d_4f^a c_1 c_2 c_3 c_4∑_i+j+k+l
    =N-5κ_ijkl^(2)(Δ· p_1)(Δ· p_3)^i(Δ· p_4)^j(Δ· p_5)^k(Δ· p_2)^l+1} +   permutations

with all permutations of the gluonic quantities (momenta, Lorentz and colour indices) to be added.



§ CONCLUSIONS


The kernels for parton evolution equations in QCD, i.e. splitting functions or the corresponding anomalous dimensions as their Mellin transforms, can be conveniently determined from the ultraviolet singularities of off-shell Green's functions with insertions of gauge-invariant twist-two spin-N operators.
The renormalization of these OMEs, though, requires the computation of unphysical counterterms for the associated set of alien operators, which effectively describe vertices of two gluons, ghosts or quarks with any number n≥ 0 of additional gluons.
The couplings of these alien operators (EOM and ghost operators) are restricted by the fundamental symmetries, particularly the gBRST relations, which reflect the gauge theory characteristics of QCD.

The set of constraints for these couplings admits explicit solutions, valid for any spin N, which can be obtained using algorithms for symbolic summation to solve the recurrence relations.
A small number of boundary conditions in these solutions can be derived from the computation of the relevant OMEs at specific fixed values of N.
In addition, we have observed that the constraints contain a hierarchy, such that couplings of alien operators with n+1 gluons can be derived from those ones with only n gluons. Thus, the basic ingredients in this bootstrap turn out to be the EOM and ghost operators with the smallest number of additional gluons at a given loop order.

We have provided results for all one-loop alien operator couplings needed in the renormalization of OMEs with physical (gauge-invariant) operators up to four loops, which represents the current frontier in splitting function computations.
This includes in particular the gluon EOM operator with five gluons attached, which is a new result. 
The all-N solutions for the couplings that we have obtained can all be related to the fundamental one-loop counterterm η(N) for the EOM and ghost operators of class I involving only two gluons or ghosts. 
We have also derived the corresponding Feynman rules and, whenever possible, compared them to those in the literature, finding full agreement.

A  Mathematica file with our results for the all-N couplings necessary for the renormalization up to four loops is made available at the preprint server <https://arxiv.org>. We note that the expressions collected in this file have the fundamental one-loop counterterm η(N) divided out.

The symmetries and the structure of the alien operators, that we have exploited in this study, are  independent of the order of perturbation theory.
Thus, we expect also analytic all-N solutions beyond one loop for the couplings of the alien operators of class II and higher. 
We leave this task to future studies.






 §.§ Acknowledgments


The Feynman diagrams in this work are drawn using FeynGame <cit.>. 

This work has been supported by the EU's Marie Sklodowska-Curie 
grant 101104792,  QCDchallenge;
the DFG through the Research Unit FOR 2926,
 Next Generation pQCD for Hadron Structure: Preparing for the EIC,
project number 40824754, DFG grant MO 1801/4-2,
the ERC Advanced Grant 101095857  Conformal-EIC; 
and by grant K143451 of the National Research, Development and Innovation Fund in Hungary.

JHEP