Parton Shower and NLO-Matching uncertainties in Higgs Boson Pair Production

We perform a detailed study of NLO parton shower matching uncertainties in Higgs boson pair production through gluon fusion at the LHC based on a generic and process independent implementation of NLO subtraction and parton shower matching schemes for loop-induced processes in the Sherpa event generator. We take into account the full top-quark mass dependence in the two-loop virtual corrections and compare the results to an effective theory approximation. In the full calculation, our findings suggest large parton shower matching uncertainties that are absent in the effective theory approximation. We observe large uncertainties even in regions of phase space where fixed-order calculations are theoretically well motivated and parton shower effects expected to be small. We compare our results to NLO matched parton shower simulations and analytic resummation results that are available in the literature.


Introduction
In the Standard Model (SM) the production of a pair of Higgs bosons at a hadron collider proceeds dominantly through the annihilation of gluons. As there is no direct coupling between the Higgs boson and gluons this process is mediated by intermediate massive quark loops. The top-quark loop contributions dominate by far, due to the large Yukawa coupling. Bottomquark loops only contribute at the 1 % level to the total cross section at leading order (LO) and can thus safely be neglected in most situations. The scattering amplitudes relevant for the calculation of the top-quark contributions up to next-to-next-to leading order (NNLO) are known in the approximation of an infinitely heavy top-quark (commonly referred to as the HEFT approximation) [11,14,13,26,25,12]. However, the validity of this approximation is questionable for Higgs boson pair production due to the large momentum transfer required to produce the Higgs bosons. Techniques to systematically improve upon it were extensively studied in [19,33,24,26,25,32]. The full result, which is exact in the mass of the top-quark, was known only to leading order (LO) [16,23,38] until recently. This is due to the complexity of computing the next-to-leading order (NLO) virtual corrections which feature two-loop, four-point integrals with both massive internal propagators and massive external lines. They have recently been calculated through the numerical evaluation of all relevant two-loop integrals as part of the full NLO calculation in [6,5].
At small transverse momenta p HH ⊥ of the Higgs boson pair, the accuracy of any fixed-order calculation is spoiled by the presence of large logarithms of the form α n s log(p HH ⊥ /m HH ) m .
They can be resummed to all orders using analytical resummation techniques which have been applied to Higgs boson pair production in [18]. Alternatively, parton shower simulations can be employed. In addition to providing a reliable transverse momentum spectrum at small p HH ⊥ , they also provide results that are fully differential in the kinematics of any soft and collinear QCD radiation. Standard techniques exist for the consistent matching of NLO fixed order calculations to parton shower simulations [20,34]. They were recently applied to Higgs boson pair production in reference [27], where the MC@NLO and POWHEG matching techniques were used to combine the fixed-order NLO calculation with the Pythia parton shower [41,40]. The results of [27] suggest that the parton shower matching can have sizeable effects not only in the region of small p HH ⊥ , but also in the region of large p HH ⊥ , where one would expect the fixed-order calculation to be reliable and the approximations inherent to parton shower simulations to break down. These effects even exceed the scale uncertainties of the fixed-order calculation.
In this publication we aim to critically assess the origin and size of the aforementioned effects and associated uncertainties. For this purpose we implemented a fully generic and process independent NLO subtraction along with the corresponding parton shower matching techniques for loop-induced processes in the Monte Carlo event generator Sherpa [21]. This allows us to perform our studies using the two different showers that are implemented in Sherpa [39,42] within the same parton shower matching framework.
This publication is structured as follows. In Section 1 we describe in detail the setup of our calculation along with a brief review of the MC@NLO matching technique and the parton showers we used for our studies. We present the results of our simulations in Section 2, focusing on the origin and size of uncertainties that are inherent to the matching technique applied. We also point out crucial differences that arise when going from the HEFT approximation to the full calculation. Our conclusions are presented in Section 3.

