Higgs boson pair production: top quark mass effects at NLO and NNLO

We compute next-to-next-to-leading order QCD corrections to the gluon-induced production cross section of Higgs boson pairs in the large top quark mass limit using the soft-virtual approximation. In the limit of infinitely-heavy top quark we confirm the results in the literature. We add two more expansion terms in the inverse top quark mass to the $M_t\to\infty$ result. Since the $1/M_t$ expansion converges poorly, we try to improve on it by factorizing the exact leading order cross section. We discuss two ways of doing that and conclude that the finite top quark mass effects shift the cross section at most by about 10\% at next-to-leading order and by about 5\% at next-to-next-to-leading order.


Introduction
In the coming years, one of the main tasks in particle physics is the understanding of the mechanism of the electroweak symmetry breaking. After the experimental determination of the Higgs boson mass, the Higgs potential is fully fixed in the Standard Model. However, it is very important to independently measure the self-coupling of the Higgs boson, which can be obtained from studying the production of Higgs boson pairs. Since the corresponding cross section is O(10 3 ) times smaller than the one for single Higgs boson production, Higgs boson pair production poses a challenging problem for the LHC, even after the luminosity upgrade around 2020.
There are a number of phenomenological analyses which investigate the possibility to extract the self coupling from cross section measurements. First studies have been performed more than 15 years ago [1][2][3]. Since the discovery of the Higgs boson there has been an increasing interest in this topic and a number of refined analyses have been performed, see, e.g., Refs. [4][5][6][7][8].
Higgs boson pairs can be produced via the fusion of two partons or vector bosons, via the radiation off vector bosons, or in association with heavy quarks. Similar to single Higgs boson production, the numerically dominant mechanism is gluon fusion although the leading order (LO) contribution is loop-suppressed. Due to the larger Yukawa coupling, the dominant contribution comes from top quark loops in the Standard Model. For this reason we concentrate in this paper on such contributions.
For the LO order corrections the exact dependence on the top quark mass and the centerof-mass energy is known [9,10]. At next-to-leading (NLO) QCD corrections have been computed for the first time more than 15 years ago [11,12] in the infinite top quark mass limit using an effective theory. Finite top quark mass effects have been investigated in Ref. [12] where a systematic expansion in the inverse top quark mass has been applied and a quantitative estimate of the quark mass effects has been provided. It has been estimated that they do not exceed O(10%) of the NLO contribution. Finite top quark mass effects have also been considered in Ref. [13] where the exact real radiation contribution is combined with the effective-theory virtual corrections. As a result, a reduction of about −10% of the cross section is obtained. We will comment in Section 3 on this issue.
Within the effective theory also next-to-next-to-leading (NNLO) contributions are available [14,15]. In this context it is interesting to note that the three-loop matching coefficient of the effective operator for two Higgs bosons and two, three or four gluons is different from the one for single Higgs boson production [16]. The results for the virtual corrections obtained in Ref. [14] have been cross-checked in Ref. [16] where the calculation has been performed without reference to the effective theory. The resummation of thresholdenhanced logarithms to next-to-next-to-leading logarithmic (NNLL) accuracy has been performed in Refs. [17,18].
The remainder of this paper is organized as follows: In the next Section we review the construction of the soft-virtual approximation for the production cross section and discuss two options to factorize the exact LO result from the higher order contributions. We argue that a factorization at the level of the differential cross section w.r.t. the Higgs boson pair invariant mass leads to more stable results. Afterwards we reconsider in Section 3 the top mass corrections at NLO. Virtual NNLO corrections including finite top quark mass effects are computed in Section 4. They are used in Section 5 to present phenomenological results for Higgs pair production up to NNLO. Section 6 contains our conclusions.

