Real-virtual corrections to Higgs boson pair production at NNLO: three closed top quark loops

TL;DR

This paper calculates specific real-virtual QCD corrections to Higgs pair production at NNLO, focusing on the interference effects involving top quark loops, providing analytic results for phase-space integrals.

Contribution

It presents the first analytic computation of real-virtual corrections to Higgs pair production involving three closed top quark loops at NNLO in QCD.

Findings

Analytic expressions for phase-space master integrals.

Expansion results around the production threshold.

Exact results using Goncharov polylogarithms.

Abstract

We compute the real-radiation corrections to Higgs boson pair production at next-to-next-to-leading order in QCD, in an expansion for large top quark mass. We concentrate on the radiative corrections to the interference contribution from the next-to-leading order one-particle reducible and the leading order amplitudes. This is a well defined and gauge invariant subset of the full real-virtual corrections to the inclusive cross section. We obtain analytic results for all phase-space master integrals both as an expansion around the threshold and in an exact manner in terms of Goncharov polylogarithms.

Equations

$$g(q_1)g(q_2)$$

$$g(q_1)q(q_2)$$

$$g(q_1)\bar{q}(q_2)$$

$$q(q_1)\bar{q}(q_2)$$

$${\sigma }_{ij\to HH+X}(s,\rho )={\delta }_{ig}{\delta }_{jg}{\sigma }_{gg}^{(0)}(s,\rho )+\frac{{\alpha }_s}{\pi }{\sigma }_{ij}^{(1)}(s,\rho )+{\left(\frac{{\alpha }_s}{\pi }\right)}^2{\sigma }_{ij}^{(2)}(s,\rho )+\dots ,$$

$$\rho$$

$$\delta$$

$${\sigma }_{gg}^{(0)}$$

$$\sum _{\lambda }{\epsilon }_{\mu }^{(\lambda )\star }(q_1){\epsilon }_{\nu }^{(\lambda )}(q_2)$$

$${\sigma }_{gg}^{(1)}$$

$${\sigma }_{qg}^{(1)}$$

$${\sigma }_{q\bar{q}}^{(1)}$$

$$\displaystyle\sigma_{ij\to HH+X}(s,\rho)\Big{|}_{n_{h}^{3}}$$

$${\sigma }_{ij}^{(2),n_h^3}$$

$$P_{ij}$$

$$P_{gg}^{(0)}(z)$$

$$P_{gq}^{(0)}(z)$$

$${\beta }_0=\frac{11}{12}{C}_A-\frac{1}{3}{T}_f{n}_l,$$

$${\sigma }_{gg,\mathrm{coll}}^{(2),n_h^3}$$

$${\int }_{1-\delta }^1\mathrm{d}z{\left[\frac{1}{1-z}\right]}_{+}{\sigma }_{gg}^{(1),n_h^3}(x/z)={\int }_{1-\delta }^1\mathrm{d}z\frac{{\sigma }_{gg}^{(1),n_h^3}(x/z)-{\sigma }_{gg}^{(1),n_h^3}(x)}{1-z}+{\sigma }_{gg}^{(1),n_h^3}(x)\ln (\delta ),$$

$${\sigma }_{gg}^{(1),n_h^3}=\sum _{n=0}^{\infty }{\delta }^{\frac{1}{2}+n}\left[c_n^{(0)}+ϵ(c_n^{(1)}-c_n^{(0)}\ln \delta )\right]+\mathcal{O}({ϵ}^2),$$

$${\sigma }_{gg,\mathrm{coll}}^{(2),n_h^3}$$

$$+c_j^{(0)}\left(\frac{{\Delta }_0}{(\frac{3}{2}+n{)}^2}-\frac{{\Delta }_1}{(\frac{5}{2}+n{)}^2}+\frac{{\Delta }_2}{(\frac{7}{2}+n{)}^2}-\frac{{\Delta }_3}{(\frac{9}{2}+n{)}^2}\right)\}$$

$$-2\frac{{\alpha }_s}{\pi }{C}_A(1-\delta )\sum _{n=0}^{\infty }\sum _{j=0}^n{\delta }^{\frac{1}{2}+n}\left\{{\tilde{c}}_j\left[\psi \left(\frac{3}{2}+n\right)-\psi (n-j+1)\right]-c_j^{(0)}{\psi }^{\prime }\left(\frac{3}{2}+n\right)\right\}$$

