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 γ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.