Factorizing the exact LO expression
We write the perturbative expansion of the partonic cross section for the production of Higgs boson pairs in the form where α s ≡ α (5) s (µ), √ s is the partonic center-of-mass energy and ij ∈ {gg, qg,qg, qq}.
Since the quark-induced channels are numerically small [11] we consider in this paper only the gg channel. We use the variable to parametrize the dependence of the cross section on the Higgs boson and top quark mass. We renormalize the top quark mass in the on-shell scheme. Furthermore, we set the factorization and renormalization scale equal to each other and write µ = µ r = µ f .
For later convenience we introduce for σ gg→HH+X and denote the sum of the first two terms on the right-hand side by σ NLO = σ LO + δσ NLO .
Finite top quark mass effects to gg → HH at NLO have been considered for the first time in Ref. [12]. The applied method is based on "reversed unitarity" [19] where, with the help of the optical theorem, the imaginary part of forward scattering amplitudes are computed to obtain the total cross section. The virtual corrections have also been computed by directly considering the gg → HH amplitude. Expansion terms up to order 1/M 12 t of the NLO contribution to σ(pp → HH) have been computed [12,20,21]. The factorization of the exact LO corrections has then been implemented at the partonic level for the total cross section using σ (1) where σ (0) gg,exact contains the exact dependence on s and ρ. In Ref. [11], where the NLO result has been computed for the first time in the heavy top quark limit, a different approach has been applied. Actually, the exact LO result has been factorized before the integration over the Higgs pair invariant mass. In this approach the (total) NLO cross section reads where 'exact' and 'exp' reminds whether exact or (in ρ) expanded results are used. Q 2 is the invariant mass squared of the Higgs boson pair. Expressed in terms of dσ/dQ 2 it is possible to re-write Eq. (4) as From the comparison of Eqs. (5) and (6) one expects that (5) leads to better results since the differential factorization (DF) in Eq. (5) results in a better-behaved integrand, in particular for Q 2 > 4M 2 t . For the virtual corrections, which are proportional to δ(s − Q 2 ), one has immediate access to the Q 2 dependence and the DF of Eq. (5) can be applied. 1 The real corrections, however, are obtained with the help of the optical theorem which directly leads to the total cross section and thus Eq. (5) can not be used. For the construction of the softvirtual (SV) approximation, which is discussed below, we need in addition to the virtual term also the contributions from soft gluon emission. Since the soft contributions are proportional to the LO Born cross section Eq. (5) can immediately be applied to the SV cross section.
In the following we discuss the construction of the SV approximation for the cross section (see also Ref. [22]). For simplicity we consider for the following schematic reasoning the total cross section σ. Note, however, that the arguments also hold for dσ/dQ 2 . In a first step we split σ according to where the two terms on the right-hand side correspond to the virtual corrections (including ultra-violet counterterms) and the real corrections (including the contributions from the factorization of initial-state singularities). They are individually divergent and only the sum is finite. In the next step we re-organize the terms on the right-hand side such that σ can be written as where "SV" refers to "soft-virtual" and "H" to "hard". This is achieved by splitting σ virt+ren into a divergent and a finite term and separating σ real+split into a (divergent) soft and (finite) hard contribution. The soft contribution is combined with σ virt+ren to obtain Σ SV such that Σ SV and Σ H are separately finite. Σ div is constructed following Ref. [23]; explicit expressions can be found in Refs. [16,22]. Note that the finite part is constructed solving the equation for Σ fin . In this paper σ virt+ren is computed including top quark mass effects. Mass effects are automatically taken into account in Σ div and Σ real+split soft since these contributions are proportional to the exact LO cross section [22]. Note that Eqs. (7), (8), (9) and (10) also hold for the differential cross section dσ/dQ 2 and thus a factorization as suggested in Eq. (5) can be performed for the soft-virtual contribution.
We adopt the notation from Ref. [22] and parametrize radiative corrections to the partonic cross section via and where the renormalization scale dependence in the strong coupling constant α s is suppressed. From Eq. (12) one obtains for the total cross section with In the second line of Eq. (14) we split G(z) into soft-virtual and hard contribution. Note that in our approach we do not have access to G H (z). Actually, at NNLO we only have G SV (z) at hand and at NLO only the heavy top expansion of Explicit results for G SV (z) can be found in Ref. [22] including higher order terms in ǫ specifying, however, the renormalization and factorization scale to √ s ≃ Q 2 . For completeness we provide the NLO and NNLO results for G SV (z) for generic µ in the limit ǫ → 0. The results read where C F = 4/3, C A = 3 are QCD colour factors, n l = 5 is the number of massless quarks, L = log(µ 2 /s) and The quantities σ (1) fin and σ (2) fin are evaluated for µ = √ s. They are obtained from dσ (1) v,fin /dt and dσ (2) v,fin /dt in Eq. (11) after integration over t.
In the next Section we apply the formalism described in this Section at NLO and compare to the results obtained in Ref. [12]. The numerical effects at NNLO are presented in Section 5 using the mass corrections computed in Section 4.