$$+2\frac{{\alpha }_s}{\pi }\left(C_A\ln \delta +{\beta }_0\right)\sum _{n=0}^{\infty }{\delta }^{\frac{1}{2}+n}{\tilde{c}}_n+\mathcal{O}(ϵ)$$

$${\Delta }_0$$

$${\Delta }_1$$

$${\Delta }_2$$

$${\Delta }_3$$

$${\tilde{c}}_j$$

$${\sigma }_{qg,\mathrm{coll}}^{(2),n_h^3}$$

$$=\frac{{\alpha }_s}{\pi }\frac{C_F}{4}\sum _{n=0}^{\infty }\sum _{j=0}^n{\delta }^{\frac{3}{2}+n}\frac{(n-j+2)!}{(n-j)!}\{c_j^{(0)}\left(\frac{1+{\delta }^2}{(\frac{3}{2}+n{)}^2}-\frac{2\delta (1+\delta )}{(\frac{5}{2}+n{)}^2}+\frac{2{\delta }^2}{(\frac{7}{2}+n{)}^2}\right)$$

$$+{\tilde{c}}_j\left(\frac{1+{\delta }^2}{\frac{3}{2}+n}-\frac{2\delta (1+\delta )}{\frac{5}{2}+n}+\frac{2{\delta }^2}{\frac{7}{2}+n}\right)\}+\mathcal{O}(ϵ).$$

$${\sigma }_{gg,\mathrm{coll}}^{(2),n_h^3}$$

$${\sigma }_{qg,\mathrm{coll}}^{(2),n_h^3}$$

$$a_s=\frac{{\alpha }_s^{(5)}({\mu }_r)}{\pi },L_{M_t^2}=\log \frac{{\mu }_r^2}{M_t^2},L_{m_H^2}=\log \frac{{\mu }_r^2}{m_H^2},L_{{\mu }_f^2}=\log \frac{{\mu }_r^2}{{\mu }_f^2}.$$

$${\sigma }_{gg,\mathrm{real}}^{(2),n_h^3}$$

$${\sigma }_{qg,\mathrm{real}}^{(2),n_h^3}$$

$${\sigma }_{q\bar{q},\mathrm{real}}^{(2),n_h^3}$$

Full text

1 Introduction

Within the Standard Model (SM) of particle physics, knowledge of the Higgs boson mass ($m_{H}$) fixes the parameters of the scalar potential, since the trilinear and quartic couplings can be expressed through $m_{H}$ and the vacuum expectation value of the scalar doublet. However, in many extensions of the SM, the trilinear and quartic couplings deviate significantly from the SM values. One way of probing the Higgs self coupling is through the production of Higgs boson pairs at the LHC. Besides the efforts undertaken on the experimental side, also precise predictions from the theoretical side are necessary to scrutinize the results of upcoming measurements. The most important production mode of a pair of Higgs bosons at hadron colliders is by top quark–mediated gluon fusion. In this channel the QCD corrections are large and higher-order computations are important.

The leading order (LO) cross section has been known with full top quark–mass dependence for more than 30 years [1, 2]. Exact next-to-leading order calculations are numerically quite challenging and have only become available fairly recently [3, 4, 5]. Analytic NLO calculations have so far only been performed for various approximations. Among them is the effective theory calculation in which the heavy top quark is integrated out [6]. This result has been extended in Refs. [7, 8] where inverse top quark mass corrections have been computed. An expansion for small top quark masses has been obtained in Refs. [9, 10, 11] and the limit of small transverse momentum is covered by the results of Ref. [12]. In Ref. [13] an approximation method for the reconstruction of the form factors has been suggested. Results from various kinematic regions are combined using conformal mapping and Padé approximation. Finite top quark mass effects to the real radiation contribution have been considered in Ref. [14].

At next-to-next-to-leading order (NNLO), the effective-theory calculation of the cross section has been performed in Refs. [15, 16, 17] and an expansion for large top quark masses has been performed in Ref. [18] in the soft-virtual approximation. Beyond the infinite top mass limit real radiations are missing so far; it is the aim of this paper to partly close this gap. Note that in the effective-theory calculation a large part of the corrections to Higgs boson pair production can be taken over from single Higgs production [19, 20, 21]; this is no longer the case once one goes beyond this approximation.

Let us mention that recently two building blocks of the next-to-next-to-next-to-leading order (N3LO) effective-theory result have become available: two-loop virtual corrections have been obtained in Ref. [22] and the four-loop matching coefficient for the effective coupling of two Higgs bosons and gluons has been computed in [23, 24].

