Probing the trilinear Higgs boson coupling in di-Higgs production at NLO QCD including parton shower effects

TL;DR

This paper provides NLO QCD predictions for Higgs pair production with variable trilinear couplings, including full top mass dependence, and compares parton shower effects, aiding in probing Higgs self-interactions.

Contribution

It introduces a publicly available NLO QCD calculation framework for Higgs pair production with variable self-coupling, incorporating full top mass effects and parton shower comparisons.

Findings

Differential results at 14 TeV show sensitivity to anomalous couplings.

Parton shower effects differ between PYTHIA and HERWIG implementations.

The method allows studying variations of the top Yukawa coupling as well.

Abstract

We present results for Higgs boson pair production with variations of the trilinear Higgs boson self-coupling at next-to-leading order (NLO) in QCD including the full top quark mass dependence. Differential results at 14 TeV are presented, and we discuss the implications of anomalous trilinear couplings as well as differences between the PYTHIA 8.2 and HERWIG 7.1 parton showers in combination with POWHEG. The implementation of the NLO QCD calculation with variable Higgs boson self-coupling is made publicly available in the POWHEG-BOX-V2 Monte Carlo framework. A simple method for using the new implementation to study also variations of the top quark Yukawa coupling is described.

Tables

Table 1. Table 1 : Total cross sections for Higgs boson pair production at full NLO QCD. The given uncertainties are scale uncertainties.

${\lambda }_{\mathrm{BSM}}/{\lambda }_{\mathrm{SM}}$	${\sigma }_{\mathrm{NLO}}@13TeV$ [fb]	${\sigma }_{\mathrm{NLO}}@14TeV$ [fb]	${\sigma }_{\mathrm{NLO}}@27TeV$ [fb]	K-factor@14TeV
-1	116.71${}_{-14.3\%}{}^{}{}_{}{}^{+16.4\%}$	136.91${}_{-13.9\%}{}^{}{}_{}{}^{+16.4\%}$	504.9${}_{-11.8\%}{}^{}{}_{}{}^{+14.1\%}$	1.86
0	62.51${}_{-13.7\%}{}^{}{}_{}{}^{+15.8\%}$	73.64${}_{-13.4\%}{}^{}{}_{}{}^{+15.4\%}$	275.29${}_{-11.3\%}{}^{}{}_{}{}^{+13.2\%}$	1.79
1	27.84${}_{-12.9\%}{}^{}{}_{}{}^{+11.6\%}$	32.88${}_{-12.5\%}{}^{}{}_{}{}^{+13.5\%}$	127.7${}_{-10.4\%}{}^{}{}_{}{}^{+11.5\%}$	1.66
2	12.42${}_{-12.0\%}{}^{}{}_{}{}^{+13.1\%}$	14.75${}_{-11.8\%}{}^{}{}_{}{}^{+12.0\%}$	59.10${}_{-9.7\%}{}^{}{}_{}{}^{+10.2\%}$	1.56
2.4	11.65${}_{-12.7\%}{}^{}{}_{}{}^{+13.9\%}$	13.79${}_{-12.5\%}{}^{}{}_{}{}^{+13.5\%}$	53.67${}_{-10.3\%}{}^{}{}_{}{}^{+11.4\%}$	1.65
3	16.28${}_{-15.3\%}{}^{}{}_{}{}^{+16.2\%}$	19.07${}_{-14.1\%}{}^{}{}_{}{}^{+17.1\%}$	69.84${}_{-12.1\%}{}^{}{}_{}{}^{+14.6\%}$	1.90
5	81.74${}_{-15.6\%}{}^{}{}_{}{}^{+20.0\%}$	95.22${}_{-11.5\%}{}^{}{}_{}{}^{+19.7\%}$	330.61${}_{-13.6\%}{}^{}{}_{}{}^{+17.4\%}$	2.14

Equations

$$M_{\lambda }$$

$$A$$

$$M_{\lambda }$$

$$\mathcal{M}=y_t^2{\mathcal{M}}_B+y_t\lambda {\mathcal{M}}_T.$$

$$|\mathcal{M}{|}^2=y_t^4\left[{\mathcal{M}}_B{\mathcal{M}}_B^{\ast }+\frac{\lambda }{y_t}({\mathcal{M}}_B{\mathcal{M}}_T^{\ast }+{\mathcal{M}}_T{\mathcal{M}}_B^{\ast })+\frac{{\lambda }^2}{y_t^2}{\mathcal{M}}_T{\mathcal{M}}_T^{\ast }\right].$$

