The Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting functions, which control the scale evolution of the parton distribution functions, are key ingredients for collider physics. Challenging precision targets posed by future experiments, for instance in the High Luminosity phase at the LHC, require us to compute the N3LO (or four-loop) splitting functions. To this end, it's promising to determine the Mellin moments of the DGLAP splitting functions, which correspond to anomalous dimensions of composite operators. The main bottleneck of this approach is the renormalisation of gluonic operators, which mix with a set of unphysical "alien" operators. In this talk, based on 2203.11181 (to appear in JHEP), I will discuss a systematic construction of all the required aliens and demonstrate how the procedure works, by computing physical anomalous dimensions for low Mellin moments at 3 and 4 loops.