Parton-shower effects in electroweak $WZjj$ production at the next-to-leading order of QCD

We present an implementation of $WZjj$ production via vector-boson fusion in the POWHEG BOX, a public tool for the matching of next-to-leading order QCD calculations with multi-purpose parton-shower generators. We provide phenomenological results for electroweak $WZjj$ production with fully leptonic decays at the LHC in realistic setups and discuss theoretical uncertainties associated with the simulation. We find that beyond the leading-order approximation the dependence on the unphysical factorization and renormalization scales is mild. The two tagging jets are furthermore very stable against parton-shower effects. However, considerable sensitivities to the shower Monte-Carlo program used are observed for central-jet veto observables.


Introduction
Vector boson scattering (VBS) processes provide particularly promising means for probing the gauge structure of the Standard Model's electroweak sector. Processes of the type V V → V V (with V denoting a W ± or Z boson) involve pure weak gauge-boson scattering contributions as well as contributions featuring Higgs exchange contributions. In the framework of the Standard Model (SM), well-behaved cross sections respecting unitarity conservation up to the highest energy scales require a subtle interplay of triple and quartic couplings of the gauge and Higgs bosons. Deviations in measurements from theoretical predictions of gauge-boson scattering processes could thus point to anomalous effects in triple and quartic gauge boson couplings, or other types of new physics in the electroweak sector.
At hadron colliders such as the CERN Large Hadron Collider (LHC), on-shell gauge boson scattering processes of the type V V → V V are not directly accessible. Instead, VBS processes are most conveniently studied by means of the purely electroweak (EW) production processes pp → V V + 2 jets, where two protons scatter by the exchange of weak gauge bosons that in general exhibit non-vanishing virtuality. In addition to the two final-state gauge bosons the scattering protons produce two so-called tagging jets that are typically located in the forward regions of the detector. Experimentally, such processes can be separated rather well from strong production processes with the same final state, because of the characteristic distribution of the tagging jets that tend to be produced at large invariant mass via the EW production mode, while they are mostly produced close in rapidity with small invariant mass in QCD-dominated background processes. The EW production mode also distinguishes itself from the QCD induced one by generating significantly fewer jets in the detector volume between the two tag jets [1]. This characteristic can be used to veto such activity and hence further reduce the QCD background. To best exploit the characteristic VBS signature, one considers fully leptonic decays of the gauge bosons, resulting in final states with four leptons that are typically located between the two tagging jets. At order O(α 6 ) such final states cannot only stem from genuine weak-boson scattering contributions involving a V ( ) V ( ) → V V → 4 leptons topology, but also from other electroweak contributions not involving resonant V bosons. Only the sum of all diagrams giving rise to a specific final state with four leptons and two jets in a gauge invariant manner is meaningful to consider.
In this work, we focus on the EW production of a final state with one neutrino, three charged leptons, and two jets, ν + − + jj (with , = µ, e), that is dominated by W + Z VBS contributions. First experimental results on this production mode were reported for an LHC center-of-mass energy of √ s = 8 TeV using a data sample with an integrated luminosity of 20.3 fb −1 by the ATLAS collaboration in Ref. [2]. Very recently, the ATLAS and CMS collaborations have presented first measurements of this production mode at √ s = 13 TeV, using data samples collected in 2015 and 2016 corresponding to integrated luminosities of 36.1 fb −1 and 35.9 fb −1 , respectively [3,4]. Based on leading-order (LO) simulations, ATLAS reports measured cross sections slightly above the theoretical prediction [3], whereas in a slightly different setup CMS finds agreement of their measurements with the SM prediction [4].
Although current measurements are limited by statistics, the situation will change in the future when more data become available. In particular it will be of great importance that the available theoretical predictions can match the expected experimental precision [5]. It is therefore necessary to have a good understanding of the accuracy of the particular theoretical prediction when comparing to experimental data. Two studies investigating the accuracy of various VBS predictions have recently been published. The first study, a detailed analysis of the various contributions to the W ± Zjj production mode at LO, was presented in Ref. [6] using the Monte-Carlo programs MoCaNLO+RECOLA [7,8], SHERPA [9,10], MadGraph5 aMC@NLO [11], and VBFNLO [12]. The fixed-order parton-level comparison revealed good agreement between the various theoretical predictions, even though only VBS-induced contributions where considered in VBFNLO while all amplitudes contributing at order O(α 6 ) have been retained in the other generators. However, when the LO generators were combined with parton-shower tools, significant scale uncertainties and larger discrepancies between the various predictions were observed, emphasizing the need for higher order QCD corrections at matrix-element level. In particular observables related to the third jet, which is not present at LO and is therefore generated purely by the parton shower, showed very large discrepancies. Such discrepancies are expected to decrease when predictions at next-to-leading order (NLO) of QCD are used instead.
A similar study of the related VBS process of W ± W ± jj production was presented in Ref. [13]. This process has received significantly more theoretical attention than the other VBS processes due to its simpler structure. As such this study did not only present results at LO, but discussed also the matching of NLO-QCD calculations with parton showers (PS). Although it was found that the description of all observables improves when using NLO+PS tools compared to LO+PS it became clear that remaining uncertainties related to the parton-shower generator can be significant and should be accounted for systematically.
While NLO-QCD corrections to VBS-induced ν + − + jj production at the LHC have been known for quite some time [14] and available in the form of a flexible partonlevel Monte-Carlo program in the VBFNLO package [12], the matching of an NLO-QCD calculation with parton-shower programs has not been presented so far for this channel. Work in that direction has already been performed for the related W + W + jj, W + W − jj, and ZZjj VBS production modes [15][16][17][18]. For the W Zjj mode, however, no public implementation was available up to now. Here we close this gap by implementing the calculation of [14] in the framework of the POWHEG BOX [19], a tool for the matching of dedicated NLO-QCD calculations with public parton-shower generators such as PYTHIA [20,21] or HERWIG [22,23], making use of the POWHEG method [24]. Doing so in a public framework, we provide the tools for performing analyses of EW W Zjj production at NLO+PS accuracy and thus avoiding theoretical uncertainties due to an inaccurate treatment of the hard scattering process as inherent to a mere LO simulation.
Technical details of our implementation are briefly summarized in Sec. 2. In Sec. 3 we present results of a detailed numerical analysis with particular emphasis on the impact different parton shower generators have on phenomenological predictions. Our conclusions are provided in Sec. 4.

