ELSEVIER

CrossMaik

Available online at www.sciencedirect.com

ScienceDirect

Nuclear Physics B 900 (2015) 412-430

www. elsevier. com/locate/nuclphysb

Higgs boson pair production: Top quark mass effects at

NLO and NNLO

Jonathan Grigoa, Jens Hoffb, Matthias Steinhausera*

a Institut für Theoretische Teilchenphysik, Karlsruhe Institute of Technology (KIT), 76128 Karlsruhe, Germany b Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany

Received 6 August 2015; accepted 21 September 2015

Available online 28 September 2015

Editor: Tommy Ohlsson

Abstract

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 Mt result. Since the 1 /Mt 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.

© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3.

1. 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 0(103)

* Corresponding author.

E-mail address: matthias.steinhauser@kit.edu (M. Steinhauser).

http://dx.doi.Org/10.1016/j.nuclphysb.2015.09.012

0550-3213/© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3.

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-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-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 corrections the exact dependence on the top quark mass and the center-of-mass energy is known [9,10]. At next-to-leading order (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 0(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 threshold-enhanced 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.

2. 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

(0) as (1) / as \2 (2)

Oij ^hh+x(s, p) = SigSjgagg)(s, p) + naij(s, p) + n) a,'J)(S, P) + "' , (1)

