NLO predictions for Higgs boson pair production with full top quark mass dependence matched to parton showers

TL;DR

This paper presents the first NLO QCD predictions for Higgs boson pair production that include full top quark mass dependence, matched with parton showers, providing more accurate theoretical tools for experimental analysis.

Contribution

It combines NLO matrix elements with parton showers for di-Higgs production including full top mass dependence, implemented in POWHEG-BOX and MadGraph5_aMC@NLO frameworks.

Findings

Shower effects are significant for certain observables like Higgs pair transverse momentum.

Differences between full NLO and HEFT approximation are larger than shower effects in distribution tails.

The results reduce theoretical uncertainties in Higgs pair production predictions.

Abstract

We present the first combination of NLO QCD matrix elements for di-Higgs production, retaining the full top quark mass dependence, with a parton shower. Results are provided within both the POWHEG-BOX and MadGraph5_aMC@NLO Monte Carlo frameworks. We assess in detail the theoretical uncertainties and provide differential results. We find that, as expected, the shower effects are relatively large for observables like the transverse momentum of the Higgs boson pair, which are sensitive to extra radiation. However, these shower effects are still much smaller than the differences between the Born-improved HEFT approximation and the full NLO calculation in the tails of the distributions.

Equations

$$x$$

$$c_{\theta }$$

$${\mathcal{V}}_b$$

$$\mathcal{N}$$

$$\displaystyle\begin{split}\mathcal{V}_{\rm fin}(\mu_{r})=&\frac{2\pi}{\alpha_{s}(\mu_{r})}\Big{(}\mathcal{V}_{b}+\mathbf{I}\otimes\mathcal{B}\Big{)}(\mu_{r})\\ &-\mathcal{B}(\mu_{r})\Big{(}C_{A}\ln^{2}\Big{(}\frac{\mu_{r}^{2}}{\hat{s}}\Big{)}+\beta_{0}\ln\Big{(}\frac{\mu_{r}^{2}}{\hat{s}}\Big{)}+\beta_{0}+2K_{g}-\frac{2\pi^{2}}{3}C_{A}\Big{)}.\end{split}$$

$${\mathcal{V}}_{\mathrm{fin}}({\mu }_r)={\mathcal{V}}_{\mathrm{fin}}({\mu }_0)\cdot \frac{\mathcal{B}({\mu }_r)}{\mathcal{B}({\mu }_0)}+C_A\mathcal{B}({\mu }_r)\left({\ln }^2\left(\frac{{\mu }_0^2}{\hat{s}}\right)-{\ln }^2\left(\frac{{\mu }_r^2}{\hat{s}}\right)\right).$$

$$R_{\mathrm{sing}}=$$

$$R_{\mathrm{reg}}=$$

$$F=\frac{h^2}{(p_T^{\mathrm{hh}}{)}^2+h^2}.$$

Full text

††preprint:

CERN-TH-2017-069

MPP-2017-41

NIKHEF-2017-020

1 Introduction

Exploring the Higgs sector is one of the major goals for the next phases of LHC experiments. In particular, the form of the Higgs potential as predicted by the Standard Model (SM) needs to be confirmed. While one important parameter of the potential, the Higgs boson mass, has been measured already to an impressive accuracy, the Higgs boson self-coupling is still only very weakly constrained. The latter can be measured for example via Higgs boson pair production in gluon fusion, which is the dominant production mechanism of Higgs boson pairs. However, the cross section is about $1000$ times smaller than that for single Higgs production, which makes the measurement very challenging even with the high luminosity upgrade of the LHC. This fact on the other hand makes this channel very interesting for New Physics searches, as the delicate cancellations between different contributions which happen in the SM are altered in most New Physics models, leading to potentially large effects.

At the LHC, the decay channel $HH\to b\bar{b}\gamma\gamma$ has so far led to the most stringent limit on the di-Higgs production cross section of $\sigma/\sigma_{SM}\leq 19$ in CMS CMS-PAS-HIG-17-008 , while the ATLAS collaboration achieved the most restrictive upper bound of $\sigma/\sigma_{SM}\leq 29$ in the $b\bar{b}b\bar{b}$ decay channel ATLAS-CONF-2016-049 .

A previous combination of various decay channels measured in the ATLAS detector led to $\sigma/\sigma_{SM}\leq 70$ Aad:2015xja . The CMS collaboration also produced new limits for resonant and non-resonant Higgs boson pair production in the $b\bar{b}VV$ channel CMS-PAS-HIG-17-006 .

On the theory side, the leading order calculation of Higgs boson pair production in gluon fusion, which proceeds via heavy quark loops, has been performed in Refs. Eboli:1987dy ; Glover:1987nx ; Plehn:1996wb . Higher order corrections were for a long time available only within the Higgs Effective Field Theory (HEFT) approximation, where the NLO corrections are calculated in the $m_{t}\to\infty$ limit, leading to point-like effective couplings of gluons to Higgs bosons. In Ref. Dawson:1998py NLO corrections were calculated in the so-called “Born-improved HEFT” approximation, where the basic HEFT result is rescaled by a factor $B_{FT}/B_{HEFT}$, $B_{FT}$ denoting the leading order matrix element squared in the full theory.