Details of the implementation
In order to implement the electroweak W Zjj production process in the framework of the POWHEG BOX we closely followed the strategy applied to the related VBS processes W + W + jj [15], Zjj [25], W + W − jj [16], and ZZjj [17]. We extracted the relevant matrix elements for electroweak W Zjj production including the fully leptonic decays of the W and Z bosons at LO (i.e. at order O(α 6 )) and at NLO-QCD from VBFNLO and adapted them to comply with the format required by the POWHEG BOX V2. We retained all approximations inherent to the calculation of [14] that forms the basis of the VBFNLO implementation. In particular, at LO we considered (anti-)quark scattering contributions giving rise to a final state with two (anti-)quarks, a neutrino, and three leptons, resulting in a ν e e + µ − µ + jj signature. While for the matrix elements squared t-channel topologies are taken into account as well as u-channel topologies, interference contributions of tchannel with u-channel diagrams are neglected. Contributions from s-channel-induced production modes are considered as part of a three-gauge-boson production process and entirely disregarded. The neglected contributions are small in the phase-space region where VBS processes are typically investigated experimentally, with widely-separated jets in the forward-and backward regions of the detector. For instance, at LO they result in only a 0.6% deviation from the full calculation in the realistic setup explicitly investigated in Ref. [6]. Throughout, the Cabibbo-Kobayashi-Maskawa matrix is assumed to be diagonal, and contributions with external top or bottom quarks are disregarded. In the following we will refer to electroweak ν e e + µ − µ + jj production in proton-proton collisions within the afore-mentioned approximations as VBS-induced W Zjj production.
Similar to the previously considered Zjj, W + W − jj, and ZZjj VBS processes, the electroweak W Zjj production cross section in the fully inclusive case exhibits various types of singularities at Born level. While such singularities can easily be avoided by the application of phase-space cuts in a fixed-order calculation, the inclusive sampling of the VBS phase space required by the POWHEG BOX calls for a more subtle treatment. We achieve that goal in a two-step procedure: First, we introduce a technical cut to remove contributions with singularities due to the exchange of a t-channel photon of low virtuality, Q 2 min ≤ 4 GeV 2 . Such contributions are entirely negligible after analysis cuts of p jet T ≥ 40 GeV are imposed on the two tagging jets, and can therefore be disregarded already at generation level. Second, an improvement in the efficiency of the phase-space integration can be achieved by applying a so-called Born-suppression factor F that dampens the integrand in singular regions of phase space. We employ a Born-suppression factor of the form with the p T,i (i = 1, 2) denoting the transverse momenta of the final-state (anti-)quarks of an underlying Born configuration Φ n , and Λ a technical damping parameter set to 10 GeV.
In addition to singularities due to the exchange of low-virtuality photons in the tchannel, in ν e e + µ − µ + jj production processes divergences can occur when the µ − µ + pair stems from the decay of a quasi on-shell photon. Since in our numerical analyses we will not consider QED shower effects that may modify the kinematics of the final-state leptons, we remove such configurations already at generation level by an explicit cut on the invariant mass of the same-type lepton pair, retaining only events with To validate our implementation, we performed several checks. We compared results for the tree-level matrix elements squared used in the POWHEG BOX with equivalent expressions obtained by MadGraph [26] for individual phase-space points and found full agreement at double-precision accuracy. Additionally, we compared parton-level results for cross sections and differential distributions at LO and at NLO-QCD with results obtained by VBFNLO, again finding full agreement within the numerical accuracy of the two calculations.
The computer code we developed will be made available in the public version of the POWHEG BOX V2, accessible at the webpage http://powhegbox.mib.infn.it/.

