@Fuks:2015hna; @deFlorian:2016sit; @deFlorian:2019app][^1], the trilinear Higgs coupling can be accessed directly in Higgs-pair production. Higgs-boson pairs are dominantly produced in the loop-induced gluon-fusion mechanism $gg\to HH$ that is mediated by top-quark loops supplemented by a per-cent-level contribution of bottom-quark loops, see Fig. \[fg:hhdia\]. There are destructively interfering box and triangle diagrams at LO with the latter involving the trilinear Higgs coupling [@Glover:1987nx; @Plehn:1996wb]. The box diagrams provide the dominant contributions to the cross section. A rough estimate of the dependence of the cross section on the size of the trilinear coupling is given by the approximate relation $\Delta\sigma/\sigma \sim -\Delta\lambda_{H^3}/\lambda_{H^3}$ in the vicinity of the SM value of $\lambda_{H^3}$. Therefore, in order to determine the trilinear coupling, the theoretical uncertainties of the corresponding cross section need to be small. Thus, the inclusion of higher-order corrections is mandatory. The QCD corrections are fully known up to next-to-leading order (NLO) [@Borowka:2016ehy; @Borowka:2016ypz; @Baglio:2018lrj] and at next-to-next-to-leading order (NNLO) in the limit of heavy top quarks [@deFlorian:2013uza; @deFlorian:2013jea; @Grigo:2014jma]. While the NLO corrections are large, the NNLO contributions are of more moderate size. Very recently, the next-to-next-to-next-to-leading order (N$^3$LO) QCD corrections have been computed in the limit of heavy top quarks resulting in 