In Refs. Frederix:2014hta ; Maltoni:2014eza , an approximation called “FTapprox” was introduced, which contains the full top quark mass dependence in the real radiation, while the virtual part is calculated in the HEFT approximation and rescaled at the event level by the re-weighting factor $B_{FT}/B_{HEFT}$.

In addition, the HEFT results at NLO and NNLO have been improved by an expansion in $1/m_{t}^{2\rho}$ in Refs. Grigo:2013rya ; Grigo:2014jma ; Grigo:2015dia ; Degrassi:2016vss , with $\rho^{\rm{max}}=6$ at NLO, and $\rho^{\rm{max}}=2$ for the soft-virtual part at NNLO  Grigo:2015dia .

The NNLO QCD corrections in the heavy top limit have been computed in Refs. deFlorian:2013uza ; deFlorian:2013jea ; Grigo:2014jma ; deFlorian:2016uhr , and they have been supplemented by an expansion in $1/m_{t}^{2}$ in Ref. Grigo:2015dia and by resummation, at NLO+NNLL in Ref. Shao:2013bz and at NNLO+NNLL in Ref. deFlorian:2015moa , leading to K-factors of about 1.2 relative to the Born-improved HEFT result.

Very recently, the full NLO corrections, including the top quark mass dependence also in the virtual two-loop amplitudes, have been calculated Borowka:2016ehy ; Kerner:2016msj and compared to previous approximations for various observables Borowka:2016ypz . The full NLO calculation was supplemented by NLL resummation in Ref. Ferrera:2016prr .

Numerous phenomenological studies of Higgs boson pair production have been performed both within and beyond the SM Dolan:2012rv ; Baglio:2012np ; Goertz:2013kp ; Gouzevitch:2013qca ; Barr:2013tda ; Dolan:2013rja ; Barger:2013jfa ; Li:2013flc ; deLima:2014dta ; Slawinska:2014vpa ; Buschmann:2014sia ; Azatov:2015oxa ; Dicus:2015yva ; Papaefstathiou:2015iba ; Grober:2015cwa ; Lu:2015jza ; Dawson:2015oha ; Ghezzi:2015vva ; Dolan:2015zja ; DallOsso:2015aia ; Dawson:2015haa ; Cao:2015oaa ; Cao:2015oxx ; Behr:2015oqq ; Li:2016nrr ; Kanemura:2016lkz ; Agostini:2016vze ; Grober:2016wmf ; Banerjee:2016nzb ; Huang:2017jws . Further, it also has been suggested recently to obtain constraints on the Higgs boson self-coupling from electroweak corrections to single Higgs boson production Gorbahn:2016uoy ; Degrassi:2016wml ; Bizon:2016wgr .

The studies of Higgs boson pair production mentioned above usually had at least one of the following drawbacks: either they are based on leading order matrix elements, while including the full top quark mass dependence, or the matrix elements include higher orders in QCD but have been performed within the infinite-top-mass approximation, which is known to fail at scales where the top quark loops are resolved Baur:2002rb ; Borowka:2016ehy ; Borowka:2016ypz .

Results for Higgs boson pair production merged to $HH+1$ jet matrix elements at leading order, with full top and bottom quark mass dependence, matched to a parton shower within HERWIG++, have been presented in Ref. Maierhofer:2013sha . The “FTapprox” Maltoni:2014eza calculation includes the matching of di-Higgs production to a parton shower Frederix:2014hta keeping the full top quark mass dependence in the real radiation, while the virtual part is calculated in the Born-improved HEFT approximation.

In this paper, we present the first combination of the full NLO calculation, including the full top quark mass dependence at two loops, with a parton shower, within both the POWHEG-BOX Frixione:2007vw ; Alioli:2010xd and the MadGraph5_aMC@NLO framework Alwall:2014hca ; Hirschi:2015iia . This allows us to compare the POWHEG Frixione:2007vw and MC@NLO Frixione:2002ik matching schemes while using the same Pythia 8 Sjostrand:2007gs ; Sjostrand:2014zea shower in both cases. We also investigate the PDF and scale uncertainties and calculate observables like the di-Higgs $p_{T}$ spectrum, $p_{T}^{hh}$, where fixed order NLO calculations cannot give a satisfactory description at low $p_{T}^{hh}$. Further, we discuss the possibility to infer the leading contribution of a full NNLO calculation from the showered results, based on a comparison of the NNLO calculation in the HEFT approximation deFlorian:2016uhr with the showered results.

2 Details of the calculation

In this section we present the details of the implementation of the calculation within the POWHEG-BOX Monte Carlo program. The results from MG5_aMC@NLO presented in the next sections are based on a similar implementation. Both codes use the same grid for the virtual two-loop amplitude discussed in Sec. 2.1. Further details about the calculation based on Born-improved HEFT and FTapprox within MG5_aMC@NLO already have been published elsewhere Frederix:2014hta ; Maltoni:2014eza .

In order to allow for comparisons and cross checks, we implemented both the effective theory as well as the full SM amplitudes at NLO. This allows to run the code in four different modes by changing the flag mtdep in the POWHEG-BOX run card. The possible choices and the corresponding calculation, as presented in the previous section, are the following:

