  Through this article, we present the two-loop massless QCD corrections to the
production of di-Higgs and di-pseudo-Higgs boson through quark annihilation in
the large top quark mass limit. Within dimensional regularisation, we employ
the non-anticommuting $\gamma_5$ and treat it under the Larin prescription. We
discover the absence of any additional renormalisation, so-called contact
renormalisation, that could arise from the short distance behaviour of two
local operators. This finding is in corroboration with the operator product
expansion. By examining the results, we discover the lack of similarity in the
highest transcendentality weight terms between these finite remainders and that
of a pair of half-BPS primary operators in maximally supersymmetric Yang-Mills
theory. We need these newly computed finite remainders to calculate observables
involving di-Higgs or di-pseudo-Higgs at the next-to-next-to-leading order. We
implement the results to a numerical code for further phenomenological studies.