Phenomenological results
In order to demonstrate the capabilities of the developed code we present phenomenological results for a few selected setups inspired by realistic experimental studies.
In each case, we consider proton-proton collisions at the LHC with a center-of-mass energy of √ s = 13 TeV. We use the PDF4LHC15 parton distribution function (PDF) set at NLO [27] as provided by the LHAPDF library [28] (LHAPDF ID 90000). For the definition of jets we employ the anti-k T algorithm [29,30] with R = 0.4, as implemented in the FASTJET 3.3.0 package [31]. The mass and width of the Higgs boson are set to m H = 125.18 GeV and Γ H = 0.00407 GeV, respectively. As EW input parameters we chose the Fermi constant, G µ = 1.16638 × 10 −5 GeV −2 , and the masses of the W and Z bosons, m W = 80.379 GeV and m Z = 91.188 GeV, respectively. The widths of the massive gauge bosons are set to Γ W = 2.085 GeV and Γ Z = 2.4952 GeV.
We provide results for a setup inspired by the ATLAS analysis of Ref.
[3], dubbed ATLAS-loose, and for a scenario with more stringent cuts, referred to as CMS-tight, following the CMS analysis of Ref.
[4]. In order not to overload the discussion in the main body of this article, in this section we focus on results in the CMS-tight scenario. For reference, we show the same distributions obtained in the ATLAS-loose scenario in Appendix A.
The CMS-tight scenario comprises the following stringent selection cuts on the jets with p T,j > 50 GeV , |y j | < 4.7 . (3) The two hardest jets satisfying these criteria are called "tagging jets", and are furthermore required to fulfill the additional requirements supplemented by the separation cuts and cuts specifically distinguishing between the hardest and second-hardest muon of each event, the electron, and the missing transverse momentum which in our calculation corresponds to the momentum of the neutrino with pseudorapidities of and cuts on the invariant masses of the µ + µ − and the three-lepton systems, respectively, In addition, the location of the three-lepton system in pseudo-rapidity, η 3 , with respect to the two tagging jets is required to fulfill We note that when considering distributions of the third jet in the following, we only consider events with p T,j 3 > 10 GeV, as softer jets cannot be reliably identified. For the discussion of the transverse-momentum distribution of the third jets, naturally no such requirement is imposed. For clarity, whenever we restrict the considered range of p T,j 3 we will specifically indicate that in the respective figures. The factorization and renormalization scales, µ F and µ R , for each phase-space point are identified with the geometric mean of the transverse momenta of the two final-state partons of the underlying Born configuration, For studying the dependence of the various results on these scales, we introduce the two parameters ξ F and ξ R such that µ F = ξ F µ 0 and µ R = ξ R µ 0 , and vary them independently. In order to explore the impact of parton-shower effects on experimentally accessible distributions and to investigate their uncertainties, we consider three different partonshower Monte-Carlo generators (SMCs): PYTHIA 6.4.28 [20], PYTHIA 8.230 [21], and HERWIG 7.1.4 [23]. For PYTHIA8 we use the Monash Tune (Tune:pp = 14) [34], for PYTHIA6 we use the Perugia 2012 M8LO tune [35] and for HERWIG7 we use the default tune. In particular, PYTHIA8 and HERWIG7 offer the possibility of a local recoil dipole shower as an alternative to their respective default global recoil showers. The local recoil schemes have been shown to describe data for deep-inelastic scattering [32] well, and have also been found to be advantageous for the theoretical description of VBS W ± W ± jj production [18] and for t-channel top-quark production [33]. Both of these processes are similar to the VBS W Zjj process from the point of view of the parton shower as they share the same color flow between the initial and final state partons at LO. In order to explore the possible advantages of a dipole shower in the description of VBS W Zjj production we In order to illustrate the scale uncertainty of our predictions, for the PYTHIA6 reference results we have performed a seven-point variation of the scale parameters ξ F , ξ R -that is the two parameters ξ F and ξ R have been varied independently from 0.5 via 1 to 2 dropping configurations where ξ F and ξ R differ by a factor larger than two. The resulting range of results is indicated by the blue shaded band in the respective distributions. We notice that differences due to the four parton shower algorithms used are of the order of or even smaller than the theoretical uncertainty due to scale variation, which amounts to about 5 %. As expected, the choice of SMC does not have a significant impact on the normalization and shape of the considered distributions. Much larger dependencies on the used SMC are however found for distributions of the non-tagging jets. We note that in the NLO-QCD calculation for VBS W Zjj pro-  Figure 2: Transverse-momentum distribution of the hardest positron (left) and azimuthal-angle separation of the two tagging jets (right) in the CMS-tight scenario. In the lower panels the ratio of the respective distribution to the PYTHIA6 reference result is shown. In each case the blue bands indicate the scale uncertainty of the PYTHIA6 simulations, statistical uncertainties are denoted by error bars. duction the third jet enters only via the real-emission corrections and thus is accounted for effectively only with LO accuracy. In the NLO+PS calculation a third jet can also stem from the parton shower. Figure 3 shows the transverse-momentum and the rapidity distributions of the third jet in the CMS-tight scenario.
While the scale uncertainty of these distributions still amounts to less than 10 %, the differences among the various SMCs can be significantly larger. In the transversemomentum distribution, the differences between predictions obtained with different versions of PYTHIA and HERWIG are very large in the low-p T region, where the fixed-order calculation is not reliable because of infrared divergent configurations. The Sudakov factor of the NLO+POWHEG calculation dampens the would-be divergences resulting in a well-behaved shape of the distribution even in the low-p T region. In the tail of the distribution the differences between the various SMCs disappear, making predictions in the range above p T,j 3 20 GeV, where jets can be identified experimentally, more reliable.
The differences between the individual SMCs very much affect also the rapidity distribution of the third jet. In particular, PYTHIA simulations with the default recoil shower tend to fill the central-rapidity region much more than HERWIG7 with its default recoil shower, while PYTHIA8 with the dipole shower resembles the HERWIG7 result. These SMC dependencies should be kept well in mind when using so-called central-jet vetoing (CJV) techniques.
An enhancement of the VBS signal with respect to QCD-dominated background processes can be achieved by exploiting the typical distribution of sub-leading jets depending on the nature of the production process. QCD-induced processes typically exhibit large jet activity in the central region of rapidity. In VBS processes, on the other hand, extra jets tend to be produced close to the tagging jets in the forward and backward regions of the detector. Thus, vetoing events with hard central jets (emerging in addition to the two tagging jets) helps to reduce the impact of unwanted background processes on the VBS signal.
The relative position of the third jet with respect to the center of the tagging-jet system is accounted for by the so-called Zeppenfeld variable, which approaches zero when the third jet is located just in the center of the two tagging jets, and one half, if it is close to one of the tagging jets, while values of z j 3 larger than 0.5 indicate that the third jet is located outside of the rapidity interval spanned by the two tagging jets. Similarly to y j 3 , this quantity helps to distinguish VBS-induced signal processes from QCD-dominated backgrounds, as much more jet activity in the rapidity region between the two tagging jets is expected for the QCD-induced production of a specific final state. Our results for the distribution of z j 3 in the CMS-tight scenario are shown in fig. 4, indicating that the third jet prefers to be close to either of the tagging jets. We observe relatively large differences between the predictions obtained with the various parton showers we considered, which reflect the differences we already observed for the rapidity of the third jet entering the definition of the z j 3 variable. This relatively large SMC uncertainty has to be kept in mind when the Zeppenfeld variable is considered as a discriminating variable between QCD-and VBS-induced processes. Up to this point, we only examined differences between predictions obtained with different parton shower algorithms, but not differences due to the settings of one and the same parton-shower program. However, it is interesting also to explore the impact of hadronization and multi-parton interactions (MPI) on predictions obtained with a specific SMC. Since our previous discussion revealed that significant differences on observables in VBS processes are to be expected mostly for distributions related to the non-tagging jets, we will restrict this comparison to those.
Interestingly, the differences due to the settings of a specific SMC are larger than one might naively expect. This is illustrated for the rapidity distribution of the third jet in fig. 5. While switching on the hadronization procedure already affects the results by an order of about 10 %, allowing for multiple parton interactions enhances the cross section even more significantly. This effect is particularly pronounced in the central-rapidity region most relevant for central-jet vetoing techniques.

