  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 ($\overline 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
($\overline{MS}$) 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 \rightarrow g g g$ and $A \rightarrow
q \overline q g$, 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 ${\cal 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.
