Parton-shower effects in Higgs production via Vector-Boson Fusion

We present a systematic investigation of parton-shower and matching uncertainties for Higgs-boson production via vector-boson fusion. To this end we employ different generators at next-to-leading order QCD accuracy matched with shower Monte Carlo programs, PYTHIA8 and HERWIG7, and a next-to-next-to-leading order QCD calculation. We thoroughly analyse the intrinsic sources of uncertainty within each generator, and then compare predictions among the different tools using the respective recommended setups. Within typical vector-boson fusion cuts, the resulting uncertainties on observables that are accurate to next-to-leading order are at the 10% level for rates and even smaller for shapes. For observables sensitive to extra radiation effects uncertainties of about 20% are found. We furthermore show how a specific recoil scheme is needed when PYTHIA8 is employed, in order not to encounter unphysical enhancements for these observables. We conclude that for vector-boson fusion processes an assessment of the uncertainties associated with a simulation at next-to-leading order matched to parton showers based only on the variation of renormalisation, factorisation and shower scales systematically underestimates their true size.


Introduction
After the discovery of a Higgs boson compatible with the prediction of the Standard Model (SM) of elementary particles at the CERN Large Hadron Collider (LHC) by the ATLAS and CMS experiments [1,2], Higgs physics has entered the era of precision physics. While all measurements completed so far consolidate the SM hypothesis, only a comprehensive analysis of the new boson's properties will reveal whether deviations from the expectation leave room for new physics in the experimentally accessible domain. The precise determination of the Higgs boson's couplings to other elementary particles, spin, and CP properties is thus of paramount importance.
A particularly clean environment for the necessary measurements at the LHC is provided by the vector-boson fusion (VBF) production mode where the Higgs boson is produced by two scattering partons in association with two hard jets (often referred to as tagging jets) in the forward and backward regions of the detector via the exchange of weak massive gauge bosons. Because of the coloursinglet nature of this mechanism, little extra jet activity occurs between the two tagging jets, in the central rapidity region of the detector. These features are of great relevance for separating the VBF signal from QCD background processes that typically exhibit entirely different jet distributions.
Precise measurements can unfold their potential only if matched by equally accurate theoretical predictions. Calculations of the highest accuracy are therefore mandatory in the analysis of VBF data obtained by the experimental LHC collaborations. We note that already now theoretical uncertainties are becoming a bottleneck in Higgs precision studies at the LHC. For instance, in the recent Higgs-combination study by the ATLAS collaboration [3], theory uncertainties are a dominant source of uncertainty in the VBF channel, exceeding statistical and experimental uncertainties. While the QCD corrections to Higgs production via VBF at next-to-leading order (NLO) accuracy have been known for almost 30 years for inclusive cross sections [4] and for almost 20 years for differential distributions with realistic selection cuts [5,6] in the form of flexible parton-level Monte Carlo programs, NLO electroweak corrections have first been presented only later on in Ref. [7] and found to be of almost the same size as the NLO-QCD corrections. Several implementations of VBF-induced Higgs-boson production in programs allowing for a matching with parton-shower programs at NLO-QCD accuracy (in the following referred to as NLO+PS accuracy) are available [8,9,10]. More recently, the fixed order next-to-next-to-leading order (NNLO) QCD corrections have been computed, again first for the fully inclusive case [11,12] and later on differentially [13,14]. These corrections have been found to be small, but not negligible, for differential distributions in the presence of VBF specific cuts. Residual scale uncertainties are tiny at this order in QCD and can be further reduced by the consideration of the next-to-next-to-next-to-leading order (N 3 LO) QCD corrections [15]. Many of the quoted QCD calculations rely on the so-called "VBF approximation", which assumes the absence of colour exchange between the two fermion lines connected by the weak gauge bosons, and neglecting the interferences among H + 2 jet final states produced via s-channel and t-or u-channel topologies, c.f. Fig. 1. At NLO accuracy, the quality of this approximation has been explicitly tested in Ref. [7] and found to be very good once VBF-specific cuts are imposed that force the two tagging jets to be well separated from each other. The impact of different kind of corrections which violate this assumption has been investigated in Refs. [16,12] and recently in Ref. [17]. In all cases, it is found to be of the order of a percent at most.
Ideally, such accurate calculations are provided in the form of public Monte-Carlo programs that can be used by the experimental collaborations directly in their analyses. To make the most of these programs it is important to understand their systematic uncertainties and limitations, for instance due to underlying approximations. In order to provide a systematic assessment of the differences and similarities between commonly used public Monte-Carlo programs designed for VBF-induced Higgs boson production at NLO+PS accuracy, in this article we perform an in-depth comparison of key observables in VBF analyses using realistic input parameters and selection cuts for the respective implementations [8,10,9,18] in the three generators MadGraph5 aMC@NLO [19,20], POWHEG-BOX [8], and HERWIG7 [21,22] VBFNLO+Herwig7/Matchbox [23,24] as well as HJets+Herwig7/Matchbox [25].
We start with a description of the three generators considered in this study in Sec. 2, describe the setup of our analyses in Sec. 3, and discuss the main results of our study in Sec. 4. We conclude with recommendations for the optimal use of the considered generators and a realistic assessment of the associated uncertainties in Sec. 5