Full text

††preprint: CERN-TH-2019-028, MPP-2019-50, ZU-TH 09/19

1 Introduction

The Higgs potential is currently the least explored part of the Standard Model (SM), measurements of the Higgs boson self-coupling(s) may therefore offer surprises. Although the Higgs boson couplings to vector bosons and third generation fermions are increasingly well measured Khachatryan:2016vau ; Aaboud:2017vzb ; ATLAS:2018doi ; Sirunyan:2018koj ; Sirunyan:2018sgc , constraints on the trilinear coupling $\lambda$ are relatively weak due to the small Higgs boson pair production cross sections Baglio:2012np ; Frederix:2014hta . Nonetheless, measurements of double Higgs production in gluon fusion, combining various decay channels, have led to impressive experimental results already Sirunyan:2018two ; Sirunyan:2018iwt ; Aad:2019uzh ; Aaboud:2018ftw , the most stringent constraints on the trilinear coupling being $-5\leq\kappa_{\lambda}\leq 12.1$ at 95% confidence level Aad:2019uzh , based on the assumption that all other couplings have SM values. Individual limits on $\kappa_{\lambda}$ based on EFT benchmarks representing a certain combination of BSM couplings which leads to characteristic kinematic distributions Carvalho:2015ttv ; Carvalho:2016rys ; Buchalla:2018yce have also been extracted Sirunyan:2018two ; Sirunyan:2018iwt . Therefore, the determination of the trilinear coupling has entered a level of precision where the assumption that the full NLO QCD corrections do not vary much with $\kappa_{\lambda}$, which has been used in the experimental analysis so far, needs to be revised. The variations of the K-factors with $\kappa_{\lambda}$ are mild in the $m_{t}\to\infty$ limit, where NLO Grober:2015cwa ; Grober:2017gut and NNLO deFlorian:2017qfk corrections have been calculated within an effective Lagrangian framework. However, it will be shown in this paper that the NLO K-factor varies by about 35% as $\kappa_{\lambda}$ is varied between $-1$ and 5 once the full top quark mass dependence is taken into account.

The question of how large or small $\kappa_{\lambda}$ can be from a theory point of view is not easy to answer in a model independent way. Recent work based on rather general concepts like vacuum stability and perturbative unitarity suggests that $|\kappa_{\lambda}|\lesssim 4$ for a New Physics scale in the few TeV range Falkowski:2019tft ; Chang:2019vez ; DiLuzio:2017tfn ; DiVita:2017eyz . More specific models can lead to more stringent bounds, see e.g. Refs. Braathen:2019pxr ; Basler:2018dac ; Babu:2018uik ; Lewis:2017dme . Recent phenomenological studies about the precision that could be reached for the trilinear coupling at the (HL-)LHC and future hadron colliders can be found for example in Refs. Cepeda:2019klc ; Li:2019uyy ; Homiller:2018dgu ; Bizon:2018syu ; Goncalves:2018yva ; Kim:2018uty ; Adhikary:2017jtu ; Corbett:2017ieo ; Alves:2017ued ; Cao:2016zob ; Azatov:2015oxa .

Higgs boson pair production in gluon fusion in the SM has been calculated at leading order in Refs. Eboli:1987dy ; Glover:1987nx ; Plehn:1996wb . The NLO QCD corrections with full top quark mass dependence became available more recently Borowka:2016ehy ; Borowka:2016ypz ; Baglio:2018lrj . The NLO results of Refs. Borowka:2016ehy ; Borowka:2016ypz have been supplemented by soft-gluon resummation at small transverse momenta of the Higgs boson pair Ferrera:2016prr and by parton shower effects Heinrich:2017kxx ; Jones:2017giv . Before the full NLO QCD corrections became available, the $m_{t}\to\infty$ limit, sometimes also called Higgs Effective Field Theory (HEFT) approximation, has been used in several forms of approximations. In this limit, the NLO corrections were first calculated in Ref. Dawson:1998py using the so-called “Born-improved HEFT” approximation, which involves rescaling the NLO results in the $m_{t}\to\infty$ limit by a factor $B_{\rm FT}/B_{\rm HEFT}$, where $B_{\rm FT}$ denotes the LO matrix element squared in the full theory. In Ref. Maltoni:2014eza an approximation called “FTapprox”, was introduced, which contains the real radiation matrix elements with full top quark mass dependence, while the virtual part is calculated in the Born-improved HEFT approximation.