Summary and conclusions
In this article we presented an implementation of VBS-induced W Zjj production in the framework of the POWHEG BOX V2.
To illustrate the capabilities of the developed code, we showed results for EW ν e e + µ − µ + jj production in realistic setups, inspired by two recently presented experimental analyses [3,4]. In particular, we investigated the impact of different parton-shower Monte-Carlo programs on experimentally accessible observables. We found the tagging jets to be rather insensitive to the specific parton-shower program used. In particular the spread in predictions using different SMCs is usually covered by the residual scale uncertainty associated with the NLO calculation. The description of subleading jets is plagued by larger uncertainties that are often outside the scale uncertainty band. The various SMCs differ widely in how they fill the volume between the two tag jets with softer emissions, as can be seen by inspecting the z j 3 distribution. We therefore strongly recommend the use of NLO+PS simulations in analyses making use of central-jet veto techniques, preferably comparing the outputs of more than one SMC, as calculations based on an LO approximation for a typical VBS signature are not capable of accurately describing the kinematics of central jets.
We also investigated the impact of hadronisation and multiple parton interactions on the third jet observables. Although these effects are more modest than the spread in SMCs, they can still amount to ∼ 10 %, a fact which should also be taken into account in central-jet veto analyses.
The computer code we developed will be made available in the public version of the POWHEG BOX V2, accessible at the webpage http://powhegbox.mib.infn.it/.