Fixed-Order NLO Calculation
For the virtual two-loop amplitude we utilize the result of reference [6,5], retaining the full finite top-quark mass effects. This amplitude was obtained by numerically evaluating all relevant 2loop 4-point Feynman diagrams with up to 4 scales. We adopt the input parameters of reference [5], with G F = 1.1663787 × 10 −5 GeV −2 , the mass of the top-quark set to m t = 173 GeV, the Higgs boson mass set to m H = 125 GeV, and their widths neglected. We also adopt the choice of reference [5] for the factorization and renormalization scales µ F = µ R = m HH /2. Perturbative uncertainties in the fixed-order part of the calculation are estimated by independently varying these scales through factors of 0.5 and 2. All studies are performed with hadronic center-of-mass energy √ s = 14 TeV. The NLO virtual amplitude is provided in the literature in the form of an interpolation grid in two Mandelstam variables, based on a fixed number of pre-computed phase-space points [27]. We extract the finite part of the UV renormalized virtual amplitude in the Conventional Dimensional Regularization (CDR) scheme with residual IR singularities subtracted according to the Catani and Seymour scheme [10], as required by the Sherpa event generator, using relations (2.5) and (2.6) of reference [27].
The leading order one-loop squared amplitudes for the Born process and real emission contributions are provided by OpenLoops [9]. For the evaluation of tensor and scalar one-loop integrals, we employ the Collier library [15], CutTools [37], and OneLOop [44,43].
For the regularization and numerical cancellation of infrared divergences in the real-emission part of the calculation we employ the dipole subtraction scheme of Catani and Seymour [10]. We have re-implemented this scheme within Sherpa in a fully generic and process-independent way for loop-induced processes. This implementation is qualitatively equivalent to the implementation in one of Sherpa's internal matrix element generators Amegic++ [31,22], apart from the fact that color-and spin-correlated amplitudes are to be provided externally through generic interfaces. Through a dedicated interface to OpenLoops and the aforementioned tools, NLO calculations for loop-induced SM processes are thus fully automated (given the availability of the virtual two-loop corrections) in Sherpa and will become available in a public version of the code. We have validated our implementation in Sherpa by comparing our results for the total cross section, for the differential Higgs boson pair invariant mass distributions, and for the differential single Higgs boson transverse momentum distributions to those published in reference [5].

Parton Showers
We consider two parton showers for matching to the fixed-order NLO calculation. Both algorithms are dipole-type showers in which QCD radiation is generated coherently off color dipoles spanned by pairs of pre-existing partons. Both showers are publicly available as part of the Sherpa event generator package.
Due to the dipole character of the parton showers, their splitting kernels can be used for the purpose of fixed-order NLO subtraction, thus simplifying the implementation of parton shower matching. The CS shower [39] directly uses the splitting kernels of the original Catani-Seymour subtraction scheme for parton evolution. The Dire shower [42] uses splitting kernels that are modified in such a way as to reproduce the collinear anomalous dimensions of the DGLAP equations. For NLO matching to the Dire shower, we use a modified version of the original Catani-Seymour subtraction scheme that reflects these changes in the splitting kernels [29]. The kernels of both showers approximate real emission amplitudes arbitrarily well in the limit of soft and collinear momenta. Away from the soft and collinear regions, however, they differ.
A further crucial difference between the two algorithms is the choice of evolution variable, which we generically denote by t in the following. The choice of evolution variable together with the shower starting scale µ PS dictates how much of the phase space away from the soft and collinear regions is available to the parton shower since the starting scale implements the following phase space constraint: For the discussion of the evolution variables we focus on the first (hardest) emission in the production of a color-neutral final state of invariant mass Q 2 = m 2 HH . We illustrate the kinematics of the first emission, producing a final state parton with momentum p j from the collision of two incoming massless partons with momenta p a and p b , in Figure 1. It is useful to consider the variables v and w, which are closely related to the standard Mandelstam variablest = (p a − p j ) 2 , Due to momentum conservationŝ +t +û = Q 2 which implies, In terms of v and w, the evolution variable in Dire is given by The functional form (4) implies that t Dire < Q 2 4 due to equation (3). It follows that for a parton shower starting scale of µ 2 PS = Q 2 4 , Dire behaves like a "power shower" in the sense that it populates the full phase space since t Dire < Q 2 4 and thus (1) is trivially fulfilled. In the CS shower, the evolution variable is given by This implies that, for a given kinematic configuration, t CSS is typically larger than t Dire , such that for a given value of µ PS the emission phase space of the CS shower is more restricted than that of Dire. It is worth noting that µ 2 PS = Q 2 4 in particular does not correspond to a "power shower" when employing the CS shower. This choice in fact severely constrains the emission phase space, since v + w can get close to 1 and thereby give rise to large values of t CSS .

