We present virtual contributions up to two loop level in perturbative Quantum Chromodynamic (QCD) to the decay of pseudo-scalar Higgs boson (A) to three gluons (g) and also to quark (q), anti-quark (q) and a gluon. With appropriate crossing, they are well suited for predicting the differential distribution of A in association with a jet in hadron colliders up to next-to-next-to-leading order (NNLO) in strong coupling constant and also for the subsequent decay of A to hadrons. We use effective field theory approach to integrate out the top quarks in the heavy top limit. The resulting theory involves two pseudo-scalar composite operators describing the interaction of A with gluons as well as with quark and anti-quark. We perform our computation in dimensional regularisation and use minimal subtraction (M S) scheme to renormalise strong coupling constant as well as the composite operators. The ultraviolet (UV) finite amplitudes contain infrared (IR) divergences that are found to be in agreement with the predictions by Catani. For both the amplitudes namely A → ggg and A → qqg, the leading transcendental terms at one and two loops are found to be identical to those in a three point form factor (FF) of the half-BPS operator in N = 4 Supersymmetric Yang Mills (SYM) theory when the QCD color factors are adjusted in a specific way. We present our results in terms of harmonic polylogs well suited for further numerical study.