the ultraviolet singularities for the numerical integration has first been applied to the NLO two-loop QCD corrections to $H\to\gamma\gamma, Z\gamma$ in Refs. [@Djouadi:1990aj; @Spira:1991tj] and later to the squark-loop contributions to $h,H \leftrightarrow
gg,\gamma\gamma$ within the minimal supersymmetric extension of the SM [@Muhlleitner:2006wx]. The method of the infrared subtraction as applied to topology 6 originates from numerical cross checks of the full NLO QCD corrections to single Higgs production in Refs. [@Spira:1995es; @Spira:1995rr; @Graudenz:1992pv; @Muhlleitner:2006wx]. The stabilization of virtual thresholds by integration by parts of the integrand has first been applied to the SUSY–QCD corrections to single Higgs production in Refs. [@Muhlleitner:2010nm; @Muhlleitner:2010zz]. The basic idea behind the integration by parts is to reduce the power of the threshold-singular denominator and in this way to stabilize the numerical integration. The treatment of the thresholds in our approach is performed by replacing the squared top mass $m_t^2$ by a complex counter part $$m_t^2 \to m_t^2 (1-i\bar\epsilon)
\label{eq:imaginary}$$ with a positive regulator $\bar\epsilon > 0$ to ensure proper micro-causality. This defines the analytical continuation of our two-loop box integrals. In the following, the parameter $\bar\epsilon$ will be kept finite in our numerical analysis, while the narrow-width limit $\bar\epsilon\to 0$ is achieved by 