The NNLO QCD corrections in the $m_{t}\to\infty$ limit have been computed in Refs. deFlorian:2013uza ; deFlorian:2013jea ; Grigo:2014jma ; deFlorian:2016uhr . These results have been improved in various ways: they have been supplemented by an expansion in $1/m_{t}^{2}$ in Grigo:2015dia , and soft gluon resummation has been performed at NNLO+NNLL level in deFlorian:2015moa . The calculation of Ref. deFlorian:2016uhr has been combined with results including the top quark mass dependence as far as available in Ref. Grazzini:2018bsd , and the latter has been supplemented by soft gluon resummation in Ref. deFlorian:2018tah .

The scale uncertainties at NLO are still at the 10% level, while they are decreased to about 5% when including the NNLO corrections. The uncertainties due to the chosen top mass scheme have been assessed in Ref. Baglio:2018lrj , where the full NLO corrections, including the possibility to switch between pole mass and $\overline{\mathrm{MS}}$ mass, have been presented.

Analytic approximations for the top quark mass dependence of the two-loop amplitudes in the NLO calculation have been studied in Refs. Grober:2017uho ; Bonciani:2018omm ; Xu:2018eos ; Davies:2018ood and complete analytic results in the high energy limit have been presented in Ref. Davies:2018qvx . The formalism of an expansion for large top quark mass has been applied recently to calculate partial real-radiation corrections to Higgs boson pair production at NNLO in QCD Davies:2019xzc .

In this work we study the dependence of total cross sections and differential distributions on the trilinear Higgs boson coupling, assuming that the BSM-induced deviations in the other couplings are at the (sub-)percent level. The study is based on results at NLO QCD with full top quark mass dependence for Higgs boson pair production in gluon fusion described in Refs. Borowka:2016ehy ; Borowka:2016ypz . While it is unlikely that New Physics alters just the Higgs boson self-couplings but leaves the Higgs couplings to vector bosons and fermions unchanged, it can be assumed that the deviations of the measured Higgs couplings from their SM values are so small that they have escaped detection at the current level of precision, for recent overviews see e.g. Refs. Cepeda:2019klc ; Dawson:2018dcd ; Brooijmans:2018xbu ; deFlorian:2016spz .

Measuring Higgs boson pair production is a direct way to access the trilinear Higgs coupling. The trilinear and quartic couplings can also be constrained in an indirect way, through measurements of processes which are sensitive to the Higgs boson self-couplings via electroweak corrections Gorbahn:2019lwq ; Nakamura:2018bli ; Borowka:2018pxx ; Bizon:2018syu ; Kilian:2018bhs ; Vryonidou:2018eyv ; Maltoni:2018ttu ; Maltoni:2017ims ; Kribs:2017znd ; Degrassi:2017ucl ; Bizon:2016wgr ; Degrassi:2016wml ; Gorbahn:2016uoy ; McCullough:2013rea . Such processes offer important complementary information, however they are more susceptible to other BSM couplings entering the loop corrections at the same level, and therefore the limits on $\kappa_{\lambda}$ extracted this way may be more model dependent than the ones extracted from the direct production of Higgs boson pairs.

For Higgs Boson pair production, due to the destructive interference in the squared amplitude between contributions containing $\lambda$ and those without the Higgs boson self-coupling (corresponding to triangle- and box diagrams, respectively, at LO), small changes in $\lambda$ modify the interference pattern and can therefore have a substantial effect on the total cross section and differential distributions.

In order to obtain a fully-fledged NLO generator which also offers the possibility of parton showering, we have implemented our calculation in the POWHEG-BOX Nason:2004rx ; Frixione:2007vw ; Alioli:2010xd , building on the SM code presented in Ref. Heinrich:2017kxx .