NLO Parton Shower Matching
In the following we will focus on the MC@NLO matching prescription [20] using the notation of [30] with no distinction between fixed-order NLO subtraction terms D (S) and parton shower matching terms D (A) since we use the parton shower splitting kernels both for parton evolution and for infrared subtraction and keep all phase space constraints explicit. We thus denote the sum of subtraction terms as a function of the real emission phase space by D(φ R ), where the real emission phase space φ R = φ B × φ 1 can be decomposed into the Born phase space φ B and an extra one-particle emission phase space φ 1 . We then define the fixed-order differential seed cross sectionsB(φ B ) and H(φ R ) in terms of the leading order (Born) term B(φ B ), the UV-subtracted virtual corrections V (φ B ), and the real-emission contributions R(φ R ) bȳ where t(φ R ) is the map from a kinematic real emission configuration to the parton shower evolution variable t. The Heaviside function Θ(µ 2 PS − t(φ R )) in (7) implements the constraint (1). For notational convenience, we will omit the explicit φ R -dependence and write t(φ R ) = t in the following. In terms of the quantities introduced above, the fixed-order total NLO cross section is given by In MC@NLO, we generate events according to where t 0 is the infrared cutoff scale of the parton shower and the modified Sudakov form factor ∆(t 0 , t 1 ), which gives the probability for no emission to occur between scales t 0 and t 1 for the first parton shower step, is given by The first line of (9) corresponds to events that have Born kinematics at the level of the fixed-order seed event with weightB (S-events). They either don't undergo any emission above the infrared parton shower cutoff scale t 0 (first term in the square bracket) or they undergo their hardest emission at some scale t between µ 2 PS and t 0 (second term in the square bracket). The second line of (9) corresponds to events with real-emission kinematics at the level of the fixed-order seed event and weight H (H-events). All events are treated further by the parton shower precisely as in the leading order case, apart from the S-events that haven't undergone any emission, which are kept as they are.
Since the square bracket in (9) integrates to 1, the total cross section and any observable that is insensitive to QCD radiation is unaltered in MC@NLO compared to the fixed-order NLO result. In fact, it can be shown that a MC@NLO event sample will reproduce the fixed-order NLO result event to order α S relative to the Born for any infrared safe observable [20]. The parametric NLO accuracy is therefore not spoiled by the parton shower matching.

Parton Shower Matching Uncertainties
As stated in the previous section, NLO parton shower matching according to the MC@NLO method preserves the parametric accuracy of the fixed-order NLO calculation. Deviations from fixed-order results can numerically be nonetheless significant [30]. Such differences reflect genuine parton shower matching uncertainties, they can be particularly prominent for observables that are sensitive to real emission configurations and thereby to the interplay between parton shower emissions and fixed-order real emission configurations. We will therefore focus on the p HH ⊥ distribution in the following section, comparing MC@NLO matched parton shower simulations to fixed-order results with both the Dire and the CS shower.
In order to formally compare the MC@NLO result to a fixed-order prediction for this spectrum, we first consider a generic observable O that is insensitive to kinematic Born configurations. For such an observable we need to take into account H-events and parton shower emissions off Sevents. At order α S relative to the Born we have where the first integral corresponds to S-events in which the parton shower has generated a nonvanishing value of O and the second integral corresponds to H-events, where a non-vanishing value of O is implied by the real-emission kinematics of the fixed-order seed event. In the tail of the distribution where we can neglect the Sudakov suppression and set ∆ = 1, we obtain after plugging in the definition of H: To order α s we haveB = B and the first integral cancels as required by the matching conditions, thus restoring the fixed-order result. This explicitly demonstrates how variations in the parton shower contributions induced by S-events are subtracted by the MC@NLO subtraction terms D in the definition of H. Numerically, however, this cancellation can be severely spoiled, potentially leading to large deviations from the fixed-order result. For the deviations to be significant the term on the first line of equation (13) must be similar in size to the fixed-order term on the second line of (13). One can therefore expect large deviations from the fixed-order calculation only if both of the following conditions are met: 1) The factorB −B dressed with the parton shower splitting kernels ( D B in (13)) is comparable in magnitude relative to the real-emission matrix elements in R. This depends on the size of the NLO corrections that enterB and on the splitting kernels in the phase space region of hard emissions.
2) The phase space of interest is accessible to the parton shower so that the first integral in (13) has support in that region. This depends on the choice of µ PS and on the shower (through the functional form of t(φ R )).
The formally sub-leading contributions originating from the parton shower matching in the first integral in (13) are, to a large extent, controlled by the choice of µ PS . To access the matching uncertainties we will therefore vary this parameter by factors 2 and 0.5. With two different parton showers at our disposal we have an additional handle on these uncertainties through the functional form of t(φ R ). The nominal choice for µ PS in the CS shower will be µ PS = m HH /2, in line with µ R and µ F . As outlined above, such a choice would open up the entire emission phase space in case of the Dire shower. Our nominal choice for the Dire shower will therefore be µ PS = m HH /4, which allows us to perform both up and downwards variations.
Based on the argument presented above, one might expect to see large parton shower contributions in the high-p T tails of other processes with large K-factors. However, it is important to note that for such effects to be visible theB − B factor must remain large, relative to the real-emission matrix element, also when multiplied by the parton shower splitting kernels. In single Higgs boson production through gluon fusion, for example, one might anticipate a large shower contribution in the tail of the Higgs transverse momentum spectrum due to the large NLO K-factor. However, in this case the parton shower splitting kernels underestimate the realemission matrix elements significantly, such that the parton shower contributions in the tail are very small in an MC@NLO-matched calculation [1,30]. For the parton showers considered in this work, this holds even when taking into account the full top mass dependence in the realemission matrix elements, which reduce the size of R by more than an order of magnitude in the tail of the transverse momentum distribution.