Revisiting NLO
In this Section we restrict ourselves to NLO and compare the results of Ref. [12] with the alternative factorization discussed in the previous section.
To obtain numerical results we employ MSTW2008 parton distribution functions (PDFs) [24] and consistently use N k LO PDFs to compute N k LO (k = 0, 1, 2) cross sections. We assume the energy of the LHC to be 14 TeV and adopt the values of the strong coupling constants α s (M Z ) that we use in our computation from the MSTW PDF fit: We renormalize the top quark in the on-shell scheme and use M t = 173.21 GeV [25].
For the Higgs boson we use m H = 125.09 GeV [26]. For numerical results shown in this Section the renormalization and factorization scales have been set to 2m H .
In Fig. 1 we show the results for the partonic cross section as a function of √ s from terms. From Fig. 1 it is evident that the hard contribution is numerically much smaller than the soft-virtual one.
In Fig. 2 we compare the solid curves from Fig. 1 (which are dotted in this plot) with the results obtained with the help of DF applied to the SV approximation and adding the hard contribution as given by the dashed lines in Fig. 1. The relative position of the lines and the colour coding is as in Fig. 1  to comparable results. However, the DF curves have their maximum at lower values of √ s and lead to a smaller cross section for larger values of √ s. Furthermore, one observes a drastic improvement in the convergence behaviour when including higher order mass corrections. In particular, for √ s = 400 GeV the difference between the infinite top mass result and the one including 1/M 12 t terms amounts to only ≈ 0.05 fb to be compared with ≈ 0.25 fb for the dashed curves.
At this point we want to stress that the splitting between the soft-virtual and hard contribution in Eq. (8) is arbitrary. In fact, the soft-virtual contribution of G(z), G SV (z), is constructed for the limit z → 1 and thus it is possible to replace G SV (z) by f (z)G SV (z) with f (1) = 1 (see, e.g., discussions in Refs. [27,28]). The hard contribution is modified accordingly such that the sum of Σ SV + Σ H does not change. One observes that for f (z) = z the soft-virtual contribution approximates the partonic NLO contribution to the cross section very well 2 such that at the hadronic level the deviations are below 2%. Based on this observation we use f (z) = z for the NNLO cross section where we only have the soft-virtual approximation at hand. Furthermore, better results are obtained by replacing L in Eq. (17) by L = log(µ 2 /Q 2 ) which is justified since √ s ≈ Q 2 in the soft limit. Note that in Ref. [15] it has been observed that the soft-virtual approximation constructed in Mellin space approximates the full (effective-theory) result with an accuracy of 2%.
It is interesting to look at the partonic K factor which is defined via Results for the two methods to factorize the exact LO term are plotted in Fig. 3 as a function of √ s where the dashed curves are already shown in Ref. [12]. One observes that DF leads to a lower K factor and that the spread among the various ρ orders is smaller. Furthermore, it is interesting to note that for DF the top quark pair threshold behaviour of the LO term is not washed out in contrast to the dashed curves. It is common to both factorization methods that there is a strong raise when approaching the threshold for Higgs boson pair production (see also discussion in Ref. [12]).  where the luminosity function is given by f g (x) are the gluon distribution functions in the MS scheme. Note that in the soft limit √ s cut is a good approximation to Q 2 . The various lines in Fig. 4 correspond to the inclusion of different orders in ρ at NLO. For convenience we show on the right end of Fig. 4 the total cross section for √ s H = 14 TeV. Note that the approximation used for the computation of the ρ n terms is not valid for large values of √ s cut (neither is the effective-theory result). However, it can be used as an estimate of the mass correction terms. Using the spread as an estimate for the uncertainty we conclude that a finite top mass induces a ±10% uncertainty on top of the infinite top quark mass result.
The lower panel of Fig. 4 shows the hadronic K factor which is obtaind from Eq. (20) by replacing σ by σ H and using NLO PDFs in the numerator and LO PDFs in the denominator. K NLO raises close to threshold, however, for √ s cut ∼ > 500 GeV one observes a flat behaviour of K NLO ≈ 1.6 (for µ = 2m H ).
Top quark mass effects to double Higgs boson production have also been considered in Ref. [13]. In the approximation used in that reference the real corrections are treated exactly, however, the infinite top quark mass approximation is used for the virtual corrections. A decrease of the cross section by about 10% due to finite top quark mass is reported.
In Fig. 5 we split 3 the partonic results for the solid lines of Fig. 1 (which corresponds to the factorization according to Eq. (4)) into virtual corrections (including the ultra-violet counterterms; dotted lines) and the real-radiation part which include the contributions from the splitting functions (dashed lines).
We observe that for √ s ∼ < 400 GeV, the region where our approximation is valid, the