A Appendix: Results in the ATLAS-loose setup
In this appendix we show plots in the ATLAS-loose scenario. The observables shown are exactly the same as the ones shown in the main body of the paper for the CMS-tight scenario.
For the ATLAS-loose scenario, we require the presence of at least two jets with transverse momenta and rapidities of p T,j > 40 GeV , |y j | < 4.5 , respectively. The two hardest jets satisfying these criteria are called "tagging jets", and are required to exhibit an invariant mass  In addition to the tagging jets each signal event is supposed to contain an e + , a µ + , and a µ − with |y | < 2.5 , where generically denotes a charged lepton and ∆R ij the separation of two particles i, j in the rapidity-azimuthal-angle plane. The muons are required to be close in invariant mass to the Z-boson resonance with while a harder transverse-momentum cut is imposed on the electron, Figure 6 shows the transverse-momentum distributions of the hardest tagging jet and of the positively charged muon, while in fig. 7 the transverse-momentum distribution of the hardest positron and the azimuthal-angle separation of the two tagging jets are illustrated for the ATLAS-loose scenario. Analogous to the case of the CMS-tight setup, parton-shower dependencies are very small for distributions that are related to the two tagging jets.
Larger differences between the various SMCs are again observed for distributions involving a non-tagging jet, such as the transverse-momentum and rapidity distributions of the third jet, shown in fig. 8, and the Zeppenfeld variable depicted in fig. 9. We note that, in general, differences between predictions obtained with various SMCs are slightly  Figure 7: Transverse-momentum distribution of the hardest positron (left) and azimuthal-angle separation of the two tagging jets (right) in the ATLAS-loose scenario. In the lower panels the ratio of the respective distribution to the PYTHIA6 reference result is shown. In each case the blue bands indicate the scale uncertainty of the PYTHIA6 simulations, statistical uncertainties are denoted by error bars.  smaller in the ATLAS-loose than in the CMS-tight setup, an effect we attribute to the more exclusive cut setup in the latter. In fig. 10 the rapidity distribution of the third jet is illustrated for different settings in the various SMCs to explore the impact of hadronization and multi-parton interactions in our predictions. In particular the impact of MPI is smaller here than in the case of the CMS-tight scenario considered in sec. 3.