Two-loop virtual corrections to Higgs pair production

We present the two-loop virtual corrections to Standard Model Higgs boson pair production via gluon fusion $gg\to HH$ in the heavy top quark limit. Based on this result, we evaluate the corresponding cross section at the LHC at 14 TeV in the next-to-next-to-leading order soft-virtual approximation. We find an inclusive K-factor of about 2.4, resulting in an increase close to 23% with respect to the previous available calculation at next-to-leading order. As expected, we observe a considerable reduction in the renormalization and factorization scale dependence.


Introduction
Recently, both ATLAS and CMS collaborations have discovered a new boson with a mass around 125 GeV [1,2] at the Large Hadron Collider (LHC). Its properties are, so far, compatible with the long sought Standard Model (SM) Higgs boson [3]. In order to decide whether this particle is indeed responsible for the Electroweak Symmetry Breaking (EWSB), it is crucial to measure its couplings to fermions and gauge bosons and to verify their proportionality to the particle masses. Furthermore, a precise measurement of the Higgs self-interaction is needed. The measurement of the Higgs self-couplings is the only way to reconstruct the scalar potential. After EWSB, the Higgs potential takes the form In the SM the trilinear and quartic self-couplings take the same value, λ = λ ′ = M 2 H /(2v 2 ), where v ≃ 246 GeV is the Higgs vacuum expectation value and M H its mass. In most new physics scenarios these couplings deviate from the SM values. Therefore, a determination of the Higgs self-interaction is necessary both to understand the EWSB mechanism and to try to distinguish the SM from other models.
The Higgs quartic coupling can be in principle studied via triple Higgs boson production. However, this cross section is too small to be measured at the LHC [4], and then a determination of its value is not possible at present time. The situation is different for the trilinear coupling λ via Higgs pair production if very high luminosities can be achieved, The possibilities of observing Higgs pair production at the LHC have been discussed in Refs. [5][6][7][8][9][10][11][12]. Though the analysis is challenging due to the smallness of the signal cross section and the large QCD background, it has been shown to be achievable at a luminosity-upgraded LHC. For example for bbγγ and bbτ + τ − final states, after the application of proper cuts, the significances obtained are ∼ 16 and ∼ 9 respectively, for √ s H = 14 TeV and L = 3000 fb −1 [8]. These are so far the most promising final states for the Higgs trilinear coupling analysis. The application of jet substructure techniques was shown to be important to further improve on the sensitivity of the discovery channels [6,7,13].
As it occurs for single Higgs [14], the dominant mechanism for SM Higgs pair production at hadron colliders is gluon-gluon fusion, mediated by a heavy-quark (mainly top) loop. The corresponding cross section has been calculated at leading-order (LO) in Refs. [15][16][17]. The nextto-leading order (NLO) QCD corrections have been evaluated in Ref. [18] in the large top-mass approximation and found to be rather large, with an inclusive K-factor close to 2, a very similar situation to the one observed for single-Higgs production at the same order [19][20][21]. Considering that the next-to-next-to-leading order (NNLO) corrections for single-Higgs are also sizable [22][23][24], it becomes essential to reach the same accuracy for double-Higgs production in order to provide precise predictions for the process.
A full NNLO calculation requires the evaluation of the corresponding amplitudes for double real radiation, real emission from one-loop corrections and the pure virtual two-loop contribution. In this article we present the explicit results for two-loop virtual corrections to the partonic process gg → HH in the heavy top quark limit. Furthermore, we combine these results with the universal formula presented in Ref. [25] to obtain the NNLO soft-virtual approximation to the cross section, as a first step towards a full NNLO calculation, and present numerical results for the cross section expected at the LHC within that approximation.