mtdep=0**:**

computation using basic HEFT,

mtdep=1**:**

computation using Born-improved HEFT,

mtdep=2**:**

computation in the approximation FTapprox (full mass dependence in the Born and in the real radiation, Born-improved HEFT for the virtual part),

mtdep=3**:**

computation in the full SM.

The corresponding modes are also available in MG5_aMC@NLO.

The leading order amplitude in the full theory and all the amplitudes in the HEFT were implemented analytically, whereas the one-loop real radiation contribution and the two-loop virtual amplitudes in the full SM rely on numerical or semi-numerical codes. Since the virtual two-loop amplitudes in the full theory are computed keeping the Higgs bosons on-shell, we assume a vanishing Higgs boson width in all the modes listed above. Higgs boson decays can be computed in the narrow width approximation by the parton shower. In the next sections we give some more details about our implementation. We will use the term “HEFT approximation”, or simply “HEFT”, for basic HEFT, while results in the Born-improved HEFT will be denoted by “B-i. HEFT”.

2.1 Virtual two-loop amplitudes

For the virtual two-loop amplitudes, we have used the results of the calculation presented in Refs. Borowka:2016ehy ; Borowka:2016ypz , which is based on an extension of the program GoSam Cullen:2014yla to two loops Jones:2016bci , using also Reduze 2 vonManteuffel:2012np and SecDec 3 Borowka:2015mxa .

The values for the Higgs boson and top quark masses have been set to $m_{h}=125$ GeV and $m_{t}=173$ GeV, such that the two-loop amplitudes only depend on two independent variables, the Mandelstam invariants $\hat{s}$ and $\hat{t}$. We have constructed a grid in these variables together with an interpolation framework, such that an external program can call the virtual two-loop amplitude at any phase space point without having to do costly two-loop integrations.

In more detail, we first transform the Mandelstam invariants $\hat{s}$ and $\hat{t}$ to new variables

[TABLE]

where $f$ can, in principle, be any strictly increasing function. Setting $f(\beta)$ according to the cumulative distribution function of the phase space points used in our original calculation, we obtain a nearly uniform distribution of these points in the $(x,c_{\theta})$ unit square. Instead of a direct interpolation of the phase space points, we chose to apply a two-step procedure: First, we generate a regular grid with a fixed grid spacing in the variables $x$ and $c_{\theta}$, where we estimate the result at each grid point applying a linear interpolation of our original results in the vicinity of the specific grid point. In a second step, we apply Clough-Tocher interpolation CloughTocher as implemented in the python SciPy package scipy . Applying this procedure reduces the size of interpolation artefacts, which we obtain due to the numerical uncertainty of our two-loop results. In Fig. 1 we test how omitting input data points influences the results of the grid interpolation. Removing 20% of the input data points changes the interpolation results by less than 0.25% for 70% of the tested points. Differences larger than 5% are obtained for 6% of the results.

We should point out that the grid is constructed from a sample of phase space points which is based on runs at $\sqrt{s}=14$ TeV. Therefore, even though the grid is not explicitly dependent on the centre-of-mass energy, one should be aware of the fact that for runs at e.g. 100 TeV, the grid may not be reliable for points with large $\hat{s}$ due to a lack of statistics in this region upon construction.

In our original calculation Borowka:2016ehy ; Borowka:2016ypz , we used Catani-Seymour dipole subtraction Catani:1997vz for the real radiation. The finite combination of the renormalized virtual amplitude $\mathcal{V}_{b}$ with the Catani-Seymour I-operator can be straightforwardly converted into the quantity ${\cal V}_{\rm{fin}}$ of Refs. Frederix:2009yq ; Alioli:2010xd , defined by

[TABLE]

For this process the colour-correlated Born squared amplitudes $\mathcal{B}_{12}$ and $\mathcal{B}_{21}$ are equal and we have $a=-2C_{A}$ and $c_{12}=c_{21}=-\beta_{0}/2-C_{A}\ln\left(\mu_{r}^{2}/\hat{s}\right)$. In the POWHEG-BOX and MG5_aMC@NLO frameworks the arbitrary scale $Q$ is chosen as $\mu_{r}$. Specifically, we obtain

[TABLE]

The grid evaluates $\mathcal{V}_{\rm fin}$ at the scale $\mu_{0}=\sqrt{\hat{s}}/2$ and the results for an arbitrary scale can be obtained from the relation

[TABLE]

The Born amplitude $\mathcal{B}$ and the colour-correlated Born amplitudes $\mathcal{B}_{ij}$ in (3) are evaluated in $D=(4-2\epsilon)$ dimensions using conventional dimensional regularization (CDR). As all formulas for the soft contributions and the collinear remnants used in the POWHEG-BOX and MG5_aMC@NLO are computed in the $\overline{\rm{MS}}$ scheme, using CDR, constructing $\mathcal{V}_{\rm{fin}}$ according to (5) ensures the treatment is consistent with that implemented in the Monte Carlo programs.

2.2 Real radiation and parton shower matching