MadGraph5 aMC@NLO
MadGraph5 aMC@NLO [19,20] is a meta-code (i.e. a code that generates codes) which makes it possible to automatically simulate arbitrary scattering processes at NLO accuracy in the strong and electroweak couplings, either at fixed order or including matching to parton showers (when one considers only corrections of strong origin), using the MC@NLO method [26]. It employs the FKS subtraction method [27,28] (as automated in MadFKS [29,30]) for the local subtraction of IR singularities. One-loop amplitudes are evaluated by switching dynamically between two integral-reduction techniques, the OPP method [31] or a Laurent-series expansion [32], and tensor-integral reduction [33,34,35]. All such techniques have been automated in the module MadLoop [36], which in turn links Cut-Tools [37], Ninja [38,39], IREGI [40], or Collier [41], together with an in-house implementation of the Open-Loops technique [42]. Uncertainties associated with factorisation and renormalisation scales or parton-distribution functions (PDFs) can be obtained without any approximation thanks to reweighting, at negligible additional CPU cost [43]. The simulation of Higgs production via VBF at NLO-QCD accuracy can be performed with the following commands: For the case of Higgs plus three jets production via VBF, one should simply add a j to the generate command, i.e.: import model loop_qcd_qed_sm_Gmu generate p p > h j j j $$ w+ w-z [QCD] output While results for the first process have been already published in Ref. [9] (although with rather old parton-shower programs), for the second they have been only briefly commented upon in Ref. [19]. In both cases, the $$ syntax forbids W ± and Z bosons to appear in s-channel propagators. Details of the approximation employed in Mad-Graph5 aMC@NLO for VBF-and VBS-type processes can be found in Ref. [44]. In this study we will consider matching to the shower Monte Carlos (SMCs) PYTHIA 8.230 [45] and HERWIG 7.1.2 [46] compiled with ThePEG 2.1.2.

POWHEG-BOX
The POWHEG-BOX [47] is a general framework for the matching of NLO calculations with parton shower programs making use of the POWHEG matching formalism [48,49]. Process-specific components have to be provided on a caseby-case basis. Higgs-boson production via VBF in association with two jets was one of the first processes being implemented in the POWHEG-BOX [8]. More recently, also code for VBF-induced Higgs production in association with three hard jets has been provided [18] being based on the matrix elements of Ref. [50] extracted from the VBFNLO code [23,24]. Both of these implementations rely on the VBF approximation. In this study we will consider matching to the SMCs PYTHIA 8.240 and HERWIG 7.1.4 compiled with ThePEG 2.1.4.

