Virtual corrections to Higgs boson pair production in the large top quark mass limit

We calculate the three-loop matching coefficient $C_{HH}$, required for a consistent description of Higgs boson pair production in gluon fusion through next-to-next-to-leading order QCD in the heavy top quark approximation. We also compute the $gg \to HH$ amplitude in $m_t \to \infty$ approximation in the full theory and show its consistency with an earlier computation in heavy-top effective theory.


Introduction
After the discovery of a Higgs boson at the LHC, detailed investigation of its properties becomes one of primary goals of ATLAS and CMS.Important among such studies is the exploration of the Higgs boson self coupling λ.In the Standard Model, this coupling is directly related to the Higgs field potential responsible for the symmetry breaking; in the broken phase, it induces couplings of three Higgs bosons between themselves.
Experimentally, information about λ is obtained from the process of Higgs boson pair production [1,2] which will be accessible after the high-luminosity upgrade of the LHC.It is well understood by now that observation of Higgs boson pair production is difficult and requires both, new ideas on how to isolate the HH signal from the background, and accurate predictions for the Higgs pair production in the Standard Model.In the past year we have witnessed significant advances in both of these directions.Indeed, building upon the early ideas of Refs.[3,4] it was suggested to study Higgs pair production in W + W − b b, γγb b, b bb b, and b bτ + τ − channels using substructure techniques [5][6][7], as well as utilize ratios of cross sections [8] for single and double Higgs production to reduce the theory uncertainty and obtain best sensitivity to Higgs boson self-couplings.It remains to be seen how these theoretical ideas will bare in real experimental searches, but the current consensus seems to be that the Higgs self-coupling can be measured with the accuracy between twenty and forty percent (see, e.g., Ref. [9]).
To interpret results of experimental measurements with this accuracy, one needs to ensure that Standard Model predictions for Higgs boson pair production are known with sufficient precision.Below we summarize the current status of theoretical computations of Higgs boson pair production in the Standard Model.The leading order predictions for gg → HH are known since long ago; they were computed in Refs.[1,2] where the exact dependence on all kinematic variables -primarily the top quark mass -has been taken into account.Improving on these results would have required the two-loop computations with massive internal (top quarks) and external (Higgs bosons) particles; currently, such computations are technically not feasible.Instead, a possible way forward is provided by studying the QCD corrections in the approximation where the top quark mass is taken to be much larger than all other kinematic invariants in the problem.Working to leading order in 1/m t expansion, one can integrate out the top quark and obtain an effective theory where Higgs bosons couple directly to gluons.Within such theory, next-to-leading order (NLO) computations for pp → HH become feasible and have been performed in Refs.[10] in m t → ∞ approximation while finite 1/m t corrections were calculated in [11].Recently, the next-to-next-to-leading order (NNLO) QCD corrections to pp → HH were computed in [12,13] in m t → ∞ approximation using the close analogy between pp → H and pp → HH production in effective theory.Soft-gluon resummations and the determination of dominant π 2 terms have been considered in [14] at next-to-next-to-leading logarithmic order.
In spite of tremendous progress with fixed order computations for double Higgs produc-tion, we note that NNLO QCD result of Refs.[12,13] is formally not complete.Indeed, at the NNLO QCD accuracy for Higgs pair production, one needs the Wilson coefficient C HH which was not available when Refs. [12,13] were written.The goal of this paper is to perform the computation of the C HH Wilson coefficient and therefore provide the last missing ingredient required to describe the Higgs boson pair production through NNLO QCD in the large-m t approximation.
Before we proceed with the computation of the Wilson coefficient, a word of caution about the validity of large-m t approximation is in order.Indeed, it is well-known that for Higgs pair production the m t → ∞ limit provides a poor description of both the total cross section and kinematic distributions.In such a situation it is far from clear that extending m t → ∞ computations to NNLO, as was, e.g., done in Refs.[12,13], is a sensible way to estimate higher order corrections to Higgs boson pair production.Understanding the validity of this approach was the primary goal of Ref. [11] 1 where it was shown that, for a properly chosen leading order cross section, the 1/m t effects at NLO are moderate, in the 15 − 20 percent range.If we assume that the same remains true at NNLO, we conclude that m t → ∞ NNLO QCD corrections can be used to provide a reliable estimate of NNLO QCD corrections with the full top quark mass dependence.
The remainder of the paper is organized as follows.In the next Section we introduce the effective Lagrangian for single and double Higgs production in gluon fusion.In Section 3 we describe the matching calculation of C HH .In Section 4 we discuss the computation of the virtual corrections to the gg → HH cross section in the full theory which serves as the cross-check of some results presented in Ref. [12].In Section 5 we present our conclusions.

