a straightforward manner. The analytical expressions can be found in Refs. [@Glover:1987nx; @Plehn:1996wb]. Working out the polarization and color sums of the matrix element of Eq. (\[eq:lomat\]), the LO partonic cross section $\hat\sigma_{LO}$ is given by $$\hat\sigma_{LO} = \frac{G_F^2\alpha_s^2(\mu_R)}{512 (2\pi)^3}
\int_{\hat t_-}^{\hat t_+} d\hat t \Big[ | F_1 |^2 + |F_2|^2
\Big]$$ with the integration boundaries $$\hat t_\pm = -\frac{1}{2} \left[ Q^2 - 2M_H^2 \mp Q^2
\sqrt{1-4\frac{M_H^2}{Q^2}} \right] \, ,
\label{eq:tbound}$$ where the symmetry factor 1/2 for the identical Higgs bosons in the final state is taken into account. The LO hadronic cross section $\sigma_{LO}$ can then be derived by a convolution with the parton densities $$\sigma_{LO} = \int_{\tau_0}^1 d\tau \frac{d{\cal L}^{gg}}{d\tau}
\hat\sigma_{LO}(Q^2 = \tau s)$$ with the gluon luminosity, given in terms of the gluon densities $g(x,\mu_F)$, $$\frac{d{\cal L}^{gg}}{d\tau} = \int_\tau^1 \frac{dx}{x} g(x,\mu_F)
g\left(\frac{\tau}{x},\mu_F\right)
\label{eq:lgg}$$ at the factorization scale $\mu_F$ and the integration boundary $\tau_0=4M_H^2/s$, where $s$ denotes the hadronic center-of-mass (c.m.) energy squared. The differential cross section with respect to the invariant squared Higgs-pair mass $Q^2$ can be obtained as $$\frac{d\sigma_{LO}}{dQ^2} = \left. \frac{d{\cal L}^{gg}}{d\tau}~
\frac{\hat\sigma_{LO}(Q^2)}{s} \right|_{\tau = \frac{Q^2}{s}} \, .
\label{eq:lodiff}$$ As can be expected from single Higgs-boson production via gluon fusion (see [@Graudenz:1992pv; @Spira:1995rr; @Harlander:2005rq; @Anastasiou:2009kn; 