proVBFH
proVBFH v1.1.2 [13] is a public parton-level Monte Carlo program for the calculation of differential distributions for VBF Higgs boson production to NNLO-QCD accuracy in the VBF approximation. It is based on POWHEG's fully differential NLO-QCD calculation for Higgs boson production in association with three jets via VBF [50,18], and an inclusive NNLO-QCD calculation [11], the latter being taken in the structure-function approximation. It achieves differential NNLO-QCD predictions through the projection-to-Born method introduced in [13]. proVBFH includes width effects for the internal W and Z bosons, and neglects fermion masses.

VBFNLO and HJets + Herwig7/Matchbox
The HERWIG7 event generator [51,21,46] features as one of its core components the Matchbox module [10], which can automatically assemble fixed-order and parton shower matched calculations with both the angular ordered [52] and dipole shower algorithms [53], using input from plugins providing matrix elements. The VBFNLO program [23,54] is interfaced as one such module, providing NLO-QCD corrections to the Hjj and Hjjj production processes in the VBF approximation. The HJets library [25] is an alternative module, providing matrix elements and NLO-QCD corrections for the full electroweak Hjj and Hjjj production processes without resorting to the VBF approximation.
In this study we consider the matching using the subtractive matching paradigm; a hard veto scale is imposed on the shower evolution to cut off the parton shower resummation at high transverse momenta. Its central value should reflect the hard transverse momenta at the process of interest, such that the shower evolution will not produce jets with significantly harder transverse momenta. A smearing is applied to the cutoff function, which we choose to be the "resummation" profile studied in more detail in [55,56]. Shower uncertainties are evaluated by varying the hard veto scale, which should reflect the bulk of the uncertainty both in the soft region and in regions which will be improved through the NLO matching.

Recoil schemes in PYTHIA8
By default PYTHIA8 employs a global recoil scheme for the generation of initial-state radiation. While this is certainly a valid approach when the underlying hard scattering does not have a colour flow between initial and final states, e.g. for colour-singlet production, it leads to inconsistencies when considering, for instance, Deep Inelastic Scattering (DIS), where the colour flow is only between an initial-state quark and a final state quark. This was discussed in Ref. [57] and a new dipole approach was introduced for initial-state radiation to better describe processes with initial-final colour flow. Since VBF can essentially be viewed as a double-DIS process where there is no QCD cross-talk between the two incoming protons, that discussion is also highly relevant here. It is known that in the VBF approximation a gluon emitted from one quark line cannot attach to the other quark line. It is therefore not very physical to distribute the recoil of such an emission over the entire event since such a prescription would destroy the relation between the kinematics and the soft radiation pattern. Instead one would expect the recoil to be along the quark line where the gluon emission took place. We therefore find it worth investigating the two different recoil schemes inside PYTHIA8 in this study. The dipole recoil scheme can be used directly with the POWHEG-BOX, whereas it is not currently possible with Mad-Graph5 aMC@NLO as the shower counterterms have been derived assuming a global recoil 1 . In the following we will therefore only show results using the dipole approach and the POWHEG-BOX. For the default (global) recoil scheme we show results obtained with both MadGraph5 aMC@NLO and the POWHEG-BOX. The inadequacy of a global-recoil scheme has been discussed for VBS processes in Ref. [44], and for Z-boson production via VBF in Ref. [59].

Input parameters
We consider proton-proton scattering at the LHC with a centre-of-mass energy of √ s = 13 TeV. For the PDFs of the proton we use an NNLO set with five massless flavours, PDF4LHC15 nnlo 100 pdfas [60], as provided by the LHAPDF6 library [61] (identifier LHAPDF ID=91200) with the corresponding strong coupling, α s (M Z ) = 0.118.
For the masses and widths of the particles entering our calculation the following values are used: As electroweak (EW) input parameters we use M W , M Z , and the Fermi constant, Other EW parameters, such as the EW coupling α and the weak-mixing angle, are computed therefrom via tree-level EW relations. The Cabibbo-Kobayashi-Maskawa matrix is assumed to be diagonal, i.e. mixing effects between different quark generations are neglected.
The renormalisation scale, µ ren , and the factorisation scale, µ fac , are identified with ξ ren µ 0 and ξ fac µ 0 , where the parameters ξ ren and ξ fac are to be varied between 1/2 and 2, and the central scale µ 0 , obtained from is computed from the mass and transverse momentum p T,H of the Higgs boson event by event. We do not include effects of hadronisation or underlying events in our simulations.