Effective Lagrangian for Higgs pair production
The leading order effective Lagrangian that describes interactions of any number of Higgs bosons with gluons in m t → ∞ limit is given by In Eq. ( 1) H and v are the Higgs boson field and the vacuum expectation value, respectively, O 1 = 1/4G a µν G µν,a , where G a µν is the gluon field strength tensor, and α s is the strong coupling constant.This Lagrangian is modified in higher orders of perturbative QCD.To account for this, we restrict Eq. (1) to describe interactions of gluons with up to two Higgs bosons,2 and write The matching coefficients C H and C HH incorporate radiative effects of top quarks that are integrated out from the Standard Model; they are given by perturbative series in the strong coupling constant.
Superscripts "0" in Eq. ( 2) indicate that operator renormalization has not yet been performed, so that both C 0 H and C 0 HH as well as matrix elements involving O 0 1 are ultraviolet divergent.Following Ref. [17] we can write C 0 This procedure leads to finite coefficient functions C H and C HH .In Eq. ( 3) we used s (µ) to denote the MS strong coupling constant in a theory with five active flavors; we will use this notation throughout the paper.We also used standard notation ) and T F = 1/2 are SU(N c ) color factors and n l = 5 is the number of massless quarks.
It is convenient to introduce the perturbative expansion of C H and C HH via (0) HH = 1.Note that equality of C H and C HH at leading order follows from the Lagrangian in Eq. (1).We have chosen to parametrize C H and C HH in terms of the five-flavor strong coupling constant.