The dependence of the K-factors on the value of $\lambda$ (and other BSM couplings) is stronger than the $m_{t}\to\infty$ limit may suggest, as shown in Ref. Buchalla:2018yce . This is particularly true for differential K-factors. For example, in the boosted regime, which is sometimes used by the experiments when reconstructing the $H\rightarrow b\bar{b}$ decay channel, Higgs bosons with a large-$p_{T}$ are involved. At large-$p_{T}$ the top quark loops are resolved and the $m_{t}\to\infty$ limit is invalid. The top quark mass corrections in the large $m_{\mathrm{hh}}$ or $p_{T}^{\mathrm{h}}$ regime are of the order of 20-30% or higher, and increase with larger centre-of-mass energy (e.g. $\sqrt{s}=27$ (HE-LHC) or 100 TeV (FCC-hh)), these corrections clearly exceed the scale uncertainties and therefore have to be taken into account.

The purpose of this paper is twofold: Based on our differential results, we discuss how the deviations from the SM, resulting from non-SM $\lambda$ values, can be identified based on the distributions for the Higgs boson pair invariant mass and Higgs boson transverse momentum distributions. In addition, we present the updated public code POWHEG-BOX-V2/ggHH, where the user can choose the value of the trilinear coupling as an input parameter. We also explain how variations of the top-Higgs Yukawa coupling can be studied using this code. Further, we compare the fixed order results to results obtained by matching the NLO calculation to a parton shower. In particular, we compare results from the Pythia 8.2 Sjostrand:2014zea and Herwig 7.1 Bellm:2017bvx parton showers and assess the parton-shower related uncertainties.

This paper is organised as follows. In Section 2 we briefly describe the calculation and give instructions for the usage of the program within the POWHEG-BOX. Section 3 contains the discussion of our results, focusing in the first part on variations of $\kappa_{\lambda}$ and in the second part on differences between showered results. We present our conclusions in Section 4.

2 Overview of the calculation

The calculation builds on the one presented in Ref. Heinrich:2017kxx and therefore will be described only briefly here.

The leading order amplitude in the full theory and all the amplitudes in the $m_{t}\to\infty$ limit 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. 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 Nason:2004rx ; Frixione:2007vw ; Alioli:2010xd , modified accordingly to compute the real corrections to the loop-induced Born amplitude. At run time the amplitudes were computed using Ninja Peraro:2014cba , golem95C Binoth:2008uq ; Cullen:2011kv and OneLOop vanHameren:2010cp for the evaluation of the scalar one-loop integrals. The stability of the amplitudes in the collinear limits has been improved by a better detection of instabilities in the real radiation and the use of the scalar four-point function from VBFNLO Arnold:2008rz ; Baglio:2014uba .

For the virtual corrections, containing two-loop amplitudes, we have used the results of the calculation presented in Refs. Borowka:2016ehy ; Borowka:2016ypz , which used 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 are only functions of two independent variables, the parton-level Mandelstam invariants $\hat{s}$ and $\hat{t}$. We have constructed a grid in these variables, based on 5291 pre-computed phase-space points, 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. We used the same setup for the grid as described in Ref. Heinrich:2017kxx and extended it in the following way: We can write the squared matrix element as a polynomial of degree two in $\lambda$,

[TABLE]

Therefore it is sufficient to know the amplitude at three different values of $\lambda$ in order to reconstruct the full $\lambda$-dependence. Choosing $\lambda=-1,0,1$ we obtain

[TABLE]

In practice we used the representation

[TABLE]

in order to get a more straightforward uncertainty estimate.

In fact, to any order in QCD, we can separate the matrix element into a piece that depends only on the top quark Yukawa coupling $y_{t}$ (“box diagrams”) and a piece that depends on the Higgs boson trilinear self-coupling $\lambda$ (“triangle diagrams”):

[TABLE]

The squared amplitude at each order can then be written as

[TABLE]

The above parametrisation makes it clear that the dependence of the cross section on both the Yukawa coupling and the Higgs boson self-coupling can be reconstructed from only the 3 terms present in Eq. (1). Of course this pattern changes once electroweak corrections, part of which have been calculated recently Bizon:2018syu ; Borowka:2018pxx , are included.

In order to allow for comparisons and cross checks, we implemented both the $m_{t}\to\infty$ limit as well as the amplitudes with full $m_{t}$-dependence 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 are the following:

mtdep=0**:**

computation using basic HEFT: all amplitudes are computed in the $m_{t}\to\infty$ limit.

mtdep=1**:**