In this paper we compute the imaginary part of the forward scattering amplitudes $ij\to ij$, where $i$ and $j$ stand for gluons and (anti-)quarks. With the help of the optical theorem one obtains the (partonic) cross section ${\rm d}\sigma/{\rm d}s$ where $\sqrt{s}$ is the center-of-mass energy of the incoming partons. Note that at the lowest order one already has to compute three-loop Feynman diagrams and at NkLO one has to consider $(k+3)$-loop diagrams. Sample Feynman diagrams can be found in Fig. 1. The virtual corrections are obtained from the contributions where exactly two Higgs bosons are cut. Note that besides the virtual corrections to the NLO 1PR diagram (Fig. 1 (e)) also diagrams such as Fig. 1 (j) have three closed top quark loops. At NLO the final state of the real radiation corrections contains two Higgs bosons and an additional parton. At NNLO one has either one or two additional partons in the final state. We refer to the former as “real-virtual” (Fig. 1 (f), (g), (h) and (k)) and the latter as “double-real” (Fig. 1 (l)).

The real-virtual corrections can be sub-divided according to the number of closed top quark loops which involve a coupling to one or two Higgs bosons. At NNLO this is either two or three, as can be seen from the Feynman diagrams in Fig. 1. We will refer to them as $n_{h}^{2}$ and $n_{h}^{3}$ contributions in the following. In this paper we consider only the $n_{h}^{3}$ contribution, with three closed top quark loops. In an asymptotic expansion in large $M_{t}$ all top quark lines are part of the so-called hard subgraphs, which means that the remaining Feynman diagrams which involve the Higgs bosons are either one- or two-loop diagrams.

At NLO, $n_{h}^{3}$ terms are only present in the virtual corrections, see Fig. 1(c). They serve as an effective LO contribution for the $n_{h}^{3}$ NNLO corrections we are interested in. In this sense, one can consider the subset of real-virtual corrections with three top quark loops as effective NLO real corrections. Thus, they share many features with the NLO real corrections and many steps of the calculation can be performed in analogy to Ref. [7].

We consider all partonic channels contributing to the cross section and compute the necessary phase-space integrals both as an expansion around the pair-production threshold [7] and also as expressions exact in $m_{H}$ and $s$. In the latter case the analytic results are expressed in terms of Goncharov polylogarithms [25]. As we will show, for all phenomenological applications it is sufficient to consider the expanded results which have a simpler mathematical structure.

While for single Higgs production top quark mass suppressed terms converge well after factoring out the exact LO cross section (see Refs. [26, 27] for the NNLO analysis), this is not the case for Higgs boson pair production [7]. However, the large top quark mass expansion can still provide valuable input for approximation methods in the full range of center of mass energies [13, 28].

The remainder of this paper is organized as follows. In the next Section we introduce our notation and comment on the techniques which have been used for the calculation. We discuss our results in Section 3 and present our conclusions in Section 4. In the Appendix we provide additional material such as details on the threshold expansion of the phase-space master integrals (Appendix A) and the calculation of the master integrals without expansion (Appendix B).

2 Preliminaries and technicalities

In this section we fix our notation and provide technical details on the calculation of the real-virtual $n_{h}^{3}$ contribution to the cross section for double Higgs production.

2.1 Cross section

We compute the cross section by applying the optical theorem to the forward scattering amplitudes

[TABLE]

where $g$ stands for gluons and $q$ and $\bar{q}$ represent generic (light) quarks and anti-quarks. $q_{1}$ and $q_{2}$ are the momenta of the in- and outgoing partons with $q_{1}^{2}=q_{2}^{2}=0$. Note the at NNLO, the partonic channel which involves different quark flavours in the initial state does not yet contribute. Of course, we only take into account diagrams which involve cuts through two internal Higgs boson propagators. Sample Feynman diagrams for the amplitudes in Eq. (1) can be found in Fig. 1. The forward scattering amplitudes depend on the top quark mass $M_{t}$, the center-of-mass energy $\sqrt{s}$ with $s=(q_{1}+q_{2})^{2}=2q_{1}\cdot q_{2}$, and the Higgs boson mass $m_{H}$.

In order to fix the notation and the pre-factors we now discuss the LO and NLO cross sections in detail. The corresponding formulae for the $n_{h}^{3}$ NNLO contributions are obtained by straightforward replacements.

We write the perturbative expansion of the partonic cross section for Higgs boson pair production as

[TABLE]