The real radiation matrix elements in the full SM were implemented using the interface Luisoni:2013cuh between GoSam Cullen:2011ac ; Cullen:2014yla and the POWHEG-BOX Frixione:2007vw ; Alioli:2010xd , modified accordingly to compute the real corrections instead of the virtual ones. The one-loop real amplitudes we generated with the new version 2.0 of GoSam Cullen:2014yla , that uses QGRAF Nogueira:1991ex , FORM Kuipers:2012rf and Spinney Cullen:2010jv for the generation of the Feynman diagrams, and offers a choice from Samurai Mastrolia:2010nb ; vanDeurzen:2013pja , golem95C Binoth:2008uq ; Cullen:2011kv ; Guillet:2013msa and Ninja vanDeurzen:2013saa ; Peraro:2014cba for the reduction. At run time the amplitudes were computed using Ninja vanDeurzen:2013saa ; Peraro:2014cba and OneLOop vanHameren:2010cp for the evaluation of the scalar one-loop integrals.

In order to avoid numerical instabilities in the one-loop real matrix elements in the limit where the additional parton becomes soft and/or collinear, a technical cut $p_{T}^{\mathrm{hh}}>10^{-3}$ $\mathrm{GeV}$ has been introduced. We carefully checked that the total cross section does not change significantly when varying the cut value.

Within MG5_aMC@NLO, the one-loop born and real amplitudes are computed using MadLoop Hirschi:2011pa , which in turn exploits CutTools Ossola:2007ax , Ninja Peraro:2014cba ; Hirschi:2016mdz or Collier Denner:2016kdg , together with an in-house implementation of the OpenLoops optimisation Cascioli:2011va . In MG5_aMC@NLO, the computation is based on event reweighting, as described in Frederix:2014hta ; Maltoni:2014eza . This functionality in MG5_aMC@NLO has since been automated as documented in Mattelaer:2016gcx . In practice, the HEFT is used to generate the events at NLO which are then reweighted to introduce the full top quark mass dependence of the one- and two-loop amplitudes.

In order to study the phenomenological impact of the two different matching schemes implemented in the POWHEG-BOX and in MG5_aMC@NLO, we compare the two results using the same parton shower. In both cases we use Pythia 8.2 with the same settings to produce showered events from both the POWHEG-BOX and the MG5_aMC@NLO results at LHE level. This means that the differences in the distributions produced by POWHEG and MG5_aMC@NLO respectively, will only be due to the corresponding matching schemes.

3 Results

In this section we present phenomenological results and compare predictions at different levels. We start presenting some consistency checks at the fixed order level and at the Les Houches event level (LHE), i.e. after the first hard emission is generated according to the POWHEG method. To assess the impact of the parton shower and estimate its capacity to include approximate higher order effects, in Section 3.3 we compare results in the basic HEFT approximation at NLO+PS with the NNLO predictions from Reference deFlorian:2016spz . Finally, in Section 3.4 NLO and NLO+PS results in the full SM are presented.

All the results we computed using the PDF4LHC15_nlo_30_pdfas Butterworth:2015oua ; CT14 ; MMHT14 ; NNPDF parton distribution functions interfaced to our codes via LHAPDF Buckley:2014ana , along with the corresponding value for $\alpha_{s}$. The masses of the Higgs boson and the top quark have been set, as in the virtual amplitude, to $m_{h}=125$ GeV, $m_{t}=173$ GeV, respectively, whereas their widths have been set to zero. As already mentioned in the previous section, we consider on-shell Higgs bosons and leave the analysis of more exclusive final states, stemming from Higgs boson decays, to future studies. Jets are clustered with the anti-$k_{\scriptscriptstyle\rm T}$ algorithm Cacciari:2008gp as implemented in the Fastjet package Cacciari:2005hq ; Cacciari:2011ma , with jet radius $R=0.4$ and a transverse momentum greater than $p_{T,min}^{\rm{jet}}=20~{}\mathrm{GeV}$. The theoretical scale uncertainty is estimated by varying the factorization scale $\mu_{F}$ and the renormalization scale $\mu_{R}$. The scale variation bands are obtained by computing the envelopes of a 7-point scale variation around the central scale $\mu_{0}=m_{\mathrm{hh}}/2$, with $\mu_{R,F}=c_{R,F}\,\mu_{0}$, where $c_{R},c_{F}\in\{2,1,0.5\}$. The extreme variations $(c_{R},c_{F})=(2,0.5)$ and $(c_{R},c_{F})=(0.5,2)$ have been omitted.

We also have varied the PDFs using the $30$ error PDFs contained in the

PDF4LHC15_nlo_30_pdfas set and found that the uncertainty due to PDF variations never exceeds 6% and therefore is well below the scale variation uncertainty. The uncertainty bands shown on our results originate from scale variations only.

We should point out that we switched off the hadronisation and the multiple interactions in the parton shower. A detailed phenomenological study including various Higgs boson decay channels as well as hadronisation effects will be left to a subsequent publication.

3.1 Comparison with previous NLO results