Selection cuts
For the simulation of VBF events we employ a set of cuts that ensure that the considered fiducial volumes can suitably be accurately described despite the approximations used in (some of) the generators of this study, such as the VBF approximation that only works in a setup that disfavours Higgs-strahlung topologies.
In order to define a H + n jets event we require the presence of at least n jets, obtained from partons via the anti-k T algorithm [62] using the FastJet package [63] with a distance parameter R. Unless specified otherwise, the value of R is set to 0.4. The thus produced jets need to exhibit a minimum transverse momentum and be located within the pseudo-rapidity range covered by the detector, The hardest two jets fulfilling this criterion are called "tagging jets". These two tagging jets are furthermore required to be located in opposite hemispheres of the detector, well separated in rapidity, and exhibit a significant invariant mass,

Numerical analysis
In the following we will present the numerical results of our study. We will first discuss uncertainties specific to the individual generators. In the second part of this section, we will compare representative predictions of the individual generators with each other.

Results from MadGraph5 aMC@NLO
We now discuss results for VBF obtained with MadGraph-5 aMC@NLO, and elaborate on effects due to the specific SMC employed and to the shower starting scale, on top of the usual estimate of theoretical uncertainties from the variation of the hard (renormalisation and factorisation) scales. As SMCs, we consider the angular-ordered HERWIG7 generator and PYTHIA8 with a global-recoil scheme. Concerning the shower starting scale Q sh , it assigns (on an event-by-event basis) the maximum hardness of the radiation that the shower can generate in terms of the specific evolution variable, and is computed from a reference shower scale µ sh . In general, one has Q sh = µ sh for the so-called H-events, while for the S-events Q sh is generated from a probability distribution of which µ sh is the upper endpoint 2 . In order to assess the sensitivity of VBF observables on the shower scale, we choose to present results where either µ sh is not changed from its default value, or where it is halved 3 . All plots, except those depicting properties of the third jet, which will be presented later, have the following layout: four histograms are displayed, with predictions obtained using PYTHIA8 (HERWIG7) in blue (red). Solid (dashed) histograms correspond to the default (halved) reference shower scale. In the inset, we show the bin-by-bin ratio over the prediction matched to HERWIG7 with nominal shower scale. A blue band, corresponding to the hard-scale variations (the renormalisation and factorisation scales are varied independently by a factor of two around the central value giving rise to a nine-point variation) is displayed for the prediction matched to PYTHIA8 with the nominal shower scale. The first observable we consider is the exclusive 4 jet multiplicity, in the left panel of Fig. 2. When looking at this figure, one should bear in mind that the two-jet bin is the only bin with genuine NLO accuracy. The threejet bin is only LO accurate, while higher multiplicities of jets are entirely due to the SMC. A consequence of this is the agreement among predictions in the two-jet bin, where predictions lie within 10% of each other, with those matched to PYTHIA8 predicting a lower rate than those with HERWIG7. In the three-jet bin, on the other hand, we observe large discrepancies, not covered by the hard-scale uncertainty: the predictions matched with PYTHIA8 exhibits a 60% excess with respect to the one matched with HERWIG7. Such a large effect is due to the global recoil scheme employed by PYTHIA8 in order to be consistent with the matching in MadGraph5 aMC@NLO, which is not suitable for VBF/VBS-type processes, c.f. our discussion in Sec. 2.5. For higher-multiplicity bins discrepancies and scale uncertainties become huge. Finally, we remark that predictions matched with PYTHIA8 display a more pronounced sensitivity on the shower starting scale, while for HERWIG7 such a dependence is very small. The next observable we consider is the transverse momentum distribution of the Higgs boson, in the right panel of Fig. 2. This observable displays an excellent agreement among all predictions, with discrepancies of few percents at most, a behaviour which is common for observables inclusive in the number of jets: indeed, the differences in the two-and three-jet bins described before tend to compensate almost exactly. We have verified that this applies for many other NLO-accurate observables, such as those related to the first and second tagging jet. As representative ones, we show the transverse momentum of the second tagging jet and the rapidity separation of the two tagging jets in Fig. 3. We remark that the dependence on the renormalisation, factorisation, and shower scales for these observables is very small, with the exception of the rapidity separation at large rapidities, comparable to the differences among predictions employing the two parton showers.
We now turn to observables related to the third jet, in particular the transverse momentum and rapidity distributions, respectively, shown in Figs. 4 and 5. In order to reach NLO accuracy also for these observables, we additionally show predictions for the production of a Higgs boson in association with three jets via VBF at NLO+PS accuracy, both matched with HERWIG7 (orange) and with PYTHIA8 (green). The line pattern (solid or dashed) has the same meaning as above. For the sake of better readability, we show two panels for each observable. In the left (right) panel, we show the four predictions for the production of a Higgs boson plus two (three) jets via VBF and the one for the production of a Higgs boson plus three (two) jets matched with HERWIG7 using the default shower scale. In the inset we show the respective ratios over the prediction for the production of a Higgs boson plus two (three) jets via VBF, matched with HERWIG7 and with nominal shower scale. The plotting range is different in the inset of the left and right panels. The predictions for the rapidity distribution of the third jet for the production of a Higgs boson plus two jets via VBF at NLO+PS show that the origin of the excess observed in the jet multiplicities when matching to PYTHIA8 mainly comes from jets in the central region, as a consequence of the global-recoil scheme. The same effect is rather flat in the transverse momentum spectrum. It is worth to observe that reducing the shower scale is not sufficient to cure this behaviour, and that the renormalisation and factorisation scale variations fail to cover differences among the shower generators. Indeed, such a behaviour is unphysical, which can be understood by looking at the predictions for the production of a Higgs boson plus three jets at NLO+PS accuracy. The difference among various predictions is now reduced to the 10% level or below, thanks of the better perturbative description of these observables. It is important to stress that the PYTHIA8 predictions still employ a global-recoil scheme, in accordance with the needs of the matching in MadGraph5 aMC@NLO. It is also worth to notice the impact of the correction (in the case of HERWIG7) when passing from an LO description (Higgs plus two jets via VBF) to an NLO one (Higgs plus three jets via VBF): while no visible effect can be appreciated in the transverse momentum spectrum, looking at the rapidity one can see how the NLO corrections tend to enhance central rapidities and deplete larger ones (|η| > 3.5).
In conclusion, supported by the results presented in this section and given the impossibility to employ PYTHIA8 in conjunction with a dipole-recoil scheme within Mad-Graph5 aMC@NLO, we strongly advise to use MadGraph5 -aMC@NLO only in conjunction with HERWIG7 for the simulation of VBF.