where $\alpha_{s}\equiv\alpha_{s}^{(5)}(\mu_{r})$ is the strong coupling constant in the five-flavour theory and $ij\in\{gg,qg,\bar{q}g,q\bar{q}\}$ denote the partonic sub-channels. We introduce the variable

[TABLE]

to parametrize the dependence of the cross section on the Higgs boson and top quark masses. For later convenience we also introduce the variable

[TABLE]

which is zero at threshold. We renormalize the top quark mass in the on-shell scheme. Additionally we use $\mu_{r}$ and $\mu_{f}$ for the renormalization and factorization scales, respectively. Note that $\mu_{f}$ is only present in the collinear counterterm, where we explicitly show the dependence.

At LO, the application of the optical theorem leads to

[TABLE]

where $1/(2s)$ corresponds to the flux factor and the factors $N_{A}=N_{c}^{2}-1$ and $(2-2\epsilon)$ originate from averaging over gluon colours and helicities, respectively. We use the notation $\widetilde{\mathrm{Im}}$ to indicate that we compute the imaginary part only of cuts which involve two Higgs bosons.

Note that the factor $1/2$ for the identical Higgs bosons is contained in $\widetilde{\mathrm{Im}}(\mathcal{M}_{gg\rightarrow gg}^{(0)})$. In our calculation, for the sum over the gluon polarizations we use

[TABLE]

As a consequence we also have to consider amplitudes with external ghost particles, which have to be subtracted from the pure gluon contribution. For example, at NLO we have

[TABLE]

Note that all contributions with one external gluon and an external ghost or anti-ghost field are equal.

In a similar manner we obtain the partonic cross sections for the $qg$ (and $\bar{q}g$) and $q\bar{q}$ channels:

[TABLE]

Note that we do not need to consider ghost-quark scattering, since this only contributes starting from N3LO.

The $n_{h}^{3}$ contributions are obtained in analogy to Eqs. (6), (7) and (8) by replacing the LO part by the (virtual) NLO corrections proportional to $n_{h}^{3}$, see Fig. 1 (c), which plays the role of an effective LO contribution. In the following we denote this part by $\sigma^{(1),n_{h}^{3}}$. The NLO contributions in the above equations have to be replaced by the NNLO $n_{h}^{3}$ amplitudes (see Figs. 1 (f), (g) and (h)), which we denote by $\sigma^{(2),n_{h}^{3}}$. We can thus write

[TABLE]

with

[TABLE]

where the ellipses in Eq. (9) stand for N3LO terms. The virtual corrections $\sigma_{ij,\rm virt}^{(2),n_{h}^{3}}$ have been computed in Ref. [18] including terms up to $\rho^{2}$. Recently that calculation has been extended up to $\rho^{4}$ [29]. The main aim of this paper is the computation of the real corrections $\sigma_{ij,\rm real}^{(2),n_{h}^{3}}$ which we discuss later in subsection 2.3. In the next subsection we discuss the collinear counterterm $\sigma_{ij,\rm coll}^{(2),n_{h}^{3}}$.

2.2 Subtracting collinear divergences

The NNLO $n_{h}^{3}$ collinear counterterm is obtained from the convolution of the NLO cross section $\sigma_{gg}^{(1),n_{h}^{3}}$ and the gluon or quark splitting functions, which are given by

[TABLE]

with

[TABLE]

where $n_{l}$ is the number of massless quarks.

In the following we concentrate on the $gg$ channel. The calculations for the (anti) quark-induced channels proceed in a similar manner. The convolution integral is given by

[TABLE]

where the factor of 2 comes from the two external gluons. The integral over the delta distribution in $P^{(0)}_{gg}$ is trivial and in the parts without plus distributions we substitute $z=1-\delta(1-\mu)$ and subsequently expand in $\delta$. The integration over $\mu$ is then straightforward. For the contribution with the plus distribution we use the relation

[TABLE]

and again expand in $\delta$ after using $z=1-\delta(1-\mu)$.

In order to present result for the collinear counterterm we parametrize the NLO contribution $\sigma_{gg}^{(1),n_{h}^{3}}$ as

[TABLE]

where $c_{n}^{(0)},c_{n}^{(1)}$ depend on $\alpha_{s},G_{F},M_{t},m_{H}$ and $\mu_{r}$. We obtain for the $gg$ channel the following infinite series representation

[TABLE]

where

[TABLE]

Our result for the $qg$ channel reads

[TABLE]

To obtain terms to $\delta^{1/2+N}$ one should evaluate the series representations up to $n=N$, and then discard any incomplete higher-order terms which are produced.