A very strong consistency check of the new implementation, which allows to test at the same time the real amplitudes and their stability, the implementation of the grid for the virtual two-loop amplitude and all the various parts of the code relevant for the NLO fixed order computation is a comparison with the previous NLO results computed in Refs. Borowka:2016ehy ; Borowka:2016ypz . A comparison at the level of individual phase space points shows that the grid interpolation slightly increases the numerical uncertainty associated to the virtual amplitude results, which were calculated with percent level precision in the previous publications. At the level of differential distributions we found excellent agreement, not only for the full NLO results, but also for the various other approximations available and for uncertainties related to scales variation.

In Fig. 2 we show a comparison between the NLO predictions obtained with the POWHEG-BOX generator and the original ones from Ref. Borowka:2016ypz . In the MG5_aMC@NLO case, where only showered results are available, we made a similar validation plot for the $m_{\mathrm{hh}}$ distribution, as it is insensitive to shower effects.

We should mention at this point that beyond $p_{T}^{\mathrm{h}}\sim 650~{}\mathrm{GeV}$, a systematic bias stemming from lack of statistics in the grid starts to develop. As a consequence, the results obtained with the grid will be systematically below the “true” results. The difference is within the statistical uncertainty up to about $p_{T}^{\mathrm{h}}\sim 750~{}\mathrm{GeV}$, and increases to about 20% at $p_{T}^{\mathrm{h}}\sim 1~{}\mathrm{TeV}$ (and $m_{\mathrm{hh}}\sim 2~{}\mathrm{TeV}$).

3.2 Comparisons at the level of Les Houches event files

Before presenting results for NLO predictions matched to the parton shower, we show comparisons of NLO curves with results at the Les Houches event (LHE) level, after the first hard emission is weighted with the Sudakov factor according to the POWHEG method. Even though the LHE level predictions still need to be showered, such a comparison allows to test the implementation and, once the results are fully showered, to disentangle the impact of the shower from the one due to the POWHEG exponentiation. For observables which are inclusive in the extra radiation, the fixed order NLO and LHE level predictions should be in perfect agreement. We show the level of agreement between the NLO and LHE curves for the Higgs-boson pair invariant mass $m_{hh}$ in Fig. 3, where a comparison is shown for predictions in the HEFT approximation and the full SM.

For observables which are directly sensitive to soft gluon radiation, like the $p_{T}^{\mathrm{hh}}$ distribution in the limit $p_{T}^{\mathrm{hh}}\to 0$, one instead expects to observe Sudakov suppression in the soft region. This can be seen in Fig. 4, where we can clearly see the suppression in the region where the fixed order NLO results become unreliable.

We also note that the LHE predictions are enhanced in the high transverse momentum region compared to the NLO curve. This is due to subleading contributions in the exponential, which in the case of large radiative corrections can become sizable, in particular for observables like $p_{T}^{\mathrm{hh}}$, where NLO is the first non-trivial order to describe the distribution. Analogous effects have already been observed in several other similar processes with large K-factors Alioli:2008tz ; Alioli:2009je ; Alioli:2016xab ; Hespel:2015zea , and we refer the interested reader to Refs. Alioli:2008tz ; Alioli:2009je for more details. We have explored the possibility to limit the amount of hard radiation which is exponentiated by changing the hdamp parameter in POWHEG. We recall that this allows to divide the contributions of the real radiation $R$ which are exponentiated in the Sudakov factor into a singular part $R_{\mathrm{sing}}$ and a regular part $R_{\mathrm{reg}}$, as follows:

[TABLE]

where the transition function $F$ is chosen to be

[TABLE]

In Fig. 5 we compare the default POWHEG setting, $h={\tt hdamp}=\infty$, with predictions where we use ${\tt hdamp}=250~{}\mathrm{GeV}$. The left plot shows predictions in the HEFT, whereas on the right we show results in the full SM. We observe that in both cases above $500~{}\mathrm{GeV}$ the LHE curve with ${\tt hdamp}=250~{}\mathrm{GeV}$ reproduces the NLO results as expected. It is interesting to study how this additional source of theoretical uncertainty is affecting other observables, especially those for which our predictions are NLO accurate. To understand this better, in Fig. 6 we show a similar comparison for $m_{\mathrm{hh}}$ (left) and the transverse momentum of a (randomly chosen) Higgs boson $p_{T}^{\mathrm{h}}$ (right), with full top quark mass dependence. The $m_{\mathrm{hh}}$ observable is completely insensitive to additional radiation, and for this reason it is unaffected by a modification of the hdamp factor. This is not true for $p_{T}^{\mathrm{h}}$, which is sensitive to the recoil against additional jet activity. For this reason we observe deviations between the NLO predictions and the LHE-level curves, the latter becoming slightly larger for harder transverse momenta. The predictions for ${\tt hdamp}=250$ are in general closer to the NLO ones over the whole kinematical range of $p_{T}^{\mathrm{h}}$. We stress however that, contrary to $p_{T}^{\mathrm{hh}}$, where the differences between the predictions for ${\tt hdamp}=\infty$ and the one for ${\tt hdamp}=250~{}\mathrm{GeV}$ reach 80% above 500 $\mathrm{GeV}$, for $p_{T}^{\mathrm{h}}$ the differences are at the 10-15% level, i.e. well within the scale uncertainties.