Top quark mass corrections at NNLO
In this Section we compute the virtual corrections to the NNLO cross section including top quark mass corrections. Afterwards we construct Σ fin as described in Section 2 and use Eq. (17) to evaluate the partonic and hadronic cross sections including 1/M 4 t correction terms.

Calculation
We have applied two methods to compute the virtual corrections. In the first one we consider the amplitudes for gg → HH up to three-loop order and in the second one the forward scattering amplitude is considered, which, after taking the imaginary part, directly leads to the total cross section. In both cases we perform an asymptotic expansion for large top quark mass. The first approach has the advantage that it is straightforward to introduce Higgs boson decays whereas the second approach can immediately be applied to real corrections.   Fig. 6. They are generated with the help of qgraf [29]. Note that in this approach no contributions with external ghosts have to be considered since we project on physical states. The transformation to FORM [30] code is done with the program q2e [31,32] and the asymptotic expansion for large top quark masses is realized with the help of exp [31,32]. After expanding the identified hard subgraphs in the small quantities one arrives at one-scale vacuum integrals up to three loops and massless one-and two-loop four-point diagrams with two massless and two massive external momenta. The vacuum integrals are computed using MATAD [33]. In the case of the four-point integrals we apply FIRE [34,35] to express them as linear combinations of master integrals. The latter are shown in Fig. 7 where we have q 2 1 = q 2 2 = 0, and q 2 3 = q 4 4 = m 2 H . Analytic results for all master integrals can be found in the literature (see, e.g., Refs. [36][37][38]). In this paper we only show results for the triangle graph in the second line of Fig. 7 since for our purpose the representations given in Refs. [36,38] are less suited. We use instead with G o ({w i }; z) are Goncharov Polylogarithms [39] with weight {w i } and argument z defined through with w i , x ∈ C and The functions of the ǫ 0 term in Eq. (23) can be expressed in terms of logarithms and dilogarithms via We have cross checked the numerical result for the ǫ 0 and ǫ 1 terms of I 1 (4) in Eq. (23) against FIESTA [40].
Using this method we have computed terms up to order 1/M 12 t at NLO [12,20,21] and terms up to order 1/M 4 t at NNLO. As an important check we have computed the 1/M 2 t corrections for general QCD gauge parameter which drops out in the final expression.
From the calculation of gg → HH one obtains in a first step results for dσ/dt. Integration over phase space then leads to dσ/dQ 2 . For the results in Section 5 this integration is performed numerically. Figure 7: One-and two-loop master integrals needed after applying asymptotic expansion to the amplitude gg → HH. All internal lines are massless, q 2 1 = q 2 2 = 0, and q 2 3 = q 4 4 = m 2 H .