In the ancillary file [30] we present expressions for $\sigma_{gg}^{(1),n_{h}^{3}}$ (from which the coefficients $c_{n}^{(0)}$ and $c_{n}^{(1)}$ can be extracted), $\sigma_{gg,\rm coll}^{(2),n_{h}^{3}}$ and $\sigma_{qg,\rm coll}^{(2),n_{h}^{3}}$ for arbitrary renormalization and factorization scale and expanded up to $\delta^{199/2}$ and $\rho^{4}$. In order to illustrate the structure of our result we provide some leading non-vanishing terms in the $\delta$ and $\rho$ expansion which are given by

[TABLE]

where $C_{A}=3$, $C_{F}=4/3$ and $T_{f}=1/2$ are colour factors. Furthermore, we have introduced the notation

[TABLE]

2.3 Workflow to compute $\sigma_{ij,\rm real}^{(2),n_{h}^{3}}$

To obtain the NNLO $n_{h}^{3}$ real-virtual contributions, we have to consider five-loop forward-scattering amplitudes with three closed top quark loops, each of which is coupled to one or two Higgs bosons. Since an exact calculation is currently not possible we perform an asymptotic expansion (see, e.g., Ref. [31]) for $M_{t}^{2}\gg m_{H}^{2},s$. As a result, we obtain products of three one-loop vacuum integrals with a two-loop integral with two or three massive Higgs boson propagators and forward scattering kinematics. For the latter we have to compute the imaginary part involving two Higgs bosons and a gluon or (anti-)quark.

The expansion in $1/M_{t}$ quickly develops a large number of terms. For this reason we pre-compute the $1/M_{t}$ expansion of the one-loop vacuum integrals with one or two external Higgs bosons and two or three external gluons. They are then inserted into the two-loop diagrams which are generated with effective Higgs-gluon vertices.

In the following we provide more technical details, the description applies to both NLO and NNLO contributions. We generate the one- and two-loop diagrams using qgraf [32] and select only the diagrams containing the relevant cuts using additional scripts. This output is then processed by q2e and exp [33, 34, 35], which generate FORM [36] code for the amplitudes and map them onto the corresponding integral families. We compute the colour factors of the diagrams using color [37]. The tadpole integrals of the “effective vertices” are evaluated using MATAD [38].

We initially define the one- and two-loop integral families without specifying forward-scattering kinematics, but rather with three independent (incoming) external momenta, $q_{1},q_{2},q_{3}$ and the relation $q_{4}=-q_{1}-q_{2}-q_{3}$. The identification $q_{2}=-q_{3}$ and $q_{1}=-q_{4}$ is applied at a later stage (see below).

We use this setup to obtain scalar expressions for each amplitude after averaging over the polarizations, spins and colours. Afterwards we decompose scalar products in the numerators in terms of denominator factors and map the scalar integrals onto 1 one-loop and 50 two-loop four-point integral families which are characterized by four and nine indices, respectively. In the two-loop case seven indices correspond to propagators and two to irreducible numerators, which we write as inverse propagators.

Next we impose the forward scattering kinematics, i.e. we set $q_{3}=-q_{2}$ and $q_{4}=-q_{1}$. This results in the propagators becoming linearly dependent. After partial fraction decomposition and identification of identical families we are left with 1 one-loop family of four propagators and 6 two-loop families of seven propagators. Graphical representations are given in Fig. 2. These integral families are suitable for IBP reduction which we perform with the help of FIRE5 [39]. We use FindRules, a built-in command of FIRE, to identify master integrals of different families and arrive at 2 one-loop and 16 two-loop master integrals which are depicted in Figs. 3 and 4.

The master integrals depend on $x=m_{H}^{2}/s$ where $0<x<1/4$. To obtain analytic results we use the powerful method of differential equations [40, 41, 42], which we derive with the help of LiteRed [43, 44]. To obtain the solutions we proceed in two ways. In the first approach we adopt the idea of Ref. [7] and expand around the threshold (i.e. $x=1/4$). We compute boundary conditions for all master integrals for $\delta\to 0$ and then use the differential equations to obtain for each master integral an expansion up to $\delta^{1019/2}$. More detail on this method can be found in Appendix A.

In the second approach we compute the master integrals without expanding in $\delta$ and express the result in terms of Goncharov polylogarithms [25]. The boundary conditions can be taken over from the first approach. We provide more details in Appendix B. As we show in Appendix B it is sufficient to use the (mathematically simpler) $\delta$-expanded results for phenomenological applications.

