In QCD the anomalous dimensions of gauge invariant operators of twist 2 play a key role, because
they control the scale dependence of the parton distribution functions. Notably, the flavour singlet
operators, such as those associated to the gluon distribution, mix under renormalisation with a
set of unphysical operators, also known as aliens. Missing this effect leads to wrong results
already for the two-loop anomalous dimensions. The correct renormalisation of gluonic operators
is an important step towards the computation of the scale evolution of flavour singlet parton
distributions, which is now required to 4 loops. Leveraging both the background field method and
an enhanced BRST symmetry, we construct the required ghost and alien operator basis up to 4
loops for arbitrary mass dimensions. Furthermore, we extract anomalous dimensions at 4 loops,
for the physical operators of mass-dimension 4 and 6, and at 3 loops for mass-dimension 8.