Amplitude gg → gg
The second method is based on the use of the optical theorem in analogy to the NLO calculation performed in Ref. [12]. This method serves as an important cross check. In the following we provide some of the technical details • The amplitudes for gg → gg are generated with the help of qgraf [29].
• In a first step about 17 million diagrams are generated. However, most of them do not contain a cut through exactly two Higgs bosons. For this reason we post-process the qgraf output and filter [21,41] the amplitudes describing the virtual corrections to gg → HH. Typical Feynman diagrams are shown in Fig. 8.
• Our in-house FORM code applies projectors (−g µν ) for each pair of external gluons which includes also non-physical degrees of freedom in the sum. Thus also contributions with ghosts in the initial state have to be considered. Note that this is in contrast to single Higgs boson production which has no virtual contributions with ghosts in the initial state. The application of the asymptotic expansion for large top quark mass leads to a factorization of the five-loop integrals into the following contributions: The mass scale in the vacuum integrals is given by the top quark. They are again computed with the help of MATAD [33]. For the remaining integrals, which depend on δ, we use the inhouse programs rows [41] and TopoID [21,41] to perform the reduction to master integrals. The latter are depicted in Fig. 10.
All three-loop master integrals factorize into a two-loop form factor contribution and the one-loop master integral of the LO calculation. From the latter only the imaginary part is needed which is well known. The results for the two-loop form factor integrals can, e.g., be found in Refs. [43].
The two-loop master integrals are more involved. A numerical calculation would probably be possible, however, we follow the approach outlined in Ref. [12] for the NLO master integrals and perform an expansion in δ = 1 − 4m 2 H /s [cf. Eq. (16)]. All integrals contain a massless sub-loop for which analytic expressions are known. The massless two-point function can be expressed in terms of Γ functions and results for the triangle with two massless external legs and the (crossed) box can be found in Ref. [37]. Analytic expressions for the triangle diagram with squared external momenta s, m 2 H , m 2 H are given in Eq. (23). After expanding in δ the remaining phase-space integration can be performed analytically. We have computed expansion terms up to order δ 10 and found agreement with the results obtained in the previous subsection which for this purpose have also been integrated analytically after performing the expansion in δ.
The approach based on the optical theorem requires special care in the treatment of the imaginary parts originating from the two-loop form factor diagrams. Such contributions either correspond to |M NLO | 2 or (M LO M ⋆NNLO + M NNLO M ⋆LO ). In the former case the two-loop integrals originate from the product of two one-loop contributions containing factors (−1+i0) ǫ and (−1−i0) ǫ , respectively, which finally leads to (−1+i0) ǫ (−1−i0) ǫ = 1.
In the other case one has (−1 . The corresponding discussion for single Higgs production can be found in Ref. [44]. Note that the approach based on the gg → HH amplitudes leads to simpler intermediate expressions. Thus, it is possible to allow for a general QCD gauge parameter ξ when computing the 1/M 2 t terms. Furthermore, also the 1/M 4 t corrections could be evaluated (for ξ = 0) whereas in the optical theorem approach only 1/M 2 t terms could be computed in Feynman gauge. However, let us mention that this approach can be used in a straightforward way to compute NNLO top quark mass effects to the hard contribution of the total cross section whereas in the approach of the previous subsection this is less obvious.
To conclude this Section let us summarize our procedure to obtain the SV corrections at NNLO. We compute virtual corrections to gg → HH including 1/M 4 t corrections. Note that we have dσ/dQ 2 | virt ∼ δ(Q 2 − s). Using Eq. (10) we construct σ fin which enters G SV in Eq. (17). The differential and total cross section is then obtained with the help of Eqs. (12) and (15).

Improving NNLO
In this Section we discuss the effect of the 1/M 2 t and 1/M 4 t terms on the NNLO cross section for the production of Higgs boson pairs. The cross section in the infinite top quark mass limit has been computed in Ref. [14]. In Ref. [16] the three-loop matching coefficient has been added, completing the NNLO prediction, and the virtual corrections from Ref. [14] have been cross checked.
The results for the virtual corrections computed in Section 4 are inserted in the formalism described in Section 2 to construct the quantity σ fin which enters Eq. (17). The result for the partonic cross section is shown in Fig. 11 as a function of the partonic center-ofmass energy where the exact LO result is compared with NLO and NNLO. At NLO and NNLO three terms in the mass expansion are shown. 4 Furthermore, at NNLO the SV approximation is shown for f (z) = z (cf. discussion in Section 3). Note that the NNLO curves peak for smaller values of √ s than at NLO and LO. As far as the top quark mass corrections are concerned the same pattern is observed as at NLO: the correction term of order ρ decreases the infinite top quark mass result which is overcompensated by the ρ 2 term resulting in a small positive correction.
In Fig. 12 we compare the NNLO-SV contribution for two different scales, µ = 2m H and µ = Q 2 . This plot furthermore shows the effect of f (z) = 1 and f (z) = z. Note that the choice f (z) = z, which we expect to better approximate the complete result, leads to an increase of the cross section.
The hadronic cross section as a function of √ s cut is shown in Fig. 13 for µ = 2m H . At LO the exact result is used and the NLO curve is based on the ρ 0 results. Using instead the ρ 6 terms leads to an upwards shift of about 5% as can be seen in Fig. 4. At NNLO three curves are shown which include terms up to ρ 0 , ρ 1 and ρ 2 . As at NLO one observes good convergence up to √ s cut ≈ 400 GeV. For higher values of √ s cut the ρ 1 curve is below and   Figure  The results in the right panel with "∞" at the bottom correspond to the prediction of the total cross section. For this plot µ = 2m H has been used. The two known mass correction terms lead to a change of the cross section by about ±2%. Assuming a similar pattern as at NLO we thus estimate that NNLO top quark mass corrections change the effective-theory result by at most ±5%.

Conclusions
We compute NLO and NNLO corrections to double Higgs boson production in gluon fusion beyond the effective-theory approach. The starting point of the calculation are full-theory Feynman diagrams. We perform an asymptotic expansion in the limit where the top quark mass is large and compute at NNLO three terms in the 1/M t expansion for the virtual corrections. They are used to construct a soft-virtual approximation for the production cross section. In the limit M t → ∞ the effective-theory result of Ref. [14] is confirmed and 1/M 2 t and 1/M 4 t terms are added. The main result of this paper is illustrated in Fig. 13 where the hadronic cross section is shown as a function of √ s cut (a technical cut on the partonic center-of-mass energy). The curves including 1/M t corrections deviate from the infinite mass result only by a few per cent which leads us to the estimate that the effective-theory result is accurate to ±5%. Analog results for the mass corrections at NLO are shown in Fig. 4 which constitutes an update of Ref. [12]. Here we estimate the uncertainty to ±10%.
We want to stress that the results obtained in this paper provide excellent approximations for small values of the partonic center-of-mass energy, say below √ s ≈ 400 GeV. Although in this region the cross section is small it is of interest since there the cross section has a characteristic behaviour. Furthermore, it is possible to use our result in this region as a benchmark for future calculations taking into account the exact dependence on M t .
The methods described in Section 4 can also be used to compute top mass corrections to the real radiation part. However, the simplifications used in Ref. [15] where results have been obtained for M t → ∞ do not apply once finite mass effects are considered. The calculation is much more challenging since significantly more Feynman diagrams contribute and more complicated master integrals have to be computed.
In this paper for the first time the effect of a finite top quark mass has been examined for the NNLO cross section for double Higgs boson production. Whereas at NLO an exact calculation is within reach this is certainly not the case at NNLO. Thus our results become particular important once our NLO approximations are compared to an exact calculation which increases the confidence in the uncertainty estimate. Furthermore, one probably can obtain a prescription to tune the approximation procedure and hence reduce the uncertainty at NNLO.