computation using Born-improved HEFT. In this approximation the NLO part is computed in the $m_{t}\to\infty$ limit and reweighted pointwise in the phase-space by the ratio of the LO matrix element with full mass dependence to the LO matrix element in HEFT.

mtdep=2**:**

computation in the approximation FTapprox. In this approximation the matrix elements for the Born and the real radiation contributions are computed with full top quark mass dependence, whereas the virtual part is computed as in the Born-improved HEFT case.

mtdep=3**:**

computation with full top quark mass dependence.

Detailed instructions on how to run the code can be found in the file manual-BOX-HH.pdf in the folder ggHH/Docs of the program.

When mtdep=3 is selected, the result of the virtual matrix element is based on a grid of pre-sampled phase-space points as described above. The phase-space points present in the grid are distributed such that they optimally sample the Standard Model (SM) Born matrix element. The same set of points is used regardless of the value of $\lambda$ selected. Due to the finite number of points present in the grid, there is an associated statistical uncertainty which amounts to $0.1\%$ on the total cross section at $14\ \mathrm{TeV}$ for $\lambda=\lambda_{\mathrm{SM}}$. However, for $\lambda\neq\lambda_{\mathrm{SM}}$ the virtual matrix element can differ significantly in shape from the SM prediction, as is apparent from examining the $m_{\mathrm{hh}}$ and $p_{T}^{\mathrm{h}}$ distributions for different values of the Higgs boson self coupling. The uncertainty associated with the use of the grid is therefore larger for non-SM values of $\lambda$. The uncertainty increases as $\lambda$ is decreased below the SM value reaching $0.6\%$ on the total cross section at $14\ \mathrm{TeV}$ for $\kappa_{\lambda}=-1$. Increasing $\lambda$ above the SM value, we obtain an uncertainty of $3\%$ on the total cross section at $14\ \mathrm{TeV}$ for $\kappa_{\lambda}=3$ and $\kappa_{\lambda}=5$. Furthermore, for differential distributions the total uncertainty is not distributed uniformly in each bin but instead increases when the shape of the matrix element most differs from the SM prediction. Focusing on the invariant mass distribution, amongst the values of the Higgs boson self-coupling considered here, the largest uncertainty is obtained for the smallest values of $m_{\mathrm{hh}}$ and $\kappa_{\lambda}=3$. The uncertainty reaches $6\%$ for the lowest bin when a $40\ \mathrm{GeV}$ bin width is used.

3 Total and differential cross sections at non-SM trilinear couplings

The results were obtained using the PDF4LHC15_nlo_30_pdfas Butterworth:2015oua ; CT14 ; MMHT14 ; NNPDF parton distribution functions interfaced to our code via LHAPDF Buckley:2014ana , along with the corresponding value for $\alpha_{s}$. The masses of the Higgs boson and the top quark have been fixed, as in the virtual amplitude, to $m_{h}=125$ GeV, $m_{t}=173$ GeV and their widths have been set to zero. The top quark mass in renormalised in the on-shell scheme. Jets are clustered with the anti-$k_{T}$ algorithm Cacciari:2008gp as implemented in the fastjet package Cacciari:2005hq ; Cacciari:2011ma , with jet radius $R=0.4$ and a minimum transverse momentum $p_{T,\mathrm{min}}^{\rm{jet}}=20$ GeV. The scale uncertainties are estimated by varying the factorisation/renormalisation scales $\mu_{F},\mu_{R}$. The scale variation bands represent scale variations around the central scale $\mu_{0}=m_{\mathrm{hh}}/2$, with $\mu_{R}=\mu_{F}=c\,\mu_{0}$, where $c\in\{0.5,1,2\}$. For the case $\lambda=\lambda_{\mathrm{SM}}$ we checked that the bands obtained from these variations coincide with the bands resulting from 7-point scale variations. The PDF uncertainties have been studied in deFlorian:2016spz and found to be in general considerably smaller than the scale uncertainties.

3.1 Total cross sections at different values of the trilinear coupling

In Table 1 we list total cross sections at 13, 14 and 27 TeV for various values of the trilinear Higgs coupling $\lambda$.

Table 1 also shows that the K-factors vary substantially as functions of the trilinear coupling. This fact is illustrated in Fig. 1, showing that the K-factor takes values between 1.56 and 2.15 if the trilinear coupling is varied between $-5\leq\kappa_{\lambda}\leq 12$.