In the following we illustrate the structure of our results and provide the leading expansion terms in $\rho$ and $\delta$ for $\sigma_{ij,\rm real}^{(2),n_{h}^{3}}$:

[TABLE]

We have managed to perform the expansion up to order $\rho^{4}$ in the Feynman gauge. We additionally perform the expansion up to order $\rho^{2}$ for a general QCD gauge parameter $\xi$, and find that all dependence on $\xi$ drops out after summing the contributions of all bare diagrams. This provides a strong check of our calculation.

3 Results

We start by presenting analytic results for the partonic cross sections $\sigma_{gg}^{(2),n_{h}^{3}}$. Combining Eqs. (19) and (22) with the virtual corrections from [18, 29] the first two expansion terms in $\rho$ and $\delta$ are given by

[TABLE]

The results for $\sigma_{ij}^{(2),n_{h}^{3}}$ for $\rho=0$ have been obtained for the first time in Ref. [16] (denoted by $\hat{\sigma}^{b}$ in that paper). Note, however, that the final phase-space integration has been performed numerically and results are presented for hadronic quantities, and therefore a direct comparison with our results is non-trivial.

In Fig. 5 we show the partonic cross sections as a function of $\sqrt{s}$. In such situations, the exact LO contribution is often factored out in order to improve the behaviour at high energies. In Fig. 5 we refrain from doing so since we want to illustrate the convergence properties below the top quark pair threshold. For convenience we repeat the well-known LO and NLO results [7].

The $gg$ initiated NNLO contribution shows a similar pattern as at LO and NLO. Up to $\sqrt{s}\approx 330$ GeV a reasonable convergence is observed when including higher order terms in $\rho$. Beyond the top threshold we have no convergence. The $qg$ channel also shows good convergence up to the top quark threshold both at NLO and NNLO. No sign of convergence is observed for the $q\bar{q}$ channel. Note, however, that the contributions from production channels with quarks in the initial state are significantly smaller than the $gg$ channel.

At NLO the $n_{h}^{3}$ contribution is numerically much smaller than the remaining parts [8]. We thus expect that also at NNLO the $n_{h}^{2}$ terms will be numerically more important.

Let us mention that based on previous experience [7, 18], we did not expect better convergence behaviour than what is shown in Fig. 5. Nonetheless, the higher order terms in $\rho$ are important ingredients for the construction of approximations. For example, at NLO, the $\rho^{3}$ and $\rho^{4}$ terms help to obtain stable results for the cross section when combining large-$M_{t}$ results with information about the threshold behaviour, using Padé approximants [13].

4 Conclusions

We compute real radiation corrections to the cross section $gg\to HH$ by applying the optical theorem to Feynman diagrams with forward scattering kinematics. We concentrate on the subset of Feynman diagrams which involve three closed top quark loops, each of which couples to either one or two Higgs bosons, the so-called $n_{h}^{3}$ contribution. With the help of an asymptotic expansion for large top quark masses the five-loop diagrams factorize into three one-loop vacuum integrals and two-loop phase-space integrals which depend on $s$ and $m_{H}^{2}$. After IBP reducing the integral families we are left with 16 master integrals, which we compute both as an expansion around the threshold and as an exact expression in $s/m_{H}^{2}$. Although the radius of convergence is limited to a small region around threshold (i.e. $s\approx 4m_{H}^{2}$) our results are useful to provide information about the reliability of the effective-theory result. Furthermore, we expect that the power-suppressed top mass corrections will prove to be useful ingredients for approximation procedures and, last but not least, they are important benchmarks for future numerical calculations.

Work to compute the remaining ($n_{h}^{2}$) contributions to the NNLO real corrections to double Higgs boson production is ongoing, and will be published in a future paper.

Acknowledgements

F.H. thanks Max Delto for fruitful discussions and acknowledges the support of the DFG-funded Doctoral School KSETA. This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 — TRR 257 “Particle Physics Phenomenology after the Higgs Discovery”.

Appendix A Threshold expansion of master integrals

In this appendix we describe the calculation of the phase-space master integrals as an expansion around the threshold, i.e., for small values of $\delta$, see Eq. (4). For simplicity, we denote the Higgs boson mass by $m$ (instead of $m_{H}$).

All 16 master integrals (cf. Fig. 4) can be expressed in the form

[TABLE]

where the momenta $p_{3}$ and $p_{4}$ correspond to massive (Higgs) particles and the integration measures are given by

[TABLE]