Results of the POWHEG-BOX
In the POWHEG-BOX an assessment of the intrinsic uncertainty related to the POWHEG matching procedure is possible by a variation of the so-called hdamp parameter. This parameter governs the splitting of the full real-emission contribution R into a singular part, R s , that enters into the Sudakov form factor and a regular part, R f , according to where p T denotes the transverse momentum of the hardest parton of the real-emission contribution and h is a parameter that can be set by the user. We explore the matching uncertainty accessible via the hdamp parameter by considering the three cases h = ∞ (i.e. no damping), h = M H , and h = m min j 1 j 2 = 600 GeV. We show plots using PYTHIA8 as the SMC with the dipole recoil strategy [57].
Naively, one would expect observables related to the hard jets that are not very sensitive to soft emission to be less affected by the choice of the hdamp parameter than distributions related to the sub-leading jets. To assess this expectation, in Fig. 6 we show examples of both types of observables in the VBF Hjj process. The invariant-mass distribution of the two tagging jets is completely insensitive to the value of h. However, the same holds true also for the transverse-momentum distribution of the third jet over the entire range considered, where larger effects might be expected. This finding clearly indicates that the VBF process considered here is quite insensitive to the actual form of the Sudakov form factor used for the POWHEG-BOX simulation. We remark that, consequently, the choice of the hdamp parameter has little impact on the numerical stability and CPU requirements of the program. We will therefore use the value h = ∞ (corresponding to hdamp = 1) as a default.
While the dependence of predictions obtained with the POWHEG-BOX on hdamp obviously is very small, another source of generator-specific uncertainty is constituted by the choice and settings of the SMC, the POWHEG-BOX is matched to. To explore this effect we present a systematic comparison of NLO+PS predictions obtained with PYTHIA8 (both default and dipole recoil scheme, c.f. Sec. 2.5), angular ordered HERWIG7, and fixed-order results at NNLO-QCD accuracy obtained with the proVBFH program. We expect only a small impact of the SMC choice on observables with little sensitivity to soft radiation effects, such as the transverse momenta of the tagging jets and related distributions. Indeed, as illustrated by Fig. 7, the transversemomentum distribution of the second tagging jet is very stable with respect to the choice of SMC, and indeed the NLO+PS simulation provides a very good approximation for the NNLO prediction. Small differences are also observed in the rapidity separation of the two tagging jets, shown in the right-hand-side of Fig. 7. We notice, however, that in this case the results obtained with the dipole recoil scheme in PYTHIA8 lie clearly above the HERWIG7 results, while the default version of PYTHIA8 resembles the HERWIG7 predictions in the region of highly separated jets, but reproduces the PYTHIA8 results in the dipole scheme for smaller rapidity separations.
Much more pronounced differences between the various SMC choices are found for distributions related to the subleading jets. Figure 8 shows the transverse-momentum distribution of the system formed by the Higgs boson and the two tagging jets, which reflects the transverse momentum of the remaining objects produced in the scattering process, in particular the non-tagging jets. Since such subleading jets in the Hjj simulation can only be accounted for by the real-emission matrix elements or parton-shower emission they are only described at leading order or partonshower accuracy. In the tail of the p T,H,j1,j2 distribution, the PYTHIA8 default results by far exceed the reference results constituted by the NNLO prediction, while no such large differences are observed in the HERWIG7 and PYTHIA8 results using the dipole recoil scheme.
A variable particularly suitable to indicate the relative position of the third jet with respect to the centre of the tagging-jet system is constituted by the so-called Zeppenfeld variable, defined as For small values of z j 3 the third jet is right in between the two tagging jets, while larger z j 3 values correspond to more peripheral configurations. The z j3 distribution helps to understand where the large differences between the various SMC simulations stem from. Obviously, the PYTHIA8 default scheme produces an abundance of radiation for small values of z j3 , i.e. in between the two tagging jets.