Two-loop virtual corrections
We present here our results on the two-loop corrections. In order to simplify the presentation we directly provide the contribution of two-loop diagrams to the corresponding partonic cross section. As usual, divergences are dealt with by using dimensional regularization with n = 4 − 2ǫ dimensions, and we use the MS renormalization scheme.
As it was mentioned before, we strictly work within the heavy top quark approximation, where the single and double-Higgs coupling to gluons is given by the effective Lagrangian where G µν represents the gluonic field strength tensor. In order to obtain the NNLO cross section for gg → HH, we need the coefficients C H and C HH up to O(α 3 S ). The first one takes the following form [26,27]: where M t is the on-shell top quark mass, µ R is the renormalization scale and N f is the number of light flavors. The coefficient C HH is known up to O(α 2 S ) [20], and coincides to that order to C H . We will write For the phenomenological results, we will assume C (2) H , where the latter is defined by the squared bracket in Eq.(3).
In Figure 1 we show a sample of the Feynman diagrams needed for the calculation, and we introduce the notation for each contribution. Since the structure of gHH and ggHH vertices is the same, the loop corrections to both Born level diagrams are proportional to the gluon form factor, and those contributions are labeled as FF(1) and FF (2). To compute these amplitudes, we rely on the two-loop gluon form factors presented in [28][29][30][31]. On the other hand, the contributions arising from diagrams with tree-level two gHH vertices (labeled as 2V (1)) and the corresponding one-loop correction to them (labeled as 2V (2)), which are of the same order in powers of the strong coupling constant as the form factor-like corrections FF(1) and FF(2) have more complex kinematics and require an explicit computation.
The NNLO virtual corrections (at the level of squared amplitudes) include the interference between FF(2)+2V(2) and LO diagrams, the squares of FF(1) and 2V (1) and their corresponding interference. The calculation was performed using the Mathematica packages FeynArts [32] and FeynCalc [33] for the generation of the diagrams and the manipulation of amplitudes, and the algorithm FIRE [34] to reduce the resulting expressions into master integrals, which are obtained from Ref. [35].
The partonic virtual corrections σ v to the cross section are obtained by integrating the squared amplitudes over the Higgs pair phase space, that is where we also include the flux factor, the average over helicities and colors of the incoming gluons and the factor for identical particles in the final state. Expanding in powers of the strong coupling α S : The renormalized NLO virtual contribution σ (1) is given by while the renormalized NNLO virtual term σ (2) can be expressed in the following general way: where we have used Catani's formula for the infrared singular behaviour of the two-loop QCD amplitudes [36][37][38], and we have defined the quantities: All the dependence on the Higgs trilinear coupling λ is embodied in the coefficient C LO . The explicit expression for the one-loop I (1) g and two-loop I (2) g insertion operators can be found in Ref. [36]. We recall here that they are functions of the dimensional regularization parameter ǫ, with poles up to 1/ǫ 2 and 1/ǫ 4 , respectively. The function f (ǫ) originates in the n-dimensional two-particle phase space, and verifies that f (0) = 1.
While the singular behaviour of the two-loop amplitudes can be anticipated, the finite contributions σ fin can only be obtained after performing the full two-loop calculation. The pole structure of our result agrees with the expressions in Eqs. (7) and (8) and the infrared-finite contributions for Higgs pair production can be cast into the form For simplicity, we set µ 2 R = s in the following expressions. We find that the one-loop contributions are given by The expansion of σ (1) fin is needed up to order ǫ 2 because of the double poles present in I (1) g . The F (1) contribution arises from the interference between FF(1) and LO, while R (1) originates from the interference of 2V(1) with the LO. The expansion up to O(ǫ 0 ) agrees with the result presented in Ref. [18].
The two-loop infrared regulated contributions take the following form: HH ) .
Here F (2) originates from the interference between the two-loop form factor-like diagrams FF (2) and LO contribution plus the square of FF(1), while V (2) arises from the square of the tree-level diagram 2V(1). The terms R (2) and I (2) combine the contributions of two interferences: 2V (2) with LO, and 2V(1) with FF(1). The Mandelstam variables are given by the expressions: while the integration limits t ± correspond to cos θ = ±1 and Q is the double-Higgs invariant mass. The last term in R (2) , originated on form factor-like contributions, vanishes if the two-loop corrections to the effective vertex ggHH are the same as those of gHH.

NNLO Soft-Virtual approximation
Expressed as in Eq.(8), the (finite parts of the) two-loop corrections are ready to be implemented in the NNLO soft-virtual (SV) approximation universal formula derived in Ref. [25]. We do not attempt for a full phenomenological analysis of the process at this level, and mostly use the SV approximation as a way to evaluate the impact of the new two-loop results in the cross section. Therefore, in Figure 2 we show the NLO and NNLO-SV K-factors for proton-proton collisions at the LHC with c.m. energy √ s H = 14 TeV, in terms of the invariant mass of the Higgs pair. Here the NNLO-SV approximation is defined by adding the pure NNLO-SV contribution to the full NLO result. At each order, we use the corresponding MSTW2008 set of parton distributions and QCD coupling [39]. The bands are obtained by independently varying the scales µ R and µ F in the range 0.5 Q ≤ µ R , µ F ≤ 2 Q, with the constraint 0.5 ≤ µ R /µ F ≤ 2. The LO cross section that normalizes the K-factors is computed at µ R = µ F = Q. We recall that we always rely on the heavy top quark limit, and that we use the SV approximation as defined in Mellin space. As mentioned before, since the coefficient C HH is still unknown we assume C H for the numerical results.
As can be seen from the plot, we find a large K-factor, with K SV NNLO = 2.37 for the total cross section, resulting in an increase of 23% with respect to the previous order (K NLO = 1.92). This value remains approximately constant along the entire Higgs pair invariant mass distribution, with the exception of the region near the threshold where the cross section is anyway very small. Despite of the still sizable corrections, it is noticeable the improvement in the perturbative expansion in the strong coupling constant, which shows the first signs of convergence at NNLO. It is only at this order than there is a (yet not very significant) overlap between two consecutive scale dependent  bands. We can also observe that the scale dependence is substantially reduced: the NNLO band results in a about a ±8% variation around the central value, more than a factor of two smaller than the corresponding NLO band.
We want to recall that in the case of single-Higgs boson production the soft-virtual approximation (compared to the full NNLO result) is known to be accurate to a few percent level. We expect it to be even better for Higgs pair production due to the larger invariant mass of the final state, which leaves less energy for extra hard radiation. In fact, we computed the NLO soft-virtual cross section, finding K SV NLO = 1.95, which differs from the full NLO result by less than 2%. In contrast, the heavy top quark approximation is not expected to be as good as for single-Higgs production since the invariant mass of the Higgs pair is not small compared to the top quark mass. Still a number of improvements can be applied to the current approximation, like keeping the exact full mass dependent LO expressions wherever they appear in the higher order expansion [18]. Future work may be directed either towards a full NNLO calculation (in the heavy top limit), or to compute subleading terms in the heavy top quark mass expansion.