The virtual corrections for $gg \to HH$ at next-to-leading order (NLO) in QCD have been efficiently approximated using a Taylor expansion in the limit of a forward kinematics. The same method has been recently applied to the calculation of a subset of the next-to-NLO corrections, which are desirable given the significant impact, at NLO, of the uncertainty due to the choice of the top mass renormalization scheme. In these proceedings, we report on the progress in the calculation of another contribution at NNLO, given by diagrams in which the two Higgs bosons couple to different top quark loops. For this contribution a naive Taylor expansion cannot be used, and we discuss an approach based on asymptotic expansions in different kinematic limits.