VBFNLO and HJets + Herwig7/Matchbox
Within the setup using the HERWIG7 interface to VBFNLO and HJets we perform the subtractive, MC@NLO-type matching and assess the uncertainties by varying the hard scale of the shower evolution as well as the factorisation and renormalisation scales of the hard process. For a detailed discussion of these uncertainties see [55,65], where VBS processes have been considered as well. We also investigate the difference between the default, angular orderedq shower, as well as the dipole-type evolution which is available as an alternative module. Since the HJets module [25] implements the calculation without any VBF approximation, we can perform a comparison to VBFNLO, which resorts to the VBF approximation that is also used in the POWHEG-BOX and MadGraph5 aMC@NLO generators. We find quite similar results of the showering in between the two HERWIG7 shower modules, as well as similar variations and stability with respect to the fixed order input. We first compare the VBF approximation for both a tight and a loose cut setup with subsequent parton showering, including the variations from the renormalisation and factorisation scales. The tight setup is defined by the cuts of Sec. 3.2, while for the loose setup we relax the selection to with all other cuts identical to the general setup. Examples are depicted in Fig. 9, where we generally find a large discrepancy between VBFNLO and HJets for the third jet in a loose setup, and a very good agreement in between the two for a tight VBF selection. Similar findings at fixed order also apply to the third jet distributions, see [66]. Within a tight VBF selection, the shower uncertainties in the NLO matched case are at the few-percent level for observables involving the hardest three jets, but can still be significant for higher jet multiplicities, something which we exemplify in Fig. 10, where we include the minimum rapidity difference of the third jet with respect to the tagging jets, defined by where x * j3 receives a minus sign if the third jet is outside the dijet window, i.e. if z j3 > 0.5. We also show the dijet invariant mass distribution.