with $E_{j}\equiv\sqrt{m^{2}+|\vec{p}_{j}|^{2}}$. The momentum $p_{5}$ corresponds to massless final-state parton and the integration measure reads

[TABLE]

The quantities $Q_{i}$ in Eq. (26) are given by

[TABLE]

We parametrise the $d$-dimensional momenta in Eq. (26) as

[TABLE]

and use the $\delta$ function in Eq. (26) to express the spacial components of $p_{4}$ as $\vec{p}_{4}=-\vec{p}_{3}-\vec{p}_{5}$. For the $d$-dimensional angular integrations we use (see, e.g., Ref. [45])

[TABLE]

We now consider Eq. (26), exploit the delta function and arrive at

[TABLE]

We can now perform the $\ell$-integration by noting that

[TABLE]

where $\cos\gamma=\cos\theta_{3}\cos\theta_{5}+\sin\theta_{3}\sin\theta_{5}\cos\phi_{5}$. Thus, we obtain

[TABLE]

where the upper limit of the $k$-integration has been determined by the $\delta$ function, and the factor $(1+[\ell+k\cos\gamma]/\sqrt{m^{2}+k^{2}+2k\ell\cos\gamma+\ell^{2}})^{-1}$ is the Jacobian of the $\delta$-function. Up to this point, the expression is exact. For convenience we also provide the propagators appearing in the $\mathcal{Q}_{i}$ in terms of $m,k,\ell,\delta$,

[TABLE]

At this point the expansion in $\delta$ is straightforward. Since $k\leq\sqrt{s\delta}/2=m\sqrt{\delta/(1-\delta)}$ and $\ell_{\delta}\leq m\delta$, which follow from Eqs. (55) and (54) respectively, we can expand in both variables. In practice we proceed as follows: after substituting Eq. (54) into Eq. (55) and making a change of integration variable $k=\xi\sqrt{s\delta}/2$ the resulting integrand is a polynomial in all integration variables $\xi,\cos\theta_{3},\cos\theta_{5}$ and $\cos\phi$.

We now show, for each master integral, three $\delta$-expansion terms for the leading coefficient in the $\epsilon$ expansion. Our results read

[TABLE]

where

[TABLE]

In principle, one can compute the series expansion of $I_{i}$ up to arbitrary order in $\delta$ using the method described above. However, the number of terms in the integrand grows rapidly, and we stopped the computation at $\mathcal{O}(\delta^{11/2})$. A more efficient approach to obtain high-order terms in $\delta$ is based on differential equations with respect to $\delta$. It turns out that for the boundary conditions required to solve the differential equations, the leading-order term of each integral is sufficient. We can compute this term without expanding in $\epsilon$ and obtain

[TABLE]

We managed to obtain all master integrals up to $\mathcal{O}(\delta^{219/2})$ within a few hours of CPU time which is sufficient for all practical purposes. Note, however, that with this approach the computation of further orders in the expansion is not difficult; we have produced expressions for the expansion up to $\mathcal{O}(\delta^{1019/2})$.

It turns out that all of the master integrals have non-integer powers of $\delta$ after expansion. An overall factor $\delta^{1/2}$ comes from the measure of the $k$ integration in Eq. (55). Terms with odd powers of $k$ in the integrand are candidates for terms with integer power of $\delta$. However, these vanish after the angular integration. Since the differential equation relates terms whose powers of $\delta$ differ by integers, we have to confirm the absence of integer powers in $\delta$ by an explicit calculation of the first few terms, as we have shown above.

Appendix B Exact master integrals

B.1 Two-loop master integrals

In order to obtain exact results for the master integrals we first transform the differential equation into a canonical form [46]. Afterwards we perform the integrations order by order in $\epsilon$ and express the analytic results in terms of Goncharov polylogarithms [25] which can be evaluated numerically using GiNaC [47, 48].

To obtain a canoncial form, we apply Lee’s algorithm [49] as implemented in the program Epsilon [50] to bring the $16\times 16$ system to normal Fuchsian form (i.e. $\epsilon$ does not yet factorize) and observe that the differential equations contain poles at $x=\{0,1/4,1,r_{1}=\exp(i\pi/3),r_{2}=\exp(-i\pi/3),-1/3\}$. The first and second poles correspond to the limits $s\rightarrow\infty$ and $s\rightarrow 4m_{H}^{2}$, respectively. The remaining four poles are only present in sectors with a third, uncut Higgs propagator, i.e. for integrals $I_{7}$ to $I_{16}$. Note that $x=1$ corresponds to the threshold for single Higgs production.