Since the uncertainty related to the value for hdamp is very tightly related to the POWHEG way of matching NLO to the parton-shower, it is important to compare these predictions with other matching schemes. We will comment more on this aspect in Section 3.4, where we compare NLO+PS predictions obtained with POWHEG and MG5_aMC@NLO.

3.3 Discussion of NNLO effects

We mentioned in the previous section that the enhancement in the tail of $p_{T}^{\mathrm{hh}}$ in the prediction at the LHE-level is due to subleading contributions in the Sudakov factor, which are intrinsically taken into account in the matching à la POWHEG. As already pointed out and discussed in Alioli:2008tz , it is therefore interesting to compare NLO+PS predictions obtained with POWHEG to the full NNLO predictions, if they are available. This is indeed the case if we restrict ourselves to results in the basic HEFT, for which differential NNLO results were computed in Ref. deFlorian:2016uhr 111We are grateful to Javier Mazzitelli for providing us the NNLO predictions shown in the comparisons of this section..

The comparison is of course meaningful only for those observables which are LO accurate in our NLO calculation and which therefore are not sensitive to the additional two-loop virtual corrections included in the HEFT NNLO predictions. In Fig. 7 we consider again the transverse momentum of the Higgs boson pair and compare the NNLO results with two different LHE-level predictions from POWHEG. On the left we keep the default setting in which ${\tt hdamp}=\infty$, on the right we set ${\tt hdamp}=250~{}\mathrm{GeV}$. In the former plot we observe a good agreement of the LHE-level curve with ${\tt hdamp}=\infty$ with the NNLO predictions in the transverse momentum range between $200~{}\mathrm{GeV}$ and $400~{}\mathrm{GeV}$. While the LHE-level result flattens out around $250~{}\mathrm{GeV}$, the NNLO result decreases slightly for larger $p_{T}^{\mathrm{hh}}$. The two theory uncertainty bands due to scale variation however largely overlap. The plot on the right shows instead that, by limiting the amount of real radiation in the Sudakov factor, the LHE-level prediction falls onto the NLO result at high $p_{T}$, and therefore cannot reproduce the NNLO behaviour.

As a further step, we can assess the impact of the parton shower, by analyzing the same observable with NLO+PS predictions showered with Pythia 8. Figure 8 shows an analogous comparison, where the NLO curves with and without shower are plotted against the NNLO predictions. We observe that the shower has a large effect on the tail of the $p_{T}^{\mathrm{hh}}$ distribution, such that the NNLO curve lies between the NLO+PS and the NLO fixed order curve for ${\tt hdamp}=\infty$. On the other hand, for ${\tt hdamp}=250\,\mathrm{GeV}$, the NLO+PS result by construction is closer to the NLO fixed order result. We should point out however that these considerations within the basic HEFT approximation may not carry over analogously to the full calculation (where NNLO predictions are not available), because it is well known that the HEFT approximation does not have the correct scaling behaviour at large transverse momenta.

3.4 NLO plus parton shower matched results

We now compare fixed order NLO results to our default POWHEG results, where we use hdamp =250 and the Pythia 8 shower. In Fig. 9 we show the Higgs boson pair invariant mass distribution and the transverse momentum distribution of a randomly chosen Higgs boson for both the fixed order and the showered calculation. As is to be expected, the invariant mass distribution is rather insensitive to the parton shower. In Fig. 10 the $p_{T}$-distributions of the harder and softer Higgs boson are shown.

We should mention at this point that the distributions of the “harder” ($p_{T}^{\mathrm{h}_{1}}$) and “softer” ($p_{T}^{\mathrm{h}_{2}}$) Higgs boson, calculated at fixed (NLO) order, are somewhat infrared sensitive if no cuts are placed on the Higgs boson transverse momenta. The reason is that, if the transverse momenta $p_{T}^{\mathrm{h}_{1}}$ and $p_{T}^{\mathrm{h}_{2}}$ are very close to each other, the available phase space for the extra radiation in the real corrections is severely restricted, leading to large logarithms which are not sufficiently balanced by the $2\to 2$ contributions. To illustrate this fact, we consider the total cross section as a function of $\Delta$, with the kinematic requirements $p_{T}^{\mathrm{h}_{1}}\geq\Delta,p_{T}^{\mathrm{h}_{2}}\geq 0$. The cross section shows an unphysical behaviour as $\Delta\to 0$, see Fig. 11: the total cross section as a function of $\Delta$ peaks around $\Delta=14~{}\mathrm{GeV}$ and then decreases for smaller values of $\Delta$, even though the available phase space for $p_{T}^{\mathrm{h}_{1}}$ is larger. This behaviour is an artifact of the fixed order calculation and is the reason why “symmetric cuts” (i.e. the same $p_{T,min}$ values for both final state particles in a $2\to 2$ calculation at NLO) should be avoided. For a more detailed discussion of this point we refer to Refs. Catani:1997xc ; Frixione:1997ks ; Fontannaz:2001nq . Here we only note that this is the reason why, with “symmetric” cuts $p_{T,min}^{\mathrm{h}_{1}}=p_{T,min}^{\mathrm{h}_{2}}=0$ and fine binning, the first bin(s) of the $p_{T}^{\mathrm{h}_{1}}$ distribution are negative at fixed order, while this behaviour is cured by the Sudakov factor, so it is absent in the LHE level and showered results.