Leading Order Results
We start our discussion with predictions obtained in the most simple setup, using leading order matrix elements for inclusive Higgs boson pair production supplemented by a parton shower. This type of simulation will be referred to as "LO+PS" in what follows. Since the transverse momentum of the Higgs boson pair is zero at leading order, any non-zero value of this observable is entirely generated by the parton shower. As a reference, we use a fixed-order prediction obtained by simulating the process p p → H H j with leading order matrix elements. Figure 2 and 3 show the result of our comparison both in the HEFT approximation and in the full theory, respectively. Comparing the full SM and the HEFT approximation, we observe qualitatively different parton shower effects. In the HEFT approximation, both parton showers significantly underestimate the fixed order spectrum in the tail of the distribution. Even if the full phase space is made available to the showers they do not reproduce the slowly falling transverse momentum spectrum predicted by the fixed-order HEFT matrix elements. To show this for the CS shower we display also LO+PS results obtained by setting the parton shower starting scale to the hadronic center-of-mass energy µ PS = √ s. In the case of the Dire shower, the full phase space is already available for µ PS = m HH /2, which corresponds to the upper edge of the uncertainty band.
In the full SM, by contrast, for large enough values of the parton shower starting scale µ PS both parton showers overestimate the fixed-order prediction. For the CS shower, this effect is restricted to smaller transverse momenta, due to the choice of evolution variable. If we lift any phase space restriction in the CS shower, by setting µ PS = √ s, we observe that in the tail of the distribution the shower overestimates the fixed-order predictions by more than an order of magnitude. The upper edge of the Dire shower uncertainty band also overestimates the fixed-order prediction, although this feature is not as pronounced as for the CS shower.
It is therefore evident that the Born matrix elements dressed with splitting kernels can strongly overestimate the real emission matrix elements in the full SM and strongly underestimate the real emission matrix elements in the HEFT. The former effect is, however, to a certain extent      limited by the phase space constraint implemented through the parton shower starting scales.
Naively, the large differences between the HEFT and the full SM simulations may seem surprising. However, the high-energy behaviour of the HEFT real emission amplitudes is unphysical because the momentum transfers vastly exceed the top-quark masses which have been integrated out in the HEFT approximation. As a result, the spectrum calculated at fixed-order falls off extremely slowly in the HEFT and the parton shower kernels thus tend to underestimate the spectrum in the tail. A similar slow fall off can been observed in the HEFT approximation for the Higgs boson transverse momentum in single Higgs boson production [4,17,7,36,8].