Comparison of different generators
Having investigated variations within the individual SMCs we now turn to a study of the three generators in the recommended default setup. A summary of the setups used with the three different generators is given in Tab. 1. Given the above discussion we show results for MadGraph5 aMC-@NLO interfaced to HERWIG7, the POWHEG-BOX interfaced to PYTHIA8 using the dipole recoil strategy, and VBFNLO+-Herwig7/Matchbox. All three generators use the VBF approximation, and have been checked to agree within statistical uncertainties when run at fixed-order (at the inclusive and differential level). Hence we expect any disagreement to arise only from differences in matching procedure and shower details rather than the fixed-order matrix elements for the hard scattering. We recall that we do not include hadronisation or underlying event effects in the comparison. In Fig. 11 we show the typical VBF observables; tagging jet rapidity separation, ∆η j1,j2 , and invariant mass, m j1,j2 , for MadGraph5 aMC@NLO (blue), POWHEG-BOX (green), and VBFNLO+Herwig7/Matchbox (orange). We also show the fixed order NNLO-QCD prediction obtained using pro-VBFH (black). The plot shows a spread in predictions of less than 10%. Both POWHEG-BOX and MadGraph5 aMC@NLO show the same shape distortion with respect to proVBFH although they have different normalisation. VBFNLO+Her-wig7/Matchbox, on the other hand, exhibits a different slope behaviour in both observables with respect to the other two generators.
There are also some differences between the three generators when considering more inclusive observables. However in this case the discrepancies are mostly due to differences in normalisations. To illustrate that point, in Fig. 12 we show the transverse momentum of the Higgs boson and of the first tagging jet in the event. All three generators agree within 10% and have very similar shapes. In particular, all three generators are comparable in shape with respect to the fixed order NNLO-QCD prediction.
Lastly we show a comparison of the Zeppenfeld variable z j 3 and the exclusive jet multiplicity in Fig. 13. We remind the reader that all three considered generators have LO accuracy for three-jet observables and pure shower accuracy for observables with more than three jets. Although there are larger differences between the generators for z j 3 , of the order of 20%, they have fairly similar shapes up to about z j 3 0.8 and, in particular, none of the predictions exhibits a large excess in the small z j 3 region. For the exclusive jet cross section it is clear that matched calculations predict a much smaller number of jets than the fixed order prediction in the three and four jet bins. They do, however, agree amongst each other at the 10% level for the 2, 3 and 4 jet rates. The discrepancy with respect to the fixed order prediction is related to soft radiation produced by the shower that is lost outside of the rather narrow jet cone.

Jet radius dependence
In this section we consider the dependence of the VBF cross sections on the jet radius R after showering, but without any hadronisation or underlying event for which we expect a parametrically different dependence on the jet radius. From parton showering, and higher order corrections in general, we expect a leading log(1/R) dependence, which has previously been studied for VBF processes in [67], and for more general processes involving hard jets the interplay with scale choices and variations at fixed order has also been investigated [68]. We show some of the results in Fig. 14. While we have not attempted to perform any fit of the R-dependence, the general pattern we see is that after parton showering leading, as well as nextto-leading order matched predictions show a similar, and significant R dependence. This dependence does not only affect the normalisation of the cross section due to the jet selection criteria, but also the shapes even for inclusive distributions like the Higgs boson transverse momentum. A comprehensive discussion of the jet radius dependence needs not only to include a study of the behaviour of NLO QCD corrections, but also to include the impact of hadronisation and multi-parton interactions. Preliminary results for investigating the jet radius dependence at NNLO have also been reported in [69].