Fig. 12 displays the transverse momentum distributions of the Higgs boson pair and of the (leading) jet. As discussed already in the context of Fig 5, the $p_{T}^{\mathrm{hh}}$ distribution diverges at fixed order for $p_{T}^{\mathrm{hh}}\to 0$, while the showered result is able to provide reliable predictions in the low $p_{T}^{\mathrm{hh}}$ region. We notice that the scale variation band is reduced in the showered result compared to the fixed order calculation. The scale uncertainties on the fixed order results are particularly large for these distributions as they are – except for the first bin – determined by the $2\to 3$ kinematics, which is described only at leading order accuracy by our calculation.

In Fig. 13 we show the difference in azimuthal angle, $\Delta\Phi^{hh}$, and the radial separation, $\Delta R^{\mathrm{hh}}=\sqrt{(\eta_{1}-\eta_{2})^{2}+(\Phi_{1}-\Phi_{2})^{2}}$, of the two Higgs bosons. We see that the unphysical behaviour for $\Delta\Phi^{hh}\to 0$ of the fixed order result is cured by the Sudakov form factor, and again the scale uncertainties of the fixed order calculation are relatively large because the tail of the distribution is predicted at the first non-trivial order. In the $\Delta R^{\mathrm{hh}}$ distribution, we observe that the shower populates the region $\Delta R^{\mathrm{hh}}<\pi$, which at fixed order is given by the $2\to 3$ component only.

Fig. 14 compares the predictions obtained with the Pythia 6 shower to the Pythia 8 results both in the basic HEFT approximation and in the full SM. It is instructive to make this comparison for hdamp=$\infty$ (left column) as well as for hdamp=250 (right column). In the basic HEFT approximation the differences between Pythia 6 and Pythia 8 are small, and setting hdamp to a finite value restores the agreement between the NLO and the NLO+PS curves at large transverse momentum. The latter is also true in the full SM. However, in the full SM, the difference between Pythia 6 and Pythia 8 is much larger, Pythia 8 showing a considerably harder spectrum in the tail of the $p_{T}^{\mathrm{hh}}$ distribution.

To conclude this section, in Fig. 15(a) we compare the full results with calculations where the underlying matrix elements are based on two approximations, either FTapprox or Born-improved HEFT. All matrix elements are combined with the same Pythia 8 shower. In order to assess the effect of the parton shower on the various approximations, we also show the fixed order results in Fig. 15(b). The broad features of these approximations remain unchanged after showering, however, as the showered results have smaller scale uncertainties, the differences between these approximations are actually enhanced if a parton shower is attached.

3.4.1 Comparison between POWHEG and MG5_aMC@NLO

In this section we compare the POWHEG results with results from MG5_aMC@NLO, the latter being based on the same grid in the invariants $\hat{s}$ and $\hat{t}$ for the virtual two-loop corrections as the POWHEG results, and based on the same Pythia 8 shower. Therefore the differences between the results can be attributed to differences in the matching scheme.

In Figs. 16 to 18 we show POWHEG results for two different values of hdamp compared to MG5_aMC@NLO results. While for the $p_{T}^{\mathrm{h}_{1}}$ and $p_{T}^{\mathrm{h}_{2}}$ distributions the differences are mostly small, they are, as to be expected, more pronounced for the distributions where the shower populates kinematic regions which are predicted at the first non-trivial order by the NLO fixed order calculation. Focusing on the comparison between the POWHEG curve with hdamp$=250$, which is our default, we can say that in general the two predictions agree well within the scale and statistical uncertainties. The small $\Delta R^{\mathrm{hh}}$ region, on the right of Fig. 18, shows the largest differences, which is not surprising as it is dominated by (multi-)jet events. We should also mention that the curve for hdamp$=250$ in these figures is close to the NLO curves by construction, as can be seen by comparing to the fixed order results shown in the previous subsection.

In Fig. 19 we vary the shower starting scale $Q_{sh}$ in MG5_aMC@NLO by a factor of two around the default value. In the latest version of MG5_aMC@NLO (version 2.5.3 onwards), the shower starting scale is picked with some probability distribution to be in the interval shower_scale_factor $\times\,[0.1\,H_{T}/2,H_{T}/2]$ with $H_{T}$ computed with Born kinematics, therefore to perform the scale variation we set the shower_scale_factor in the run card to 0.5, 1 and 2.

The $m_{\mathrm{hh}}$ distribution can be considered as a control plot to demonstrate that, as expected, this has no effect on the $m_{\mathrm{hh}}$ distribution. In contrast, in the $p_{T}^{\mathrm{hh}}$ distribution, the differences due to variations of the matching scale start to exceed the scale uncertainties towards larger $p_{T}^{\mathrm{hh}}$ values.