Direct calculation of matching coefficients
Since C H and C HH are matching coefficients between full and effective theories, it is convenient to derive them as follows: compute amplitudes of any physical process that depends on one or both of them in full and effective theories and adjust C H and C HH in such a way that the two amplitudes agree.Of course, to determine C H and C HH independently, we need to consider two, rather than one, physical processes; we choose them to be i) Higgs boson production in gluon fusion gg → H and ii) Higgs boson pair production in gluon fusion gg → HH.The amplitude of the first process depends on C H .The amplitude of the second process depends on both C H and C HH .
We begin with the computation of C H and consider the process g(q 1 )g(q 2 ) → H with q 2 1 = q 2 2 = 0 and q 1 • q 2 = m 2 H /2. We are interested in the behavior of this process in the limit q 1 ∼ q 2 ∼ m H ≪ m t where the scattering amplitude can be computed in both full and effective theories.The requirement that the two amplitudes are equal up to power-suppressed terms reads lim ( We now study this equation order-by-order in QCD perturbation theory.At leading order, the amplitude in the full theory is given by the one-loop triangle gg → H diagram which can be Taylor expanded in external gluon momenta.The amplitude in the effective theory follows from the Lagrangian Eq. ( 1) and reads Upon equating full and effective theory amplitudes, we find C H = −α s /(3π), which is the first term in the expansion of the result in Eq. ( 4).
At NLO, the situation changes for the following reasons.On one hand, loop corrections to gg → H amplitudes in the effective theory appear.On the other hand, Taylor expansion of gg → H amplitude in small momenta and the Higgs mass no longer gives correct full theory amplitude even in the limit q 1 ∼ q 2 ∼ m H ≪ m t since non-analytic dependencies on s and m 2 H do, in general, appear.To cure these problems, the large-mass expansion procedure [19] is applied to Feynman diagrams that contribute to the full theory amplitude.The large-mass expansion splits all loop momenta into soft k ∼ q 1 ∼ q 2 ∼ m H and hard k ∼ m t and allows systematic Taylor expansions of integrands in both of these regimes.Scaling of loop momenta determines scaling of integrals since d d k| soft ∼ s d/2 and d d k| hard ∼ m d t .Since the gg → H amplitude necessarily involves at least one loop of top quarks, only one of the two loop momenta can be soft.For the NLO amplitude in full theory this implies3 We note that hard part of the amplitude A full is obtained by Taylor expansion of integrands of loop integrals in powers of q 1,2 /m t and m H /m t ; therefore, to obtain A hard NLO only two-loop vacuum integrals need to be computed.On the contrary, the soft part of the amplitude requires computation of integrals of form-factor type which depend on external soft kinematic parameters.When quantum corrections are computed in the effective theory, only soft contributions are generated.Therefore Since we are interested in C H which, by construction, can not depend on s, Eqs. ( 5), ( 7) and ( 8) can be matched provided that In Eq. ( 9), ζ 0 3 is the decoupling constant of the gluon field (cf.Refs.[17,20]), which is needed for the (on-shell) wave function renormalization of external gluons induced by the top quark loops.
The result shown in Eq. ( 9) allows us to obtain the matching coefficient C H by ignoring all loop corrections to gg → H amplitude in the effective theory and by computing Taylor expansion of relevant diagrams in q 1,2 /m t and m H /m t in the full theory.Extension of the above discussion to NNLO is straightforward.We write and solve for C H order by order in the strong coupling constant α s .
Before we show the (known) result for C H , we would like to make a few technical remarks.First, we note that it may be inconvenient to deal with external gluon polarization vectors (cf.Eq. ( 6)) in multi-loop computations.If so, one can use an appropriate projection operator to avoid them.A convenient choice, that respects transversality of the gluon polarization vectors, is which transforms the leading order amplitude in Eq. ( 6) into Second, we note that we first renormalize the top quark mass on-shell, and the strong coupling α s in the MS scheme with six active flavors.We then apply the two-loop decoupling relations to transform α s to α s .We note that in this relation the O(ǫ) terms have to be kept at one-loop order since the two-loop term of C 0 H has an 1/ǫ pole whereas the one-loop term is finite.The finite result for C H , obtained via C 0 H /Z O 1 is given by [17,21,22] where α s = α  We are now in position to extend the above discussion in such a way that the computation of C HH becomes possible.To this end, we choose the gluon fusion process where two Higgs bosons are produced, g(q 1 )g(q 2 ) → H(q 3 )H(q 4 ).We then apply the same reasoning as for the single Higgs boson production and compare amplitudes for A gg→HH computed in full and effective theories assuming that q 1 ∼ q 2 ∼ q 3 ∼ q 4 ∼ m H ≪ m t .We note, however, that there is a subtlety in this case that is related to the fact that pairs of Higgs bosons can not only be produced through the ggHH operator but also through one or two ggH operators in effective theory, see Fig. 1.This can occur in two different ways.For example, already at leading order, double Higgs production in the full theory receives contributions from a box diagram gg → HH and from a triangle diagram gg → H * where the virtual Higgs boson splits into a HH pair.The second contribution has nothing to do with the matching coefficient C HH .Our master formula that is based on equating amplitudes in full and effective theories automatically takes care of this since an identical contribution is also generated in the effective theory through a local interaction vertex ggH.Hence, diagrams with intermediate off-shell Higgs bosons cancel exactly between full and effective theory amplitudes so that at leading order only gg → HH box diagram in the full theory is needed to obtain the Wilson coefficient C HH .Similar subtleties occur in higher orders, see, e.g., Fig. 1(c).Nevertheless, separation of loop momenta into soft and hard and the understanding that effective theory loops are always soft allows us to consider only hard contributions in the full theory and equate them directly to products of matching coefficients and various tree amplitudes in the effective theory.We therefore obtain the following generalization of Eq. ( 10) valid in the case of Higgs pair production When writing Eq. ( 14) we introduced labels 1PI and 1PR, to denote one-particle reducible and one-particle irreducible contributions in both full and effective theory.Moreover, we separated various one-particle reducible contributions on both sides of Eq. ( 14) into those that involve and do not involve the triple Higgs boson coupling λ.We also note that these one-particle reducible contributions contain poles in soft kinematic parameters, so that it is more appropriate to talk about Laurent rather than Taylor expansion of full theory amplitudes in Eq. ( 14).However, all kinematic poles cancel exactly between the left-and the right-hand side of Eq. ( 14), as required by the consistency of effective theory.
We note that Eq. ( 14) can be immediately used for the computation of the matching coefficient C HH since this is the only unknown quantity there.However, before doing that, it is important to realize that Eq. ( 14) can be significantly simplified.Indeed, as the immediate generalization of the leading order discussion in the previous paragraph, we observe the exact matching between one-particle reducible contributions to Eq. ( 14) caused by nonvanishing triple Higgs boson coupling; this allows us to remove A eff tree,1PR,λ =0 and A full 1PR,λ =0 from both sides of Eq. ( 14).It is natural to think that further simplifications are possible.For example, it is easy to imagine that Z 2 O 1 C 2 H A eff tree,1PR,λ=0 and A full 1PR,λ=0 should match exactly on the two sides of the equation and can be removed.Indeed, this is what happens through two loops but the two contributions do not match exactly at three loops leaving a remainder that gets re-absorbed into C HH matching coefficient.Finally, we want to point out that all calculations have been performed for arbitrary gauge parameter ξ which drops out in the final result, a strong check of the correctness of our calculation.
The final result for C HH that we obtain can be summarized as follows.Using the parametrization of C H and C HH in Eq. ( 4), we find

HH , ∆
(2) where n l is the number of massless quarks.We note that the difference between C (2)

HH and C
(2) H is significant.Indeed, for n l = 5 and µ = m t , we find H ≈ 6.15, (16) which implies that C (2) H ≈ 1.8.We note that in the computation of Refs.[12,13] it was assumed that 0 < C (2) H ; Eq. ( 16) shows that our result for C (2) HH is within this interval but close to its upper boundary.The numerical effects on C HH = C H on the cross section is investigated in Section 5.

4
Virtual corrections to gg → HH production at NNLO In the previous Section we computed the matching coefficient C HH by comparing hard contributions in the full theory and tree contributions in the effective theory.In this way, we only had to compute vacuum bubble integrals to obtain C HH .However, we can calculate the full gg → HH amplitude in m t → ∞ approximation if we account also for soft contributions in the full theory.Then we obtain the NNLO virtual corrections to gg → HH amplitude independent of effective theory computations.
How difficult is it to compute soft contributions through NNLO for the double Higgs production?It turns out that it is not so hard.Indeed, since we have to deal with at most three-loop diagrams in the full theory and since at least one of those three loops has to be hard, the most complicated soft integrals that need to be computed are twoloop three-point functions and one-loop four-point functions with all internal and two external lines massless.All such integrals are known which means that we can obtain full gg → HH amplitude from the full theory.
We consider production of the Higgs boson pair in gluon collisions g(q 1 )+g(q 2 ) → H(q 3 )+ H(q 4 ) and introduce Mandelstam variables s = (q 1 + q 2 ) 2 = (q 3 + q 4 ) 2 , t = (q 1 − q 3 ) 2 = (q 2 − q 4 ) 2 and u = (q 1 − q 4 ) 2 = (q 2 − q 3 ) 2 .Gluons and Higgs bosons are on the mass shell, q 2 1,2 = 0 and q 2 3,4 = m 2 H .We write virtual contributions to gg → HH differential cross section as where again α s = α s (µ).The leading order cross section in Eq. ( 17) can be written as where We note that C LO is the sum of two leading order contributions to Higgs boson pair production cross section associated with box and triangle diagrams and that ǫ-dependent factors in Σ LO originate from the d-dimensional two-particle phase space, and the average over gluon polarizations.The higher order ǫ terms in Eq. ( 18) differ from such terms in Ref. [12] since the matching coefficients used in [12] are strictly four-dimensional.We can emulate this effect in our calculation and reproduce the results of Ref. [12].Similar comments also apply to the NLO results given below.
Representative Feynman diagrams contributing to the amplitude A gg→HH can be found in Fig. 2. We compute the differential cross sections using large-mass expansion [19] with the help of the C++ program exp [23] that factorizes all integrals into hard ( soft (two-loop three-point and one-loop four-point) integrals.As we already noticed, all such integrals can be computed in a straightforward way.Once this is done, we obtain perturbative results for the virtual corrections to the gg → HH cross section through NNLO in the heavy top approximation.
We note that, since virtual corrections are computed in the full theory, the results are made ultraviolet finite by means of standard renormalization procedure.In particular, no matching computations are required.Therefore, by comparing the result of the full theory computation with Ref. [12], one can independently verify the effective theory computations reported there and, at the same time, check the consistency of C HH computation described in the previous Section.
To present the results for virtual corrections, we follow the standard practice and isolate infrared-divergent pieces using Catani's representation of scattering amplitudes [24].For ultraviolet finite gg → HH scattering amplitude, we write The two operators I (1,2) g depend on QCD color factors C A , C F and n l T F , the Mandelstam variable s and the dimensional regularization parameter ǫ.In the limit ǫ → 0, Figure 3: One-(a) and two-loop (b) form-factor contributions which lead to F (1) and F (2) .Multiplying (c) and (d) with the LO amplitude leads to R (1) and R (2) .V (2) is obtained from squaring contribution (c).
develop 1/ǫ 2 and 1/ǫ 4 singularities, respectively.On the other hand, A (1,2),fin contributions to NLO and NNLO amplitudes are finite.The exact form of I (1,2) g operators can be found in Refs.[24,25]; we do not reproduce them here.Using the representation of scattering amplitude Eq. ( 20), we write the virtual contributions to gg → HH cross sections as dσ It follows from Eq. ( 21) that all divergent contributions are proportional to either leading or NLO cross sections.Since the leading order cross section has already been given in Eq. (18), it is sufficient to provide results for finite NLO and NNLO contributions.
In the following we present our results in a way which allows for a simple comparison with Ref. [12].Contributions to gg → HH amplitude split naturally into two classes -one that corresponds to only one effective vertex (ggH or ggHH; they occur after shrinking the vacuum bubbles to a point) and the other one that involves two ggH vertices.
In the former case, all soft contributions are reducible to three-point functions and are proportional to the leading order amplitude C LO .Diagrams with two effective vertices start to contribute at NLO and the corresponding one-loop corrections are needed at NNLO.For convenience we show sample diagrams up to NNLO in Fig. 3 where also the notation for the individual contributions is introduced.Following this classification, we write the finite contribution to the one-loop cross section as dσ where the first term in square brackets is the contribution of diagrams with a single effective vertex and the second term is the contribution of all diagrams with two effective vertices.We perform a similar decomposition at NNLO and write dσ where the new element V (2) is the contribution of NLO diagrams with two effective vertices [cf.Fig. 3(c)] squared.
In addition to soft contribution described so far, hard contributions also enter equations ( 22) and (23).They can be computed directly using full theory diagrams without resorting to separating these hard contributions into C H and C HH .We can then combine C H and C HH results described in the previous Section with the effective theory computation reported in Ref. [12] and compare the result with the full m t → ∞ computation described in this Section.The two results agree which provides a good consistency check for both, the effective theory computation and the calculation of the C HH Wilson coefficient reported in the previous Section.

Conclusions
We computed the three-loop Wilson coefficient of a G 2 H 2 operator that describes interactions of two Higgs bosons with gluons in the approximation that the top quark mass is infinitely large.This is the last missing ingredient that is required to perform consistent NNLO QCD computation of Higgs pair production in the large-m t limit.Our main result -the three-loop contribution to the Wilson coefficient C HH -is given in Eq. (15).We have also computed virtual corrections to Higgs pair production in gluon fusion in the full theory using asymptotic expansions in the inverse top quark mass and verified consistency of our C HH computation with the calculation of gg → HH virtual corrections within the effective field theory [12].
An interesting feature of the computed three-loop corrections is that they break the equality C H = C HH that persists through two-loops.Therefore, their main effect is to change the relative contributions of the box and triangle diagrams to double Higgs production.Since box and triangle contributions cancel exactly at the threshold for producing the two Higgs bosons, the relatively small difference between C H and C HH gets kinematically amplified. 5Indeed, using the relation between Higgs boson self-coupling, the vacuum expectation value and the Higgs boson mass 2λ where ∆ (2) HH from Eq. ( 15) is used.The strong kinematic enhancement at the threshold s = 4m 2 H is evident.Numerically, assuming m H = 125 GeV and α s = 0.11, the correction to cross section for gg → HH computed using C H = C HH approximation amounts to 6.4 percent at √ s = 270 GeV and 1.7 percent at √ s = 400 GeV.The change in the total hadronic cross section pp → HH amounts to 1%, compared to the case C (2) HH .While all these corrections are quite moderate, the change in threshold behavior is interesting and is qualitatively different from a relatively uniform enhancement of lower-order cross sections provided by soft QCD effects.

( 5 )
s (µ) is the MS coupling constant defined in the theory with n l = 5 massless flavors and m t is the pole mass of the top quark.

Figure 1 :
Figure 1: Effective-theory diagrams with ggH and ggHH operators contributing to the double Higgs boson production.

Figure 2 :
Figure 2: Sample Feynman diagrams contributing to the amplitude A gg→HH .
2 v = m 2 H , we write the relative correction as dσ C H =C HH − dσ C H =C HH dσ C H =C HH