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

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.


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 → bbbb has so far led to the most stringent limit on the di-Higgs production cross section of σ/σ SM ≤ 29 in ATLAS [1], while the CMS collaboration achieved the most restrictive upper bound of σ/σ SM ≤ 28 in the bbτ τ decay channel [2]. A previous combination of various decay channels measured in the ATLAS detector led to σ/σ SM ≤ 70 [3]. The CMS collaboration also produced new limits for resonant and non-resonant Higgs boson pair production in the bbV V channel [4].
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. [5][6][7]. 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 → ∞ limit, leading to point-like effective couplings of gluons to Higgs bosons. In Ref. [8] NLO corrections were calculated in the so-called "Born-improved HEFT" approximation, where the basic HEFT result is rescaled by a factor B F T /B HEF T , B F T denoting the leading order matrix element squared in the full theory.
In Refs. [9,10], an approximation called "FT approx " 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 F T /B HEF T .
The NNLO QCD corrections in the heavy top limit have been computed in Refs. [12,[15][16][17], and they have been supplemented by an expansion in 1/m 2 t in Ref. [13] and by resummation, at NLO+NNLL in Ref. [18] and at NNLO+NNLL in Ref. [19], 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 [20,21] and compared to previous approximations for various observables [22]. The full NLO calculation was supplemented by NLL resummation in Ref. [23].
Numerous phenomenological studies of Higgs boson pair production have been performed both within and beyond the SM . Further, it also has been suggested recently to obtain constraints on the Higgs boson self-coupling from electroweak corrections to single Higgs boson production [54][55][56].
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 [20,22,57].
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. [58]. The "FT approx " [10] calculation includes the matching of di-Higgs production to a parton shower [9] 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 [59,60] and the MadGraph5_aMC@NLO framework [61,62]. This allows us to compare the POWHEG [59] and MC@NLO [63] matching schemes while using the same Pythia 8 [64,65] shower in both cases. We also investigate the PDF and scale uncertainties and calculate observables like the di-Higgs p T spectrum, p hh T , where fixed order NLO calculations cannot give a satisfactory description at low p hh T . 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 [17] with the showered results.

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 FT approx within MG5_aMC@NLO already have been published elsewhere [9,10].
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 FT approx (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".

Virtual two-loop amplitudes
For the virtual two-loop amplitudes, we have used the results of the calculation presented in Refs. [20,22], which is based on an extension of the program GoSam [66] to two loops [67], using also Reduze 2 [68] and SecDec 3 [69].
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ŝ andt. 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ŝ andt to new variables where f can, in principle, be any strictly increasing function. Setting f (β) 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 θ ) 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 θ , 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 [70] as implemented in the python SciPy package [71]. Applying this procedure reduces the size of interpolation artefacts, which we obtain due to the numerical uncertainty of our two-loop results.
We should point out that the grid is constructed from a sample of phase space points which is based on runs at √ 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ŝ due to a lack of statistics in this region upon construction.
In our original calculation [20,22], we used Catani-Seymour dipole subtraction [72] for the real radiation. The finite combination of the renormalized virtual amplitude V b with the Catani-Seymour I-operator can be straightforwardly converted into the quantity V fin of Refs. [60,73], defined by For this process the colour-correlated Born squared amplitudes B 12 and B 21 are equal and we have a = 2C A and c 12 = c 21 = β 0 /2. In the POWHEG-BOX framework the arbitrary scale Q is chosen as µ r . Specifically, we obtain (2.5) The grid evaluates V fin at the scale µ 0 = √ŝ /2 and the results for an arbitrary scale can be obtained from the relation The Born amplitude B and the colour-correlated Born amplitudes B ij in (2.3) are evaluated in D = (4 − 2 ) dimensions using conventional dimensional regularization (CDR). As all formulas for the soft contributions and the collinear remnants used in the POWHEG-BOX are computed in the MS scheme, using CDR, constructing V fin according to (2.5) ensures the treatment is consistent with that implemented in the POWHEG-BOX.

Real radiation and parton shower matching
The real radiation matrix elements in the full SM were implemented using the interface [74] between GoSam [66,75] and the POWHEG-BOX [59,60], 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 [66], that uses QGRAF [76], FORM [77] and Spinney [78] for the generation of the Feynman diagrams, and offers a choice from Samurai [79,80], golem95C [81][82][83] and Ninja [84,85] for the reduction. At run time the amplitudes were computed using Ninja [84,85] and OneLOop [86] 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 hh T > 10 −3 GeV has been introduced. We carefully checked that the total cross section does not change significantly when varying the cut value.
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.

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 [87]. 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 [88][89][90][91] parton distribution functions interfaced to our codes via LHAPDF [92], along with the corresponding value for α 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 T algorithm [93] as implemented in the Fastjet package [94,95], with jet radius R = 0.4 and a transverse momentum greater than p min T = 20 GeV. The theoretical scale uncertainty is estimated by varying the factorization scale µ F and the renormalization scale µ R . The scale variation bands are obtained by computing the envelopes of a 7-point scale variation around the central scale 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.

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. [20,22]. 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. 1 we show a comparison between the NLO predictions obtained with the POWHEG-BOX generator and the original ones from Ref. [22].  (b) m hh in the full SM  first hard emission is weighted with the Sudakov factor according to the POWHEG method. Despite the LHE level predictions being unphysical, this 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. 2, 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 hh T distribution in the limit p hh T → 0, one instead expects to observe Sudakov suppression in the soft region. This can be seen in Fig. 3, 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 hh T , 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 [96][97][98], and we refer the interested reader to Refs. [96,97] 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 sing and a regular part R reg , as follows: where the transition function F is chosen to be In Fig. 4 we compare the default POWHEG choice setting, h = hdamp = ∞, with predictions in which we set hdamp = 250 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 GeV the LHE curve with hdamp = 250 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. 5 we show a similar comparison for m hh (left) and the transverse momentum of a (randomly chosen) Higgs boson p h T (right), with full top quark mass dependence. The m 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 h T , 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 hdamp = 250 are in general closer to the NLO ones over the whole kinematical range of p h T . We stress however that, contrary to p hh T , where the differences between the predictions for hdamp = ∞ and the one for hdamp = 250 GeV reach  (b) p h T in the full SM 80% above 500 GeV, for p h T 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.

Discussion of NNLO effects
We mentioned in the previous section that the enhancement in the tail of p hh T 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 [96], 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. [17] 1 . 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. 6 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 hdamp = ∞, on the right we set hdamp = 250 GeV. In the former plot we observe a good agreement of the LHE-level curve with hdamp = ∞ with the NNLO predictions in the transverse momentum range between 200 GeV and 400 GeV. While the LHE-level result flattens out around 250 GeV, the NNLO result decreases slightly for larger p hh T . 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 7 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 hh T distribution, such that the NNLO curve lies between the NLO+PS and the NLO fixed order curve for hdamp = ∞. On the other hand, for hdamp = 250 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.

NLO plus parton shower matched results
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 h 1 T ) and "softer" (p h 2 T ) 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 h 1 T and p h 2 T 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 → 2 contributions. To illustrate this fact, we consider the total cross section as a function of ∆, with the kinematic requirements p h 1 T ≥ ∆, p h 2 T ≥ 0. The cross section shows an unphysical behaviour as ∆ → 0, see Fig. 8: the total cross section as a function of ∆ peaks around ∆ = 14 GeV and then decreases for smaller values of ∆, even though the available phase space for p h 1 T 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 → 2 calculation at NLO) should be avoided. For a more detailed discussion of this point we refer to Refs. [99][100][101]. Here we only note that this is the reason why, with "symmetric" cuts p h 1 T,min = p h 2 T,min = 0 and fine binning, the first bin(s) of the p h 1 T 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. T distribution diverges at fixed order for p hh T → 0, while the showered result is able to provide reliable predictions in the low p hh T 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 bindetermined by the 2 → 3 kinematics, which is described only at leading order accuracy by our calculation.
In Fig. 12 we show the difference in azimuthal angle, ∆Φ hh , and the radial separation,    ∆R hh = (η 1 − η 2 ) 2 + (Φ 1 − Φ 2 ) 2 , of the two Higgs bosons. We see that the unphysical behaviour for ∆Φ hh → 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 ∆R hh distribution, we observe that the shower populates the region ∆R hh < π, which at fixed order is given by the 2 → 3 component only.
In Fig. 13a we compare the full results with calculations where the underlying matrix elements are based on two approximations, either FT approx 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. 13b. 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.
Because of the fact that for the p hh T distribution, the tail is predicted at the first nontrivial order, the effect of the shower on this distribution is rather large, exceeding a factor of two beyond p hh T ∼ 300 GeV, as shown in Fig. 14. However, as can also be seen from  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 hh T > 500 GeV. Fig. 14 also shows that FT approx does a good job in this case, as the tail of the p hh T distribution is determined by the real radiation. The FT approx curve still lies above the full result though, because of the differences in the virtual part that enters theB function in POWHEG, which determines the overall normalisation for the shower.
Finally, in Fig. 15 we compare results obtained with the Pythia 6 shower to the Pythia 8 results (in the basic HEFT approximation) and find that the differences are small.

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ŝ andt 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.  T and p h 2 T 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= ∞, which is the default, we can say that in general the two predictions agree well within the scale and statistical uncertainties. The small ∆R hh region 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 MG5_aMC@NLO the shower starting scale is picked with some probability distribution to be in the interval shower_scale_factor × [0.1 √ŝ , √ŝ ], therefore to perform the scale variation we set the shower_scale_factor in the run card to 0.5, 1 and 2.
The m hh distribution can be considered as a control plot to demonstrate that, as expected, this has no effect on the m hh distribution. In contrast, in the p hh T distribution, the

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ŝ andt 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 h T , 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 hh T , ∆Φ hh and ∆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 hh T distribution, around p hh T ∼ 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 hh T 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 hh T ∼ 400 GeV, worsening towards higher p hh T values. As expected, the FT approx 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 hh T . 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.
acknowledge support and resources provided by the Max Planck Computing and Data Facility (MPCDF).