Because of the fact that for the $p_{T}^{\mathrm{hh}}$ distribution, the tail is predicted at the first non-trivial order, the effect of the shower on this distribution is rather large, exceeding a factor of two beyond $p_{T}^{\mathrm{hh}}\sim 300$ GeV, as shown in Fig. 20. However, as can also be seen from Fig. 20, the differences due to the shower are still much smaller than the difference between the full calculation and the Born-improved HEFT approximation, which is off by an order of magnitude for $p_{T}^{\mathrm{hh}}>500$ GeV. Fig. 20 also shows that FTapprox does a good job for this observable, as the tail of the $p_{T}^{\mathrm{hh}}$ distribution is determined by the real radiation. In the POWHEG case, the FTapprox curve still lies above the full result because the differences in the virtual part enter the $\bar{B}$ function in POWHEG, which determines the overall normalisation for the shower.

Finally, we compare in Fig. 21 the fixed order result to showered results using different values for hdamp in POWHEG and for the shower starting scale $Q_{\rm{sh}}$ in MG5_aMC@NLO. The new shower starting scale in MG5_aMC@NLO is picked in some interval with $H_{T}/2$ as its maximum as stated above, while the old shower starting scale was picked in the interval $[0.1\,\sqrt{\hat{s}},\sqrt{\hat{s}}]$. One can observe that with the new shower starting scale in MG5_aMC@NLO, $Q_{\rm{sh}}=Q_{\rm{def}}^{\rm{new}}$, the showered results match onto the fixed order curve at large values of $p_{T}^{\mathrm{hh}}$, while the latter is not the case for POWHEG with hdamp=$\infty$ and MG5_aMC@NLO with the old default shower starting scale.

4 Conclusions

We have presented the combination of the full NLO prediction for Higgs boson pair production, including the top quark mass dependence at two loops, with a parton shower. This has been implemented within two frameworks, POWHEG-BOX and MG5_aMC@NLO, using the same Pythia 8.2 shower in both cases. Individual phase-space points of the two-loop amplitude, which depends only on the two independent kinematic invariants $\hat{s}$ and $\hat{t}$ once the top-quark and Higgs boson masses are fixed, have been used to create a grid and combined with an interpolation framework, such that a value for the amplitude can be obtained at any phase space point without re-evaluating the loop integrals.

We find that the impact of the parton shower on the transverse momentum distribution of one Higgs boson, $p_{T}^{\mathrm{h}}$, is quite small and that the features of the various approximations that have appeared previously in the literature are preserved by the shower.

The impact of the shower on the $p_{T}^{\mathrm{hh}}$, $\Delta\Phi^{hh}$ and $\Delta R^{hh}$ distributions is fairly large, as these are the distributions where the tail is predicted at the first non-trivial order in the fixed order calculation. In the tail of the $p_{T}^{\mathrm{hh}}$ distribution, around $p_{T}^{\mathrm{hh}}\sim 400$ GeV, the showered NLO results are larger than the fixed order results by more than a factor of two, within both POWHEG and MG5_aMC@NLO. This feature is also present if Pythia 6 is used instead of Pythia 8, and if we vary the shower starting scale in MG5_aMC@NLO. However, the differences due to the shower in the $p_{T}^{\mathrm{hh}}$ distribution are still much smaller than the discrepancy between the showered full calculation and the showered Born-improved HEFT approximation, the latter overshooting the full result by an order of magnitude around $p_{T}^{\mathrm{hh}}\sim 400$ GeV, worsening towards higher $p_{T}^{\mathrm{hh}}$ values. As expected, the FTapprox results, which include the full mass dependence in the real radiation, behave very similar to the full calculation in the (real radiation dominated) tails of the distributions like $p_{T}^{\mathrm{hh}}$.

In summary, we observe that the inclusion of the full mass dependence in general has a more important impact on the distributions relevant to Higgs boson pair production than effects coming from different shower matching schemes, variations of the shower starting scales, different parton showers or different PDFs. A detailed study of hadronisation effects and Higgs boson decays will be performed in a subsequent publication.

The POWHEG version of the code for Higgs boson pair production developed for this work is publicly available in the POWHEG-BOX V2 package, under the User-Processes-V2/ggHH/ directory, and will become available also in the newer POWHEG-BOX RES version in the User-Processes-RES/ggHH/ folder. All the information can be found at the web page http://powhegbox.mib.infn.it. The implementation in MG5_aMC@NLO is not part of the public release yet, but the customised code can be obtained by contacting the authors.

We hope that making two-loop results available in the form of a grid included in public Monte Carlo programs, as done in this work, will open the door to further developments in this direction.

Acknowledgements

We would like to thank our collaborators Sophia Borowka, Nicolas Greiner, Johannes Schlenk and Tom Zirke for their contribution to the double Higgs fixed order NLO project and Simone Alioli for interesting discussions and suggestions concerning the POWHEG-BOX. We also thank Javier Mazzitelli for providing us the NNLO HEFT data for comparisons. We also would like to thank Fabio Maltoni and Paolo Torrielli for useful discussions and Carlo Oleari for comments on the manuscript. G.L. thanks the Max Planck Institute for Physics for hospitality during several phases of this work. This research was supported in part by the Research Executive Agency (REA) of the European Union under the Grant Agreement PITN-GA2012316704 (HiggsTools). E.V. acknowledges support from the Netherlands Foundation for Fundamental Research of Matter (FOM) under project number 15PR3224. We gratefully acknowledge support and resources provided by the Max Planck Computing and Data Facility (MPCDF).