NLO Results
We start the discussion of NLO-matched parton shower simulations with the results of a HEFT treatment, shown in Figure 4. As discussed in Section 2.1 (and shown in Figure 2) the combination of Born matrix elements and parton shower splitting kernels strongly undershoot the full real-emission matrix elements when employing the HEFT approximation. We therefore expect the tail of the distributions to converge to the fixed-order result both for the CS shower and the Dire shower. As shown in Figure 4, this is indeed the case. In addition to a more precise description of the tail (compared to the LO+PS type simulations) we observe a reduction of the parton shower starting scale uncertainties. The individual variations of H and S-event contributions are of order one for some values of p HH ⊥ but cancel to a large extent in the sum. Moving on to the discussion of results in the full SM, we remind the reader of our findings in the corresponding LO+PS type simulations. In the full SM, the parton shower splitting kernels in combination with Born matrix elements overestimate the real-emission matrix elements (see Figure 3). The parton shower effects in the tail of the p HH ⊥ distribution are therefore large. As shown in Figure 5, this also holds at NLO. The parton shower matched results converge to the fixed order result in the tail for nominal choices for µ PS . Upward variation of µ PS , however, (indicated by the upper edge of the blue uncertainty bands) lead to parton shower effects of up to +100 % even in the tail of the distribution. As shown in the lower panels of Figure 5, the excesses in the tail are indeed generated by parton shower emissions off S-events. The extent of these effects is limited by the phase space available to the parton shower, as determined by the choice of µ PS and the functional form of the evolution variable. We observe that results generated using the Dire shower, particularly for larger values of µ PS , have a different shape than those generated with the CS shower.
For the MC@NLO algorithm, by construction, the large parton shower effects in the tail should be cancelled to first order in α S . As outlined in Section 1.4, any mis-cancellation is due to a numerically large discrepancy between B andB. We demonstrate this explicitly in Figure 6, where we show modified Dire MC@NLO with B substituted forB, leading to a complete cancellation of the first integral in (13). This procedure eliminates large parts of the excess in the tail independently of µ PS , as anticipated. Variations in S-and H-event contributions remain large, as shown in the lower panels of Figure 6, but they cancel in the sum. The procedure of replacingB with B would, of course, spoil the NLO accuracy of any inclusive observable but allows us here to demonstrate the origin of the discrepancy between the showered and fixed-order results in the tail of the p HH ⊥ distribution.
In the HEFT approximation one may naively expect effects of similar size in the tail of the distributions. As demonstrated using a LO+PS simulation and as shown in Figure 2, however, the fixed-order real emission contributions completely dominate in this region. The bulk of the contributions in the tail are hence generated by the second integral of (13). As a result, the relative impact of parton shower effects in the tail remains small as shown in Figure 4.

Comparison to the Literature
In Figure 7 we compare our results for the p HH ⊥ spectrum to the NLO parton shower matched results presented in reference [27]. These results were obtained with the Pythia 8 shower [41,40] interfaced to MadGraph5 aMC@NLO [3,28] and POWHEG BOX [2] for matching according to the MC@NLO method and the POWHEG method [34], respectively. In MadGraph5 aMC@NLO the nominal value of µ PS is set randomly in the interval [0.1H T /2, H T /2], where H T is the sum of the transverse masses of the Higgs bosons. For the simulations based on MadGraph5 aMC@NLO and Pythia we show uncertainty bands that were obtained by varying the nominal parton shower starting scale by factors of 2 and 0.5. The POWHEG matching prescription can be recovered within the MC@NLO framework by setting the parton shower starting scale to the collider energy µ PS = √ s and by setting D = R [30]. Therefore, lacking a natural equivalent to µ PS in the POWHEG framework, we compare only to nominal POWHEG predictions produced with the hdamp parameter set to 250 GeV, as described in [27].
Focusing on the region p HH ⊥ > 100 GeV, we note that all MC@NLO predictions considered here are generally compatible within the uncertainty bands. However, the agreement between the nominal results of our simulations and the fixed-order result is much better in the tail. The POWHEG results exhibit a very large excess in the tail that is not covered by the uncertainty bands of our Sherpa predictions. Similar discrepancies between MC@NLO and POWHEG have been observed in the context of other processes [1,30] and can be attributed to the large phase space available to the parton shower as a result of setting µ PS = √ s and the numerically large discrepancy betweenB and B [35]. As described in Section 1.2, the former can be achieved in Dire by setting µ PS = m HH /2, which is represented by the upper edge of the uncertainty band around the Dire prediction. We note that the shape of this curve is in fact most similar to the POWHEG prediction in the tail of the distribution.
Comparing the different uncertainty bands themselves, we observe large differences. The shape of uncertainty bands obtained with MadGraph5 aMC@NLO and with the CS shower are somewhat similar with a peak around 300 GeV but the size of the uncertainty band around the MadGraph5 aMC@NLO result is much larger throughout. The uncertainties on the Dire prediction describe a more evenly shaped band.
Differences in the region of small transverse momenta are not fully covered by the µ PS -variation bands. Although we expect these variations to be indicative of NLO-matching uncertainties, we do not expect them to cover all parton shower uncertainties. These include, but are not limited to, ambiguities in the choice of a kinematic recoil scheme and ambiguities in the choice of the renormalization scales for the strong coupling in the splitting kernels.
In Figure 8 we show a comparison to the calculation of reference [18] which employed analytic next-to-leading-log (NLL) resummation techniques instead of parton showers. We observe good agreement within the uncertainties except near the peak region and the region around p HH 100 GeV where we find discrepancies of about 5 % that are not fully covered by our uncertainty bands. Taking into account the resummation scale uncertainties on the analytic results (not shown in Figure 8), which are of the order of 3 % in the peak region and 10 % near p HH ⊥ = 100 GeV [18], we consider the observed agreement satisfactory.