3.2 Differential cross sections

In Fig. 2 we show the $m_{\mathrm{hh}}$ distribution for various values of $\kappa_{\lambda}=\lambda_{\mathrm{BSM}}/\lambda_{\mathrm{SM}}$. The ratio plots show the ratio to the result with $\lambda_{\mathrm{SM}}$. A characteristic dip develops in the $m_{\mathrm{hh}}$ distribution around $\kappa_{\lambda}=2.4$, which is the value of maximal destructive interference between diagrams containing the trilinear coupling (triangle-type contributions) and “background” diagrams (box-type contributions). Therefore we provide results for a denser spacing of $\kappa_{\lambda}$ values around this point.

In Fig. 3 we show the transverse momentum distributions $p_{T}^{h}$ of one (any) Higgs boson for different $\kappa_{\lambda}$ values. The dip for $\kappa_{\lambda}\sim 2.4$ is still present, however much less pronounced than in the $m_{\mathrm{hh}}$ distribution.

Fig. 4 demonstrates the effect of variations of the top quark Yukawa coupling $y_{t}$ on the $m_{\mathrm{hh}}$ and $p_{T}^{\mathrm{h}}$ distributions, where $\kappa_{\lambda}$ is fixed to the SM value. Using eq. (5), it is apparent that $y_{t}$ variations can be obtained from appropriate $\kappa_{\lambda}$ variations with the same code. For example, $\sigma(y_{t}=1.2,\kappa_{\lambda}=1)=(1.2)^{4}\,\sigma(y_{t}=1,\kappa_{\lambda}=1/1.2)$.

3.3 Discussion of parton shower related uncertainties

In this section we show distributions for NLO results matched to a parton shower, focusing mostly on the transverse momentum of the Higgs boson pair. For this distribution NLO is the first non-trivial order, and therefore it is particularly sensitive to differences in the treatment of radiation by the parton shower. We compare the Pythia 8.2 Sjostrand:2014zea and Herwig 7.1 Bellm:2017bvx parton showers, applied directly to the POWHEG Les Houches events (LHE). In the Herwig case, we also compare the default shower (the angular-ordered $\tilde{q}$-shower) with the dipole shower. In addition, we assess the uncertainties stemming from the matching and show results where the Herwig shower scale parameter HardScale is varied. For all shower algorithms considered, the default tune of the corresponding version is used. Multiple-parton interactions (MPI) and hadronisation are switched off. The hdamp parameter in POWHEG is set to ${\tt hdamp}=250$ GeV.

In general, observables that are inclusive in the additional radiation, like the transverse momentum of one (any) Higgs boson, $p_{T}^{h}$, show little sensitivity to the details of the parton showering, as can be seen from Fig. 5(a), showing the fixed-order NLO prediction, as well as the Pythia 8.2 (PP8) and both Herwig 7.1 showers (angular-ordered PH7-$\tilde{q}$, and PH7-dipole). In contrast, Fig. 5(b) displays the distribution of the distance $\Delta R^{\mathrm{hh}}=\sqrt{(\eta_{1}-\eta_{2})^{2}+(\Phi_{1}-\Phi_{2})^{2}}$ between the two Higgs bosons. There, the Sudakov exponent and the parton shower effectively resum the fixed-order prediction in the region where the two Higgs bosons are close to a back-to-back configuration, and the parton shower increases the fixed-order real radiation contribution in the region $\Delta R^{\mathrm{hh}}<\pi$.

In Figs. 6(a) and 6(b), the transverse momentum $p_{T}^{\mathrm{hh}}$ of the Higgs boson pair system is shown for the fixed-order and parton-showered predictions, at $\kappa_{\lambda}=1$ and $\kappa_{\lambda}=2.4$. In all cases, the Pythia and Herwig showers agree very well in the small-$p_{T}^{\mathrm{hh}}$ range, but start to deviate already at $p_{T}^{\mathrm{hh}}\sim 100$ GeV. While both Herwig showers give very similar results and reproduce the fixed-order calculation at high-$p_{T}^{\mathrm{hh}}$, the Pythia shower produces much harder additional radiation and the ratio to the fixed-order result plateaus at $\sim 2.0$ over the remaining range. We should mention that rather large differences between Pythia 8.2 and Herwig 7.1 showers matched to POWHEG also have been found studying top quark pair production Ravasio:2018lzi . The origin of the large NLO parton shower matching uncertainties affecting certain observables in Higgs boson pair production have previously been studied in literature Jones:2017giv . For the SM result, the excess at large $p_{T}^{\mathrm{hh}}$ produced when using POWHEG with Pythia 8.2 was found to be due to additional hard sub-leading jets generated purely by the shower Bendavid:2018nar .