Recommendations and conclusions
In this work we performed a quantitative investigation of parton-shower and matching uncertainties for the production of a Higgs boson plus two jets via VBF. The relevance of such a study is supported by the fact that, already in analyses based only on part of the data taken during Run II of the LHC, for VBF the dominant source of uncertainties are theoretical ones. Improving on Higgs analyses in the VBF channel thus crucially requires a quantitative understanding of the tools used for the simulation of Higgs production via VBF.
In the study of matching uncertainties, we found that, within a single generator and SMC, theoretical uncertainties estimated by the usual renormalisation and factorisation scale variations, possibly supplemented by variations in a variable that controls the shower hardness (shower starting scale for MadGraph5 aMC@NLO or hdamp for the POWHEG-BOX), turn out to be small, hardly above the fewpercents figure. This also applies to the hard shower scales variations in HERWIG7, which can become more significant if properties of the third jet are probed. However we showed that the differences among predictions obtained with different SMCs can be more significant, easily exceeding the aforementioned estimate of theory uncertainties. For observables described at NLO-QCD accuracy, these differences are at the 10% level. However, they are mostly due to normalisation effects, while shapes of distributions are described to an even better accuracy when the various NLO+PS programs and the NNLO result are compared. For LO-accurate observables, differences turn out to be much larger, but not always physical. A prominent example is the description of third-jet observables when PYTHIA8 is employed with a global-recoil scheme, which gives a huge enhancement in the central-rapidity region. Such an enhancement has been proven to be unphysical by looking at an NLO-accurate description of the same variable, where it disappears. Taking this fact into account, uncertainties for third-jet observables can be quantified in the 20% domain.
As a consequence, we recommend against using PYTHIA8 with a global-recoil scheme for VBF. Instead one should change the recoil scheme to the dipole one when this is compatible with the matching (i.e. with POWHEG-BOX). When this is not (yet) the case (with MadGraph5 aMC@NLO) one should use an entirely different SMC like HERWIG7, which performs the matching internally and uses recoil schemes which respect the colour flow information of the hard process either through the initial conditions to the angular ordered evolution in case of the defaultq-shower or the nature of the alternative dipole shower algorithm, which lead to comparable results.
We conclude that, within the typical VBF phase space, all the programs considered in this study yield reliable results. However, we remind the reader that because of the VBF approximation used in most of the considered generators, valid predictions can only be expected after appropriate selection cuts are employed. As far as VBF Higgs processes are concerned, the HJets plugin to HERWIG7 can provide accurate predictions for H + 2 jet and H + 3 jet final states at NLO QCD without resorting to the VBF approximation and we have used this as an explicit check to demonstrate good agreement within the VBF selection region.
We also stress that a comprehensive study of uncertainties for VBF predictions necessarily needs to include the effects of multi-parton interactions, colour reconnection and hadronisation. The impact of these effects will vary largely with the jet radius and need to be confronted with the perturbative variations in order to obtain a global picture. We leave such a study to future work.
We have included the RIVET [70] analysis used in this study with the ancillary files of the arXiv submission for anyone interested in reproducing our results.  Fig. 3. Same as in Fig. 2, for the transverse momentum of the second tagging jet (left) and for the rapidity separation of the two tagging jets.  Fig. 5. Same as in Fig. 4, for the rapidity of the third jet.  Fig. 7. Transverse-momentum of the second tagging jet (left) and separation of the two tagging jets in pseudorapidity (right) within the cuts of Eqs. (3)-(4) at NNLO, and at NLO+PS accuracy using the POWHEG-BOX matched with HERWIG7 and PYTHIA8 using two different recoil schemes. No hadronisation effects are taken into account.   Fig. 14. The jet radius dependence illustrated for the pseudorapidity difference between the tagging jets, and the transverse momentum of the Higgs boson. Inclusive quantities also show a significant dependence on the jet radius due to selection criterion involving jets.