Other Observables
Having discussed the Higgs boson pair transverse momentum at length, we briefly discuss parton shower effects on a number of other observables.
Lorentz invariant observables that depend only on the momenta of the Higgs bosons are not affected by the parton shower. The kinematics of the Higgs boson pair is only altered by emissions off dipoles spanned by two initial state partons. The recoil generated by such an emission affects the Higgs boson pair only through a Lorentz boost. Lorentz-invariant quantities like the Higgs boson pair invariant mass are therefore not affected by the parton shower. Our simulations were checked by inspecting the Higgs boson pair invariant mass distributions, we observe agreement within the statistical uncertainties of well below one percent. Figure 9 shows the leading jet transverse momentum p j ⊥ . As opposed to the Higgs boson pair, the parton emitted in the hardest emission off an S-event and the final state parton in an H event is affected strongly by secondary emissions because the kinematics are altered by final-final and final-initial dipoles. Such emissions decorrelate the leading jet and the Higgs boson pair momenta. Parton shower effects are therefore qualitatively different. The effect of  the parton shower are generally moderate and the spectrum remains compatible with the fixedorder prediction within the scale uncertainties. It is worth noting that even the full µ PS -variation bands remain within the uncertainty bands of the fixed-order calculation. This picture drastically changes when considering observables that are more sensitive to high multiplicity final states. As an example we consider the differential H T distribution, defined as the scalar sum of jet transverse momenta: where the index i labels all jets in the respective event. In a parton shower event, the total energy is typically distributed among a larger number of jets than in a fixed-order calculation. Their scalar contributions to H T are added, giving rise to larger values of H T than for the fixed-order NLO calculation. We show this effect in Figure 9. Comparing the Dire and CS shower predictions, we note that the uncertainty bands overlap but that the shower starting scale variations are much larger for the CS shower.
In Figure 10 we show the azimuthal separation between the Higgs bosons ∆φ HH . At leading order, the momenta of the Higgs bosons are perfectly correlated due to momentum conservation. Only in events with additional radiation can one observe a non-trivial distribution of the azimuthal separation between the Higgs bosons. As shown in Figure 10, parton shower corrections to the fixed-order result are mostly covered by the fixed-order uncertainties except in the region of ∆φ HH = π which corresponds to back-to-back configurations and which is sensitive to soft QCD emissions.
Also shown in Figure 10 is the transverse momentum of a randomly chosen Higgs boson. The effect of a parton shower emission on the transverse momentum of a given Higgs boson is random, either decreasing or increasing its value. If the distribution was completely flat, any parton shower effects would therefore average out. Since the distribution is falling, the parton shower effects of increasing the transverse momenta of low-p ⊥ Higgs bosons is not counter-balanced by the effect of decreasing the transverse momenta of high-p ⊥ Higgs bosons, thus inducing a slope relative to the fixed-order result. This effect is small but clearly visible in Figure 10.

Conclusions
We have presented a study of NLO parton shower matching uncertainties in Higgs boson pair production through gluon fusion at the LHC. We assessed these uncertainties by matching the fixedorder NLO calculation to two dipole shower algorithms in the Sherpa event generator according to the MC@NLO framework. The interplay between fixed-order real emission contributions and parton shower emissions was studied in detail through variations of the parton shower starting scale. We find large matching uncertainties that exceed the fixed-order uncertainties even in regions of phase space where the fixed-order calculation is well motivated and where parton shower matching effects are expected to be small. Our nominal predictions are in good agreement with the fixed-order result in these regions, however. A comparison to MC@NLO matched results from the literature revealed qualitative differences which are, nevertheless, compatible within the large uncertainties. We observe larger differences in a comparison to POWHEG predictions in the tail of the transverse momentum spectrum, where POWHEG overestimates the fixed-order spectrum by a factor of 2. We find reasonable agreement throughout between our results and those obtained through analytic resummation techniques.