With the Herwig default shower, systematic uncertainties can be estimated by varying the maximal transverse momentum allowed for shower emissions, by changing the so-called hard scale $\mu_{Q}$. We apply a factor $c_{Q}=\{0.5,\,2.0\}$ on the central hard shower scale, separately for all variations of the factorisation/renormalisation scales $\mu_{R,F}$. Fig. 7 shows the $p_{T}^{\mathrm{hh}}$ and $\Delta R^{\mathrm{hh}}$ distributions as examples of the SM case, $\kappa_{\lambda}=1$, and underlines their sensitivity to changes in the shower hard scale. Quantitatively, the hard scale variations inflate the sole factorisation/renormalisation scale uncertainties by a factor of two in the regions where the Herwig 7.1 and Pythia 8.2 showers were in disagreement (see Figs. 5(b) and 6). If the envelope of all scale variations, including the hard shower scale, was to be taken as a theoretical systematic uncertainty, the resulting uncertainty would be of the order of 50% in these bins. It would be enlightening to further study parton shower (and non-perturbative) effects, in the particular context of Higgs boson pair production at NLO, as well as for loop-induced colour singlet production in general, and try to reduce discrepancies among the different algorithms.

4 Conclusions

We have presented results for Higgs boson pair production in gluon fusion at full NLO QCD for non-standard values of the trilinear Higgs boson coupling $\lambda$. We have also shown how results with a modified top quark Yukawa coupling can be produced with the same code.

We have demonstrated that the dependence of both the total and the differential K-factors on the value of $\lambda$ is stronger than the $m_{t}\to\infty$ limit may suggest. The total cross section is a quadratic polynomial in $\lambda$, with a minimum around $\kappa_{\lambda}\approx 2.4$, which is present both at LO and NLO with full top quark mass dependence, stemming from destructive interference of diagrams with and without a trilinear Higgs coupling. The $m_{\mathrm{hh}}$ distribution shows a dip around this minimum, which is to lesser extent also visible in the transverse momentum distribution of one of the Higgs bosons. We have assumed in our study that modifications of the Higgs couplings to other particles are small and can be increasingly well constrained by other processes. Nonetheless, it should be kept in mind that a dip in the $m_{\mathrm{hh}}$ distribution could also originate from other effective couplings, for example an effective $t\bar{t}HH$ coupling, while $\kappa_{\lambda}=1$ Buchalla:2018yce .

We have also combined our NLO QCD results with the Pythia 8.2 and Herwig 7.1 parton showers. In the Herwig 7.1 case we employed both the default shower (the angular-ordered $\tilde{q}$-shower) and the dipole shower. We observed that for distributions particularly sensitive to the additional radiation, the parton showers exhibit a somewhat different behaviour. While both Herwig 7.1 showers generate comparable results and perform as expected in the NLO regime, the Pythia 8.2 shower produces harder radiation, for example in the tail of the $p_{T}^{\mathrm{hh}}$ distribution. Varying the shower hard scale in Herwig 7.1 on top of $\mu_{R},\mu_{F}$ variations leads to uncertainty bands which approximately cover these differences. However, the parton shower uncertainties can then become sizeable and even surpass the fixed-order scale uncertainties.

The POWHEG version of the code for Higgs boson pair production including the possibility to vary the trilinear coupling and the top quark Yukawa coupling is publicly available in the POWHEG-BOX-V2 package at the website http://powhegbox.mib.infn.it, in the User-Processes-V2/ggHH/ directory.

Acknowledgements

We would like to thank Gerhard Buchalla, Alejandro Celis, Matteo Capozi, Stephan Jahn and Emmanuele Re for helpful discussions. This research was supported in part by the COST Action CA16201 ‘Particleface’ of the European Union, and by the Swiss National Science Foundation (SNF) under grant number 200020-175595. We gratefully acknowledge resources provided by the Max Planck Computing and Data Facility (MPCDF).