where as = a(5)(^), .Js is the partonic center-of-mass energy and ij e {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

p = Th, (2)

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 \x = \xr = \xf. For later convenience we introduce for agg^HH+X

ic r,\ „LO , x NLO , O NNLO , m

Ogg^HR+x(s, p) = a + So + So + ..., (3)

and denote the sum of the first two terms on the right-hand side by aNLO = aLO + SoNLO.

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/2 of the NLO contribution to a(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

Y^ ^NLO n a(1) LCgg,np

o(l) _o(0) a(N) a (N) — "gg,exp _ n=0 (4)

"gg,N "gg^act^gg , Agg ^ (0) N , W

gg,exp „LO pn

where og°jexact contains the exact dependence on s and p.

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

s d (0) dogg,exp

,r(1) I H n2 dggg,exact dQ2

o(g) = / dQ --"7m , (5)

gg ^ dQ2 dCT(0)

J UQ dogg,exp

4m^ dQ2

where 'exact' and 'exp' remind whether exact or (in p) expanded results are used. Q2 is the invariant mass squared of the Higgs boson pair. Expressed in terms of do/dQ2 it is possible to re-write Eq. (4) as

o(°) S do(i)

(1) "gg, exact / ^2 dogg,exp

°gg = (0) dQ ■ (6)

ogg, exp J dQ 4m2„

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 Q2 > 4M2.

For the virtual corrections, which are proportional to S(s — Q2 ), one has immediate access to the Q2 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) cannot be used. For the construction of the soft-virtual (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 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 a. Note, however, that the arguments also hold for da/dQ2. In a first step we split a according to

a _a virt+ren + a real+split

where the two terms on the right-hand side correspond to the virtual corrections (including ultraviolet 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 a can be written as

a _ £sv + ^h '

where "SV" refers to "soft-virtual" and "H" to "hard". This is achieved by splitting avirt+ren into a divergent and a finite term and separating areal+split into a (divergent) soft and (finite) hard contribution. The soft contribution is combined with avirt+ren to obtain £SV such that

real+split

^SV _ ^div + ^fin + ^soft ,

real+split LH _ ^hard •

Esv and Eh are separately finite. Edv is constructed following Ref. [23]; explicit expressions can be found in Refs. [22,16]. Note that the finite part is constructed solving the equation

virt+ren

_ £fin + ^div .

for Efin. In this paper avirt+ren is computed including top quark mass effects. Mass effects are automatically taken into account in Ediv and £soft+split 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 da/dQ2 and thus a factorization as suggested in Eq. (5) can be performed for the soft-virtual contribution.

To be more specific we present the NLO and NNLO differential version of Eq. (10) which reads [22,16]

(1) v, fin

(2) v,fin

I (D I

(1) v, fin

I (1) I

I (1) 2

f(2) I

1 In fact, for the virtual corrections Eqs. (4) and (5) are equivalent.

where t = (q1 - q3)2 with q1 (q3) being the incoming (outgoing) momentum of a gluon (Higgs boson), a(0) = a LO

and explicit expressions for the operators I(1'2) can be found in Ref. [23]. We adopt the notation from Ref. [22] and parametrize radiative corrections to the partonic cross section via

= aLOzG(z).

G(z) = 5(1 - z) + a^G(1)(z) + (a) G(2)(z) + ..

= Gsv(z) + Gh(z),

where the renormalization scale dependence in the strong coupling constant as is suppressed. From Eq. (12) one obtains for the total cross section

--f <2 = /

4rnl 1-5

dzaLO(zs)G(z),

5 = 1 -■

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 GH(z). Actually, at NNLO we only have Gsv(z) at hand and at NLO only the heavy top expansion of /!_5 dzaLO(zs) GH(z) is available to us.

Explicit results for Gsv(z) can be found in Ref. [22] including higher order terms in e specifying, however, the renormalization and factorization scale to ^fs ~ yQ2. For completeness we provide the NLO and NNLO results for Gsv(z) for generic / in the limit e ^ 0. The results read

G<iv=D-

-2n2 a,

(1) fin LO

- 4caldo + 8CD

4n2 2 2 11n2 2

-C2 L2 +- -2

3 A + 108

t(2) fin

(1) fin ,LO

( 1 2n2

+CM -1 - —

(1) fin

+ CA - +

1 11n2

- 38f3

+ CA — +

'607 517n2 n4

- - 407^3 \

' 2n 2 a™ — fi + ni\-81 -

37 + 18

+ D0{ L [ - — CA + 3CAnn + CAni 27 9

56 4n2

134 10n2

CA( + 3

+ Ca —n - 4-

fin TLO

404 22n2

+ CA> + 9

+ 78ft + Dx

2 2 /'44 2 8 16CAL2 + L( — CA - 3CAni

268 20n2

+ Ca \ - — ni + 8

"fin TLO

44 2 2 8 - - CA - 48CA L + 3 CAni

+ 32Ca D3.

where CF = 4/3, CA = 3 are QCD color factors, nl = 5 is the number of massless quarks, L = log(^2/s) and

D-1 = S(1 - z)

The quantities o1(n) and a fin are evaluated for ^ = */s. They are obtained from daVoi^/dt and

do,( fin/di 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.

Dn>0 =

logn(1 - z)

3. 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 LO PDFs to compute N 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 constant as(MZ) that we use in our computation from the MSTW PDF fit:

aLO = 0.13939, asNLO = 0.12018, asNNLO = 0.11707. (19)

We renormalize the top quark in the on-shell scheme and use Mt = 173.21 GeV [25]. For the Higgs boson we use mH = 125.09 GeV [26]. For numerical results shown in this section the renormalization and factorization scales have been set to 2mH.

In Fig. 1 we show the results for the partonic cross section as a function of from Fig. 6 of Ref. [12], see solid lines. The splitting of these results into soft-virtual and hard contributions is shown as dotted and dashed curves, respectively. For the complete result and soft-virtual approximation the infinite-top quark result corresponds to the (red) second line from below. The lowest line includes in addition the 1 /Mt2 corrections. The curves including higher order mass corrections are above the infinite-top quark result. The topmost curves include 1 /M}2 terms. The hard contribution shows a quite flat behavior above ^Js > 400 GeV and exhibits smaller shifts when including top mass corrections. Furthermore, the infinite top mass result corresponds to the

sv+h -

sv-----

■I 0.4

Vï[GeV]

Fig. 1. Partonic cross section as a function of the center-of-mass energy including various orders in the inverse top quark mass. The dotted and dashed curves show the breakup of the complete result (solid lines) into soft-virtual and hard contribution. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 2. Partonic cross section as a function of the center-of-mass energy including various orders in the inverse top quark mass. The dashed and solid curves correspond to the factorization for the total and differential cross section, respectively. The color coding is taken over from Fig. 1.

topmost curve and the lowest curve includes 1 /M}2 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 color coding is as in Fig. 1. For lower values of the two approaches lead to comparable results. However, the DF curves have their maximum at lower values of s/s and lead to a smaller cross section for larger values of . Furthermore, one observes a drastic improvement in the convergence behavior when including higher order mass corrections. In particular, for = 400 GeV the difference between the infinite top mass result and the one including 1 /Mt12 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), GSV(z), is constructed for

sv+h-----

sv(df)+h -

300 400 500 600 700 800 v^ [GeV]

\fs [GeV]

Fig. 3. Partonic NLO K factor for the factorization performed at the level of the total (dashed) and differential (solid) cross section.

the limit z ^ 1 and thus it is possible to replace Gsv(z) by f(z)GSV(z) with f(1) = 1 (see, e.g., discussions in Refs. [27,28]). The hard contribution is modified accordingly such that the sum of Ssv + ^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 well2 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/Q2) which is justified since Vs ^ vQ2 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

aLO+5aNLO

K NLO = -aLO-• (20)

Results for the two methods to factorize the exact LO term are plotted in Fig. 3 as a function of Vs 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 p orders is smaller. Furthermore, it is interesting to note that for DF the top quark pair threshold behavior 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]).

Fig. 4 shows the hadronic cross section aH for Higgs boson pair production including NLO corrections as a function of ^/S^ut which is a technical upper cut on the partonic center-of-mass collision energy. It is introduced via

an(sn,scut) = J dr ^) (r)a(rsH)0(scut - tsh), (21)

4m2H/sn

2 Note that the corresponding curves are not shown in Fig. 2.

Fig. 4. NLO hadronic cross section and K factor as a function of „Jscut.

where the luminosity function is given by 1 1

(dd^f) T) = / dxi y dx2fg(xi)fg(x2)S(r - X1X2). (22)

fg(x) are the gluon distribution functions in the MS scheme. Note that in the soft limit ^/Su is a good approximation to Q2. The various lines in Fig. 4 correspond to the inclusion of different orders in p at NLO. For convenience we show on the right end of Fig. 4 the total cross section for ^JSH = 14 TeV. Note that the approximation used for the computation of the pn terms is not valid for large values of VScut (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 obtained from Eq. (20) by replacing a by aH and using NLO PDFs in the numerator and LO PDFs in the denominator. KNLO raises close to threshold, however, for ^/Su > 500 GeV one observes a flat behavior of KNLO « 1.6 (for x = 2mH).

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 split3 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 countert-erms; dotted lines) and the real-radiation parts which include the contributions from the splitting functions (dashed lines).

We observe that for < 400 GeV, the region where our approximation is valid, the top quark mass corrections to the real radiation part (upper dashed curves) reduce the infinite top result by about 10%, in agreement with the observations of Ref. [13]. On the other hand, the top mass

3 Note the individual terms contain 1/e poles which cancel in the sum. In Fig. 5 the finite contributions are plotted.

260 280 300 320 340 360 380 400 y/s [GeV]

Fig. 5. Splitting of partonic cross section (solid lines) into real (upper dashed lines) and virtual (lower dotted lines) contributions (see text for details).

effects to the virtual contribution leads to a positive shift as compared to the effective-theory result. Summing real and virtual corrections leads to an overall positive effect from top mass corrections as can be seen by the solid curves, see also Fig. 2. Note that up to *fs & 400 GeV top mass corrections are dominated by the 1 /Mf terms.

4. 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 as described in Section 2 and use Eq. (17) to evaluate the partonic and hadronic cross sections including 1 /Mf correction terms.

4.1. 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.

4.1.1. Amplitude gg ^ HH

NNLO calculations require corrections up to three-loop order to the process gg ^ HH. Typical contributing Feynman diagrams at LO, NLO and NNLO are shown in 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.

vssafioiv

3>.....

jjsBsar

IOOQOQOOJ

iQOOQOQOJ

innq^on"

i8Sj&8!!S>

Fig. 6. One-, two- and three-loop Feynman diagrams contributing to the process gg ^ HH. Solid lines refer to top quarks, curly lines to gluons and dashed lines to Higgs bosons.

Ï2 Ï3

Fig. 7. One- and two-loop master integrals needed after applying asymptotic expansion to the amplitude gg ^ HH. All

2 2 2 4 2

internal lines are massless, q1 = q^ = 0, and qj = q^ = mH.

The latter are shown in Fig. 7 where we have q2 = q2 = 0, and qj = qj = m2H. Analytic results for all master integrals can be found in the literature (see, e.g., Refs. [36-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

= -L( V-

Go(-1; x)Go(0; y) - Go(0; y)G0(-1/y; x) + 2G0(-1, 0; x)

- 2Go(-1/y, 0; x) + e

- 2inGo(-1, 0; x) - Go(-1/y; x)Go(0,0; y)

+ Go(-1; x)( - inGo(0; y) + Go(0, 0; y)^ + 2inGo(-1/y, 0; x)

+ Go(0; y)(i*Go(-1/y; x) + Go(-1, -1; x) + 2Go(-1, 0; x) - Go(-1,-1/y; x) + Go(-1/y,-1; x)

- 2G0(-l/y, 0; x) - Go(-1/y, -1/y; xj) + 2G0(-1, -1, 0; x) + 4G0(-1, 0, 0; x)

- 2Go(-1, -1/y, 0; x) + 2Go(-1/y, -1,0; x)

- 4Go(-1/y, 0,0; x) - 2Go(-1/y,-1/y, 0; x)

+ O(e2)\ , (23)

1 1 + ^S 4m2

y=* • x = ttvt S=1 - (24)

Go({wi}; z) are Goncharov Polylogarithms [39] with weight {wi} and argument z defined through

Go(wi,w2,...,w„; x) = dt-Go(w2,...,w„; t), (25)

J t - wx

with wi, x e C and

Go(0n; x) = -1 log" x . (26)

The functions of the term in Eq. (23) can be expressed in terms of logarithms and diloga-rithms via

Go(0; y) = log(y), Go(-1; x) = log(1 + x) , Go(-1/y; x) = log(1 + xy), Go(-1, 0; x) = log(x) log(1 + x) + Li2(-x) , Go(-1/y, 0; x) = log(x) log(1 + xy) + Li2(-xy). (27)

We have cross checked the numerical result for the e0 and e1 terms of I1(4) in Eq. (23) against FIESTA [40].

Using this method we have computed terms up to order 1 /M}2 at NLO [12,20,21] and terms up to order 1 /M4 at NNLO. As an important check we have computed the 1 /M2 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 da/dt. Integration over phase space then leads to da/dQ2. For the results in Section 5 this integration is performed numerically.

4.1.2. 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

JJÛfifi, vQflQSLj

IgpftO/ vQflfi$L>

■nnnn ' \ onnoi

Lfioflo» vûflûsij—*—LÛÛJ îaa,

iQgfiJb vûQûflJ

, ^TSSSSÎT^

vDfifiSH

OOOWtt

BSSSS3.

fiSSISO

vQfiSSLi

Fig. 8. LO, NLO and NNLO Feynman diagrams needed for the forward scattering amplitude gg ^ gg. Solid lines refer to top quarks, curly lines to gluons and dashed lines to Higgs bosons. At NNLO only virtual contributions are shown. Wavy lines denote the cuts.

00000000't (b)

Fig. 9. Resulting Feynman diagrams after shrinking the top quark loops to a point according to the rule of asymptotic expansion. The blob represents one-loop vacuum integrals.

and filter [21,41] the amplitudes describing the virtual corrections to gg ^ HH. Typical Feynman diagrams are shown in Fig. 8.

• FORM [42,30] code is then generated by passing the output via q2e [31,32], which transforms Feynman diagrams into Feynman amplitudes, to exp [31,32].

• Our in-house FORM code applies projectors (—gxv) 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:

1. four one-loop vacuum integrals times one-loop integrals which are already present at LO, see Fig. 9(a);

2. three one-loop vacuum integrals times two-loop integrals, see Fig. 9(b);

3. two one-loop vacuum integrals times three-loop integrals, see Fig. 9(c);

4. two-loop times one-loop vacuum integrals times two-loop integrals;

5. three-loop times one-loop vacuum integrals times one-loop integrals.

Fig. 10. Phase space master integrals occurring in the amplitude for gg ^ gg. The first line contains LO and NLO integrals. The integrals in line two and three are needed for the NNLO virtual corrections. Single and double lines represent massless and Higgs propagators, respectively. Double lines with gray-shaded interspace (last two diagrams in first row) correspond to Higgs boson propagators which shall not be cut. A cross marks an inverse propagator.

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 S, we use the in-house 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 Ref. [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 S = 1-4m2H/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 T 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, m2H, m2H are given in Eq. (23). After expanding in S the remaining phase-space integration can be performed analytically. We have computed expansion terms up to order S10 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 S.

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 |^NLO|2 or (^LO^*NNLO + ^NNLO^*LO). In the former case the two-loop integrals originate from the product of two one-loop contributions containing factors (-1 + i0)e and (-1 - i0)e, respectively, which finally leads to (-1 + i0)e (— 1 - i0)e = 1. In the other case

one has (-1 + i0) + (-1 - i0) = 1 - 4n2e2 + O(e3). 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 f when computing the 1 /Mt2 terms. Furthermore, also the 1 /Mf corrections could be evaluated (for f = 0) whereas in the optical theorem approach only 1 /M2 terms could be computed in Feynman gauge. However,

I °-4

250 350 450 550 650 750

yfs [GeV]

Fig. 11. Comparison of the LO, NLO and NNLO contributions to the partonic cross section. At LO the exact result is shown and at NLO and NNLO the first three terms in the large-Mt expansion are shown. For all curves the NNLO-value for as is used. For the renormalization and factorization scale we use | = 2mH.

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 /Mf corrections. Note that we have da/dg2|virt ~ S(Q2 - s). Using Eq. (10) we construct afin which enters GSV in Eq. (17). The differential and total cross section is then obtained with the help of Eqs. (12) and (15).

5. Improving NNLO

In this section we discuss the effect of the 1 /Ml and 1 /Mf 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 afin which enters Eq. (17). The result for the partonic cross section is shown in Fig. 11 as a function of the partonic center-of-mass 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 *Js 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 p decreases the infinite top quark mass result which is overcompensated by the p2 term resulting in a small positive correction.

In Fig. 12 we compare the NNLO-SV contribution for two different scales, x = 2mH and i = \fQ2. 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.

4 In this plot we only include mass corrections up to order p2 at NLO to have a direct comparison with NNLO.

/i = 2mn

f\ n =

A NNLO-SV-----

J• \ z • NNLO-SV "

300 400 500 600 700 800 y/i [GeV]

Fig. 12. NNLO partonic cross section for different scales and for f(z) = 1 (dotted) and f(z) = z (solid). Only the p0 result is shown.

Fig. 13. Hadronic LO, NLO and NNLO-SV (with f(z) = z) cross sections as a function of -,/scut. For their evaluation the respective value of as is used. At LO the exact result and at NLO only the p0 term is shown. At NNLO the p0, p1 and p2 results are plotted. The results in the right panel with "to" at the bottom correspond to the prediction of the total cross section. For this plot \i. = 2m^ has been used.

The hadronic cross section as a function of ^/Su is shown in Fig. 13 for ^ = 2mH. At LO the exact result is used and the NLO curve is based on the p0 results. Using instead the p6 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 p0, p1 and p2. As at NLO one observes good convergence up to ^/scut ^ 400 GeV. For higher values of ySu the p1 curve is below and the p2 curve above the infinite top quark mass result leading to a ±2.5% effect for the total cross section on the rightmost part of the plot. To be conservative, we thus estimate that the NNLO top quark mass effects lead to an uncertainty of ±5%. Note that the NNLO corrections amount to about 20% of the LO result.

Fig. 14 shows the hadronic K factor at NNLO which is defined by

lo , x^nlo , £ NNLO\ | NNLO \aH + b°H + baH jlNNLOpdfs

KH =--LO-' (28)

&H lLOpdfs

NLO|„o z • NNLO-SV|„o z • NNLO-SV|pi z • NNLO-SV|„2

250 350 450 550 650 750 oc ^/Scцt [GeV]

Fig. 14. Hadronic NLO (dotted) and NNLO (solid) K factor as a function of ^scut.

Table 1

Total hadronic cross section at LO, NLO and NNLO-SV including top quark mass effects using | = 2mH and f(z) = z.

oh [fb] K(N)NLO XNNLO [%]

LO 22.7 - -

LO+NLO|po 36.4 1.60 -

LO+NLO|po +NNLO|po 39.7 1.75 0

LO+NLO|p0 +NNLO|p1 38.7 1.70 -2.5

LO+NLO|p0 +NNLO|p2 40.5 1.78 +2.0

as a function of V?cut. For comparison also the NLO result from Fig. 4 is shown as dotted curve using the Mt ^^ result. One observes that the various p orders lead to similar results for KHNNLO. Furthermore, there is a strong raise close to threshold which is due to the steeper behavior of the NNLO cross section as can be seen in the inlay of Fig. 11. For higher values of ^/scut, in particular for the total cross section, K^NLO approaches 1.7-1.8.

Results for the total cross section at LO, NLO and NNLO are summarized in Table 1 for i = 2mH. At NLO only p0 terms are included in the analysis whereas at NNLO p0, p1 and p2 terms are considered. In this way we can estimate the top mass effects of the NNLO term. Besides the cross section also the K factor is shown. At NNLO we use f(z) = z and we furthermore show

the relative deviation due to 1 /Mfn terms defined through

|pn - 8a

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%.

6. 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 Feyn-man 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 /Mt expansion for the virtual corrections. They are used to construct a soft-virtual approximation for the production cross section. In the limit Mt — to the effective-theory result of Ref. [14] is confirmed and 1/Mf and 1/M4 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^ut (a technical cut on the partonic center-of-mass energy). The curves including 1 /Mt 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 « 400 GeV. Although in this region the cross section is small it is of interest since there the cross section has a characteristic behavior. Furthermore, it is possible to use our result in this region as a benchmark for future calculations taking into account the exact dependence on Mt.

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 Mt —> to 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.

Acknowledgements

We would like to thank Kirill Melnikov for providing to us the analytic result for ¡\(4) in Eq. (23), for many useful discussions and for carefully reading the manuscript. This work is supported by the Deutsche Forschungsgemeinschaft through grant STE 945/2-1 and by KIT through its distinguished researcher fellowship program. Parts of this work were supported by the European Commission through contract PITN-GA-2012-316704 (HIGGSTOOLS). J.G. would like to express a special thanks to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and support.

References

[1] A. Djouadi, W. Kilian, M. Muhlleitner, P.M. Zerwas, Eur. Phys. J. C 10 (1999) 45, arXiv:hep-ph/9904287.

[2] U. Baur, T. Plehn, D.L. Rainwater, Phys. Rev. D 67 (2003) 033003, arXiv:hep-ph/0211224.

[3] U. Baur, T. Plehn, D.L. Rainwater, Phys. Rev. D 69 (2004) 053004, arXiv:hep-ph/0310056.

[4] M.J. Dolan, C. Englert, M. Spannowsky, J. High Energy Phys. 1210 (2012) 112, arXiv:1206.5001 [hep-ph].

[5] A. Papaefstathiou, L.L. Yang, J. Zurita, Phys. Rev. D 87 (2013) 011301, arXiv:1209.1489 [hep-ph].

[6] J. Baglio, A. Djouadi, R. Grober, M.M. Muhlleitner, J. Quevillon, M. Spira, J. High Energy Phys. 1304 (2013) 151, arXiv:1212.5581 [hep-ph].

[7] F. Goertz, A. Papaefstathiou, L.L. Yang, J. Zurita, J. High Energy Phys. 1306 (2013) 016, arXiv:1301.3492 [hep-ph].

[8] M. Gouzevitch, A. Oliveira, J. Rojo, R. Rosenfeld, G.P. Salam, V. Sanz, J. High Energy Phys. 1307 (2013) 148, arXiv:1303.6636 [hep-ph].

[9] E.W.N. Glover, J.J. van der Bij, Nucl. Phys. B 309 (1988) 282.

[10] T. Plehn, M. Spira, P.M. Zerwas, Nucl. Phys. B 479 (1996) 46;

T. Plehn, M. Spira, P.M. Zerwas, Nucl. Phys. B 531 (1998) 655 (Erratum), arXiv:hep-ph/9603205.

[11] S. Dawson, S. Dittmaier, M. Spira, Phys. Rev. D 58 (1998) 115012, arXiv:hep-ph/9805244.

[12] J. Grigo, J. Hoff, K. Melnikov, M. Steinhauser, Nucl. Phys. B 875 (2013) 1, arXiv:1305.7340 [hep-ph].

[13] F. Maltoni, E. Vryonidou, M. Zaro, J. High Energy Phys. 1411 (2014) 079, arXiv:1408.6542 [hep-ph].

[14] D. de Florian, J. Mazzitelli, Phys. Lett. B 724 (2013) 306, arXiv:1305.5206 [hep-ph].

[15] D. de Florian, J. Mazzitelli, Phys. Rev. Lett. 111 (2013) 201801, arXiv:1309.6594 [hep-ph].

[16] J. Grigo, K. Melnikov, M. Steinhauser, Nucl. Phys. B 888 (2014) 17, arXiv:1408.2422 [hep-ph].

[17] D.Y. Shao, C.S. Li, H.T. Li, J. Wang, arXiv:1301.1245 [hep-ph].

[18] D. de Florian, J. Mazzitelli, arXiv:1505.07122 [hep-ph].

[19] C. Anastasiou, K. Melnikov, Nucl. Phys. B 646 (2002) 220, arXiv:hep-ph/0207004.

[20] J. Grigo, J. Hoff, K. Melnikov, M. Steinhauser, PoS RADCOR 2013 (2013) 006, arXiv:1311.7425 [hep-ph], 2013.

[21] J. Grigo, J. Hoff, arXiv:1407.1617 [hep-ph].

[22] D. de Florian, J. Mazzitelli, J. High Energy Phys. 1212 (2012) 088, arXiv:1209.0673 [hep-ph].

[23] S. Catani, Phys. Lett. B 427 (1998) 161, arXiv:hep-ph/9802439.

[24] A.D. Martin, W.J. Stirling, R.S. Thorne, G. Watt, Eur. Phys. J. C 63 (2009) 189, arXiv:0901.0002 [hep-ph].

[25] K.A. Olive, et al., Particle Data Group Collaboration, Chin. Phys. C 38 (2014) 090001.

[26] G. Aad, et al., ATLAS CMS Collaborations, Phys. Rev. Lett. 114 (2015) 191803, arXiv:1503.07589 [hep-ex].

[27] C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, B. Mistlberger, Phys. Lett. B 737 (2014) 325, arXiv:1403.4616 [hep-ph].

[28] F. Herzog, B. Mistlberger, arXiv:1405.5685 [hep-ph].

[29] P. Nogueira, J. Comput. Phys. 105 (1993) 279.

[30] J. Kuipers, T. Ueda, J.A.M. Vermaseren, J. Vollinga, Comput. Phys. Commun. 184 (2013) 1453, arXiv:1203.6543 [cs.SC].

[31] R. Harlander, T. Seidensticker, M. Steinhauser, Phys. Lett. B 426 (1998) 125, arXiv:hep-ph/9712228.

[32] T. Seidensticker, arXiv:hep-ph/9905298.

[33] M. Steinhauser, Comput. Phys. Commun. 134 (2001) 335, arXiv:hep-ph/0009029.

[34] A.V. Smirnov, V.A. Smirnov, arXiv:1302.5885 [hep-ph].

[35] A.V. Smirnov, Comput. Phys. Commun. 189 (2014) 182, arXiv:1408.2372 [hep-ph].

[36] T.G. Birthwright, E.W.N. Glover, P. Marquard, J. High Energy Phys. 0409 (2004) 042, arXiv:hep-ph/0407343.

[37] R.K. Ellis, G. Zanderighi, J. High Energy Phys. 0802 (2008) 002, arXiv:0712.1851 [hep-ph].

[38] F. Chavez, C. Duhr, J. High Energy Phys. 1211 (2012) 114, arXiv:1209.2722 [hep-ph].

[39] A.B. Goncharov, Math. Res. Lett. 5 (1998) 497, arXiv:1105.2076 [math.AG].

[40] A.V. Smirnov, Comput. Phys. Commun. 185 (2014) 2090, arXiv:1312.3186 [hep-ph].

[41] J. Hoff, Methods for multiloop calculations and Higgs boson production at the LHC, Dissertation, KIT, 2015.

[42] J.A.M. Vermaseren, arXiv:math-ph/0010025.

[43] T. Gehrmann, T. Huber, D. Maitre, Phys. Lett. B 622 (2005) 295, arXiv:hep-ph/0507061.

[44] A. Pak, M. Rogal, M. Steinhauser, J. High Energy Phys. 1109 (2011) 088, arXiv:1107.3391 [hep-ph].