Let us first have a closer look at the sub-matrix which corresponds to $I_{1}$ and $I_{2}$. After Fuchsification the matrix residue at $x=1/4$ for the subsystem reads

[TABLE]

which indicates that one of the eigenvalues is half-integer for $\epsilon\rightarrow 0$. This implies the appearance of the square root $\sqrt{1-4x}$ in the alphabet, which is due to the two-particle branch cut. We rationalize this root by introducing a new variable $y$, defined as

[TABLE]

After rationalization we managed to find a canonical basis for the first 14 master integrals.

The homogeneous part of the differential equation for $I_{15}$ contains another residue with half-integer eigenvalue, namely $1/(y^{4}-3y^{3}+5y^{2}-3y+1)$. After re-scaling $I_{15}$ we are able to bring the differential equations into a form in which we can also factor out $\epsilon$ and the homogeneous parts of $I_{15}$ and $I_{16}$ only contain single poles in $y$. The root $\sqrt{1-6y+7y^{2}-6y^{3}+y^{4}}$ appears now in the inhomogeneous contributions. Note, however, that only the leading terms in $\epsilon$ of these two integrals are needed, which means that we do not have to iteratively integrate over the square root.

To simplify the integration of the differential equations and the manipulation of the amplitudes, we do not partial fraction the polynomials

[TABLE]

which appear in the denominators of the differential equations. Rather, we use the following integration kernels

[TABLE]

in our final expressions, where $n=1,2$ and $k=1,\ldots,4$. Note that for the numerical evaluation we can partial fraction the quadratic and quartic kernels and rewrite the integrals as a sum of Goncharov polylogarithms, i.e.

[TABLE]

where the $r_{i}$ are the roots of $P_{2}$ and the $s_{i}$ are the roots of $P_{4}$. While all iterated integrals over the kernels in Eq. (63) are real-valued for $-1\leq y\leq 0$, the individual Goncharov polylogarithms on the r.h.s. of Eq. (64) are not.

We expand all master integrals in $\epsilon$ up to the order needed for the finite NNLO part, which means that some master integrals are expanded up to the $\epsilon^{2}$ and for others only the $1/\epsilon$ pole is required. Let us mention that the $\mathcal{O}\left(\epsilon^{2}\right)$ terms of the master integrals $I_{9}$ and $I_{12}$ contain Goncharov polylogarithms up to weight 4. However, only the difference $I_{9}-I_{12}$, which only contains weight-3 Goncharov polylogarithms, is needed up to $\mathcal{O}\left(\epsilon^{2}\right)$. The sum $I_{9}+I_{12}$ is only needed up to the linear $\epsilon$ term and contains at most weight-3 terms. Therefore the cross section can be expressed in terms of Goncharov polylogarithms up to weight 3.

In the ancillary file [30] we provide analytic expressions for all masters integrals, both expanded in $\delta$ up to order $\delta^{219/2}$ and expressed in terms of Goncharov polylogarithms. For the latter we do not give separate expressions for $I_{9}$ and $I_{12}$ but for the combinations $I_{9}-I_{12}$ and $I_{9}+I_{12}$.

We are now in a position to compare the exact results for the master integrals with the $\delta$-expanded expressions. In Fig. 6 we show, for two typical examples, the exact result (solid curve) and results expanded up to various orders in $\delta$ (dashed curves). We plot the expressions as a function of $\delta$ but suppress the threshold region where very good agreement is found. The expressions expanded up to $\delta^{10}$ start to deviate from the exact curve above $\delta\approx 0.8$. Agreement up to $\delta\approx 0.9$ is observed if 20 expansion terms in $\delta$ are included and for 100 terms we reach $\delta\approx 0.97$. Note that $\delta=0.9$ corresponds to $\sqrt{s}\approx 800$ GeV where the parton distribution functions are already quite small. For many applications it is therefore sufficient to work with the $\delta$-expanded expressions.

During preparation of this manuscript, the Mathematica package PolyLogTools [51] was made available. We were able to use it to expand the exact expressions for our master integrals in $\delta$ to order $\delta^{11/2}$, and found full agreement with our expansions of Appendix A.

B.2 One-loop master integrals

The two one-loop master integrals (see Fig. 3) have been computed in Ref. [7] as an expansion in $\delta$. We have solved the system of differential equations using the same approach as at two loops (see above) and obtained the following results which are exact in $y$ (see Eq. (61)):

[TABLE]

where ${\cal N}$ is given in Eq. (58).