The coupling of a pseudo-scalar Higgs boson to gluons is mediated through a heavy quark loop. In the limit of large quark mass, it is described by an effective Lagrangian that only admits light degrees of freedom. In this effective theory, we compute the three-loop massless QCD corrections to the form factor that describes the coupling of a pseudo-scalar Higgs boson to gluons. Due to the axial anomaly, the pseudo-scalar operator for the gluonic field strength mixes with the divergence of the axial vector current. Working in dimensional regularization and using the 't Hooft-Veltman prescription for the axial vector current, we compute the three-loop pseudo-scalar form factors for massless quarks and gluons. Using the universal infrared factorization properties, we independently derive the three-loop operator mixing and finite operator renormalisation from the renormalisation group equation for the form factors, thereby confirming recent results in the operator product expansion. The finite part of the three-loop form factor is an important ingredient to the precise prediction of the pseudo-scalar Higgs boson production cross section at hadron colliders. We discuss potential applications and derive the hard matching coefficient in soft-collinear effective theory.