integral [@Grigo:2013rya] and the integrand [@Grigo:2015dia]. This indicated sizable mass effects in the virtual two-loop corrections alone. In addition, the large top-mass expansion has been extended to the virtual NNLO QCD corrections resulting in 5% mass effects estimated on top of the NLO result [@Grigo:2015dia]. The large-top-mass expansion of the NLO QCD corrections has been used to perform a conformal mapping of the expansion parameter and to apply Padé approximants. In this way, an approximation of the full calculation has been achieved for $Q$ values up to about 700 GeV [@Grober:2017uho]. Another approximation builds on an expansion in terms of a variable that dominantly corresponds to the transverse momentum of the Higgs bosons. The results of this approach show good agreement with the full calculation for $Q$ values up to about 900 GeV [@Bonciani:2018omm]. Analytical results are also available in the large-$Q$ limit [@Davies:2018qvx]. The latter have recently been combined with the numerical results of Refs. [@Borowka:2016ehy; @Borowka:2016ypz] for the full QCD corrections [@Davies:2019dfy]. In the following, we will discuss the details of our NLO calculation.

Virtual corrections \[sc:virtuals\]
-----------------------------------

Typical diagrams of the two-loop virtual corrections are shown in Fig. \[fg:gghhvirt\]. They can be arranged in three different classes: 