Understanding the scale dependence of parton distribution functions is vital for precision physics
at hadron colliders. The well-known DGLAP evolution equation relates this scale dependence to
the QCD splitting functions, which can be calculated perturbatively in terms of the anomalous
dimensions of leading-twist gauge-invariant operators. The computation of the latter in general
requires one to take into account contributions of gauge-variant (or alien) operators. In this
talk, we discuss the systematic study of these alien operators at arbitrary spin. Specifically, using
generalized BRST symmetry relations, we derive the one-loop couplings and Feynman rules of the
aliens necessary to perform the operator renormalization up to four loops in QCD. This provides
an important step towards the determination of the four-loop splitting functions which will be of
significant phenomenological importance at future colliders.