Precise predictions for same-sign W-boson scattering at the LHC

Vector-boson scattering processes are of great importance for the current run-II and future runs of the Large Hadron Collider. The presence of triple and quartic gauge couplings in the process gives access to the gauge sector of the Standard Model (SM) and possible new-physics contributions there. To test any new-physics hypothesis, sound knowledge of the SM contributions is necessary, with a precision which at least matches the experimental uncertainties of existing and forthcoming measurements. In this article we present a detailed study of the vector-boson scattering process with two positively-charged leptons and missing transverse momentum in the final state. In particular, we first carry out a systematic comparison of the various approximations that are usually performed for this kind of process against the complete calculation, at LO and NLO QCD accuracy. Such a study is performed both in the usual fiducial region used by experimental collaborations and in a more inclusive phase space, where the differences among the various approximations lead to more sizeable effects. Afterwards, we turn to predictions matched to parton showers, at LO and NLO: we show that on the one hand, the inclusion of NLO QCD corrections leads to more stable predictions, but on the other hand the details of the matching and of the parton-shower programs cause differences which are considerably larger than those observed at fixed order, even in the experimental fiducial region. We conclude with recommendations for experimental studies of vector-boson scattering processes.


Introduction
Vector-boson scattering (VBS) at a hadron collider usually refers to the interaction of massive vector bosons (W ± , Z), radiated by partons (quarks) of the incoming protons, which in turn are deflected from the beam direction and enter the volume of the particle detectors. As a consequence, the typical signature of VBS events is characterised by two energetic jets and four fermions, originating from the decay of the two vector bosons. Among the possible diagrams, the scattering  Fig. 1 Sample tree-level diagrams that contribute to the process pp → μ + ν μ e + ν e jj at order O α 6 . In addition to typical VBS contributions (left), this order also possesses s-channel contributions such as decay chain (middle) and tri-boson contributions (right) process can be mediated by a Higgs boson. The interaction of longitudinally polarised bosons is of particular interest, because the corresponding matrix elements feature unitarity cancellations that strongly depend on the actual structure of the Higgs sector of the Standard Model (SM). A detailed study of this class of processes will therefore further constrain the Higgs couplings at a very different energy scale with respect to the Higgs-boson mass, and hint at, or exclude, non-Standard Model behaviours. The VBS process involving two same-sign W bosons has the largest signal-to-background ratio of all the VBS processes at the Large Hadron Collider (LHC): evidence for it was found at the centre-of-mass energy of 8 TeV [1][2][3], and it has been recently measured at 13 TeV as well [4]. Presently, the measurements of VBS processes are limited by statistics, but the situation will change in the near future. On the theoretical side, it is thus of prime importance to provide predictions with systematic uncertainties at least comparable to the current and envisaged experimental precision [5,6].
W + W + scattering is the simplest VBS process to calculate, because the double-charge structure of the leptonic final state limits the number of partonic processes and total number of Feynman diagrams for each process. Nonetheless, it possesses all features of VBS at the LHC and is thus representative of other VBS signatures. Therefore, it is the ideal candidate for a comparative study of the different simulation tools.
In the last few years, several next-to-leading-order (NLO) computations have become available for both the VBS process [7][8][9][10][11][12][13] and its QCD-induced irreducible background process [13][14][15][16][17]. All these VBS computations rely on various approximations, typically neglecting contributions which are expected to be small in realistic experimental setups [12,18]. Recently, the complete NLO corrections to W + W + have been evaluated in Ref. [19], making it possible for the first time to study in detail the quality of the VBS approximations at NLO QCD. 1 This article starts with the definition of the VBS process in Sect. 2, before describing the approximations of the various computer codes in Sect. 3. In Sect. 4 a leading-order (LO) study of the different contributions which lead to the production of two same-sign W bosons and two jets is performed. In the same section predictions for VBS from different tools are compared at the level of the cross section and differential distributions. The comparison is extended to the NLO corrections to VBS in Sect. 5. The effect of the inclusion of matching LO and NLO computations to parton shower (PS) is discussed in Sect. 6. Finally, Sect. 7 contains a summary of the article and concluding remarks.

Definition of the process
The scattering of two positively-charged W bosons with their subsequent decay into different-flavour leptons can proceed at the LHC through the partonic process: pp → μ + ν μ e + ν e jj + X. (1) At LO, this process can proceed via three different coupling-order combinations: O α 6 , O α 2 s α 4 , and O α s α 5 . The first, commonly referred to as EW contribution or VBS, 2 receives contributions from Feynman diagrams such as those depicted in Fig. 1: in addition to genuine VBS contributions (left diagram), it also features s-channel contributions with non-resonant vector bosons (centre diagram) or from triple-boson production (right diagram). Note that s-, t-, and u-channel contributions are defined according to the quark lines. The s-channel denotes all Feynman diagrams where the two initial-state partons are connected by a continuous fermion line, while for the t-and u-channel the fermion lines connect initial state quarks to final state quarks. The u-channel refers to contributions with crossed fermion lines with respect to t-channel, which appears for identical (anti-) quarks in the final state. The s-channel contributions play a particular role in the study of the various contributions in Sect. 4.1. When using approximations, care must be taken that only gauge-invariant subsets are considered to obtain physically meaningful results. We discuss the commonlyused possible choices in detail in the next section.
The second coupling combination of order O α 2 s α 4 corresponds to diagrams with a gluon connecting the two quark lines, and with the W bosons radiated off the quark lines. Because of the different colour structure, this contribution features a different kinematic behaviour than VBS. Nonetheless it shares the same final state, and therefore constitutes an irreducible background to the EW process.
Finally, the third contribution of order O α s α 5 is the interference of the two types of amplitudes described above. It is non-zero only for those partonic sub-processes which involve identical quarks or anti-quarks. Such a contribution is usually small (3%) within typical experimental cuts [19].
In the rest of this article, the notations LO or NLO(-QCD) without any specification of coupling powers refer to the contributions at order O α 6 and O α s α 6 , respectively.
In experimental measurements, special cuts, called VBS cuts, are designed to enhance the EW contribution over the QCD one and to suppress the interference. These cuts are based on the different kinematical behaviour of the contributions. The EW contribution is characterised by two jets with large rapidities as well as a large di-jet invariant mass. The two W bosons are mostly produced centrally. This is in contrast to the QCD contribution which favours jets in the central region. Therefore, the event selection usually involves rapidity-difference and invariant-mass cuts for the jets. Note that, as pointed out in Ref. [19], when considering full amplitudes, the separation between EW and QCD production becomes ill-defined. Hence, combined measurements which are theoretically better defined should be preferably performed by the experimental collaborations at the LHC.

Theoretical predictions for VBS
We now discuss the various approximations which are implemented in computer programs for the EW contribution at order O α 6 . Since we are mostly interested in the scattering of two W bosons, which includes the quartic gauge-boson vertex, it may appear justified to approximate the full process by considering just those diagrams which contain the 2 → 2 scattering process as a sub-part. However, this set of contributions is not gauge invariant. In order to ensure gauge invariance, an on-shell projection of the incoming and outgoing W bosons should be performed. While this can be done in the usual way for the time-like outgoing W bosons, the treatment of the space-like W bosons emitted from the incoming quarks requires some care. Following Refs. [22,23] these W-boson lines can be split, the W bosons entering the scattering process can be projected on-shell, and the emission of the W bosons from the quarks can be described by vectorboson luminosities. Such an approximation is usually called effective vector-boson approximation (EVBA) [24][25][26].
An improvement of such an approximation consists in considering all t-and u-channel diagrams and squaring them separately, neglecting interference contributions between the two classes. These interferences are expected to be small in the VBS fiducial region, as they are both phase-space and colour suppressed [12,18]. The s-channel squared diagrams and any interferences between them and the t/u-channels are also discarded. This approximation is often called t-/uapproximation, VBF, or even VBS approximation. We adopt the latter denomination in the following. This approximation is gauge-invariant, a fact that can be appreciated by considering the two incoming quarks as belonging to two different copies of the SU(3) gauge group.
A further refinement is to add the squared matrix element of the s-channel contributions to the VBS approximation.
The approximations performed at LO can be extended when NLO QCD corrections to the order O α 6 are computed. The VBS approximation can be extended at NLO in a straightforward manner for what concerns the virtual contributions. For the real-emission contributions special care must be taken for the gluon-initiated processes. The initial-state gluon and initial-state quark must not couple together, otherwise infrared (IR) divergences proportional to s-channels appear, which do not match with the ones found in the virtual contributions. The subset of diagrams where all couplings of the initial state gluon to initial state quark are neglected forms a gauge-invariant subset, with the same argument presented above.
A further refinement is to consider the full real contributions, which include all interferences, and part of the virtual. In particular one can consider only one-loop amplitudes where there is no gluon exchange between the two quark lines and assuming a cancellation of the IR poles.
When considering the full NLO corrections of order O α s α 6 , besides real and virtual QCD corrections to the EW tree-level contribution of order O α 6 also real and virtual EW corrections to the LO interference of order O α s α 5 have to be taken into account. Since some loop diagrams contribute to both types of corrections, QCD and EW corrections cannot be separated any more on the basis of Feynman diagrams, and the cancellation of IR singularities requires the inclusion of all of them [19].

Description of the programs used
In the following, the codes employed throughout this article and the approximations implemented in each of them are discussed: -Phantom [27] is a dedicated tree-level Monte Carlo for six-parton final states at pp, pp, and e + e − colliders at orders O(α 6 ) and O(α 2 s α 4 ) including interferences between the two sets of diagrams. It employs complete tree-level matrix elements in the complex-mass scheme [28][29][30] computed via the modular helicity formalism [31,32]. The integration uses a multi-channel approach [33] and an adaptive strategy [34]. Phantom generates unweighted events at parton level for both the SM and a few instances of beyond the Standard Model (BSM) theories.
-WHIZARD [35,36] is a multi-purpose event generator with LO matrix-element generator O'Mega. For QCD amplitudes it uses the colour flow formalism [37].
For NLO QCD calculations [38], where Whizard is in the final validation phase, it provides FKS subtraction terms [39,40], while virtual matrix elements are provided externally by OpenLoops [41] or Recola [42,43]. Furthermore, WHIZARD can automatically provide POWHEG matching to parton shower [44]. WHIZARD allows to simulate a huge number of BSM models as well, in particular in terms of higher-dimensional operators for VBS processes including means to provide unitarity limits [45,46]. -The program Bonsay [47] consists of a general-purpose Monte Carlo integrator and matrix elements taken from different sources: Born matrix elements are adapted from the program Lusifer [48], which have been generalised to calculate also real matrix elements. Virtual matrix elements have been calculated using an in-house matrixelement generator. One-loop integrals are evaluated using the Collier library [49,50] Same-sign W + W + j j production is provided via the pro-cess ID 250. The corresponding s-channel contributions are available separately as tri-boson processes with semileptonic decays, with process IDs between 401 and 492. For the final state studied in this article, only W + W + W − production with a hadronically decaying W − boson, process ID 432, can contribute. These can simply be added on top of the VBS contribution. Interferences between the two are therefore neglected. The usage of leptonic tensors in the calculation, pioneered in Ref. [7], thereby leads to a significant speed improvement over automatically generated code. Besides the SM, also a variety of new-physics models including anomalous couplings of the Higgs and gauge bosons can be simulated.
MC) is an automatic meta-code (a code that generates codes) which makes it possible to simulate any scattering process including NLO QCD corrections both at fixed order and including matching to parton showers, using the MC@NLO method [59]. It makes use of the subtraction method by Frixione, Kunszt and Signer (FKS) [39,40] (automated in the module MadFKS [60,61]) for regulating IR singularities. The computations of one-loop amplitudes are carried out by switching dynamically between two integral-reduction techniques, OPP [62] or Laurent-series expansion [63], and tensor-integral reduction [64][65][66]. These have been automated in the module MadLoop [67], which in turn exploits CutTools [68], Ninja [69,70], IREGI [71], or Collier [50], together with an in-house implementation of the OpenLoops optimisation [41]. Finally, scale and PDF uncertainties can be obtained in an exact manner via reweighting at negligible additional CPU cost [72]. The simulation of VBS at NLO-QCD accuracy can be performed by issuing the following commands in the program interface: set complex_mass_scheme import model loop_qcd_qed_sm_Gmu generate p p > e+ ve mu+ vm j j QCD=0 [QCD] output With these commands the complex-mass scheme is turned on, then the NLO-capable model is loaded, 3 finally the process code is generated (note the QCD=0 syntax to select the purely-EW process) and written to disk. No approximation is performed for the Born and realemission matrix elements. Only strongly-interacting particles circulating in the loops are generated for the virtual matrix element. The version capable of computing both QCD and EW corrections will be released in  [42,43]. It can compute arbitrary processes at the LHC with both NLO QCD and EW accuracy in the SM. This is made possible by the fact that Recola can compute arbitrary processes at tree and one-loop level in the SM. To that end, it relies on the Collier library [49,50] to numerically evaluate the one-loop scalar and tensor integrals. In addition, the subtraction of the IR divergences appearing in the real corrections has been automatised thanks to the Catani-Seymour dipole formalism for both QCD and QED [73,74]. The code has demonstrated its ability to compute NLO corrections for high-multiplicity processes up to 2 → 7 [75,76]. In particular the full NLO corrections to VBS and its irreducible background [19,77] have been obtained thanks to this tool. One key aspect for these high-multiplicity processes is the fast integration which is ensured by using similar phase-space mappings to those of Refs. [28,48,78]. In MoCaNLO+Recola no approximation is performed neither at LO nor at NLO. It implies that, also contributions stemming from EW corrections to the interference are computed.
We conclude this section by summarising the characteristics of the various codes in Table 1. In particular, it is specified whether

Input parameters
The hadronic scattering processes are simulated at the LHC with a centre-of-mass energy √ s = 13 TeV. The NNPDF 3.0 parton distribution functions (PDFs) [79] with five massless flavours, 4 NLO-QCD evolution, and a strong coupling constant α s (M Z ) = 0.118 5 are employed. 6 Initial-state collinear singularities are factorised according to the MS scheme, consistently with what is done in NNPDF.
For the massive particles, the following masses and decay widths are used: The measured on-shell (OS) values for the masses and widths of the W and Z bosons are converted into pole values for the gauge bosons (V = W, Z) according to Ref. [81], The EW coupling is renormalised in the G μ scheme [51] according to with and where M 2 V corresponds to the real part of the squared pole mass. The numerical value of α, corresponding to the choice of input parameters is The Cabibbo-Kobayashi-Maskawa matrix is assumed to be diagonal, meaning that the mixing between different quark generations is neglected. The complex-mass scheme [28][29][30] is used throughout to treat unstable intermediate particles in a gauge-invariant manner. The central value of the renormalisation and factorisation scales is set to defined via the transverse momenta of the two hardest jets (identified with the procedure outlined in the following), event by event. 7 This choice of scale has been shown to provide stable NLO-QCD predictions [12].
Following experimental measurements [1][2][3]82], the event selection used in the present study is: -The two same-sign charged leptons are required to fulfil cuts on transverse momentum, rapidity, and separation in the rapidity-azimuthal-angle separation, -The total missing transverse momentum, computed from the vectorial sum of the transverse momenta of the two neutrinos, is required to be p T,miss > 40 GeV.
7 By default, the renormalisation and factorisation scales employed in the Powheg-Box slightly differ from the ones defined in Eq. (7), as the momenta of two final-state quarks in the underlying Born event are employed instead of those of the two hardest jets.
-QCD partons (light quarks and gluons) are clustered together using the anti-k T algorithm [83], possibly using the FastJet implementation [84], with distance parameter R = 0.4. We impose cuts on the jets' transverse momenta, rapidities, and their separation from leptons, VBS cuts are applied on the two jets with largest transverse momentum, unless otherwise stated. In particular, we impose a cut on the invariant mass of the di-jet system, as well as on the rapidity separation of the two jets, if not explicitly stated otherwise. -When EW corrections are computed, real photons and charged fermions are clustered together using the anti-k T algorithm with radius parameter R = 0.1. In this case, leptons and quarks are understood as dressed fermions.

Contributions
In the present section, the cross sections and distributions are obtained without applying the VBS cuts on the variables m jj and |Δy jj |, Eq. (11). In Table 2, the cross sections of the three LO contributions are reported. The EW, QCD, and interference contributions amount to 57%, 37%, and 6% of the total inclusive cross section, respectively. The QCD contribution does not possess external gluons due to charge conservation. Thus the diagrams of order O(g 2 s g 4 ) only involve gluon exchange between the quark lines. This results in a small contribution even if the VBS cuts have not been imposed. The interference between EW and QCD contributions is small, due to colour suppression, but not negligible.
In Fig. 2, these three contributions are shown separately and summed in the differential distributions in the di-jet invariant mass m jj and the rapidity difference |Δy jj |. For the di-jet invariant-mass distribution (left), one can observe that the EW contribution peaks around an invariant mass of about Table 2 Cross sections at LO accuracy for the three contributions to the process pp → μ + ν μ e + ν e jj, obtained with exact matrix elements. These results are for the set-up described in Sect. 3.3 but no cuts on m jj and |Δy jj | are applied. The uncertainties shown refer to the estimated statistical errors of the Monte Carlo integrations GeV. This is due to diagrams where the two jets originate from the decay of a W boson (see middle and right diagrams in Fig. 1). Note that these contributions are not present in calculations relying on the VBS approximation as these are only s-channel contributions. The EW contribution becomes dominant for di-jet invariant mass larger than 500 GeV. The same holds true for jet-rapidity difference larger than 2.5 (right). This justifies why cuts on these two observables are used in order to enhance the EW contribution over the QCD one. In particular, in order to have a large EW contribution, rather exclusive cuts are required.
This can also be seen in Fig. 3 where the three contributions are displayed as double-differential distributions in the di-jet invariant mass and jet rapidity difference. Again, it is clear that the region with low di-jet invariant mass should be avoided in VBS studies as it is dominated by tri-boson contributions. This motivates in particular the choice of the cut m jj > 200 GeV for our LO inclusive study in Sect. 4.2.

Inclusive comparison
In Fig. 4, ratios for double-differential cross sections in the variables m jj and |Δy jj | are shown. 8 Two plots are displayed: the ratios of the |t| 2 +|u| 2 and |s| 2 +|t| 2 +|u| 2 approximations over the full calculation. In the first case, the approximation is good within ± 10% over the whole range apart from the low invariant-mass region at both low and large rapidity difference. The low rapidity-difference region possesses remnants of the tri-bosons contribution that have a di-jet invariant mass 8 In Fig. 4, the level of the accuracy of the predictions in each bin is around a per mille. around the W-boson mass. It is therefore expected that the |t| 2 + |u| 2 approximation fails in this region. The second plot, where the |s| 2 + |t| 2 + |u| 2 approximation is considered, displays a better behaviour in the previously mentioned region. The full calculation is approximated at the level of ± 5% everywhere apart from the region where |Δy jj | < 2.

Comparison in the fiducial region
In Table 3, we report the total rates at LO accuracy at order O(α 6 ) obtained in the fiducial region described in Sect. 3.3. Two things should be highlighted here: first, despite the different underlying approximations, the two most-distant predictions (Powhegbox and MG5_aMC) are only 0.7% apart. This simply means that the details of the various VBS approximations have an impact below 1% at the level of the fiducial cross section at LO for a typical phase-space volume used by experimental collaborations. This is in agreement with the findings of Refs. [12,18]. Second, the four complete predictions (Whizard, Phantom, MG5_aMC, and MoCaNLO+Recola) are not in statistical agreement. While we have checked the point-wise agreement of the matrixelement, we cannot exclude other reasons for the disagreement, for example a non-representative (i.e. too-aggressive) estimate of the Monte Carlo uncertainty or a non-perfect mapping of the six-body phase-space. However, the level of ambiguity is at the 0.5% level, which we deem satisfactory compared to the larger differences observed at NLO or when including matching to parton shower.
In Fig. 5, we show the distributions in the invariant mass (left) and the rapidity difference (right) of the two tagging jets which are key observables for VBS measurements. In both  Fig. 4 Ratios for double-differential distributions in the variables m jj and |Δy jj | at LO i.e. order O(α 6 ) of approximated squared amplitudes over the full matrix element. The approximated squared amplitudes are computed as |A| 2 ∼ |t| 2 + |u| 2 (left) and |A| 2 ∼ |s| 2 + |t| 2 + |u| 2 (right). The cuts applied are the ones of Sect. 3.3 and no cuts on m jj and |Δy jj | are applied We compare three different predictions at NLO QCD: the VBS approximation implemented in Bonsay (dubbed |t| 2 + |u| 2 ), the VBS approximation with the s-channel contributions from VBFNLO (dubbed |s| 2 + |t| 2 + |u| 2 ), and the full computation. The full computation employs exact matrix elements meaning that t/u/s interferences, factorisable and  non-factorisable QCD corrections, as well as EW corrections to the order O(α s α 5 ) are included. The total cross sections within the above-mentioned kinematic cuts are shown in Table 4. The |t| 2 +|u| 2 approximation for NLO QCD predictions is lower by about 6% than the full calculation. The inclusion of s-channel diagrams improves the approximate prediction, leading to an excess at the 3% level.
These differences are more evident in differential distributions. In Fig. 6, we show the differential distributions in the di-jet invariant mass m jj and rapidity separation |Δy jj |. For large m jj and large |Δy jj |, as expected, the VBS approximation is performing well and its s-channel extension agrees with the full calculation within 10%. This is in contrast with the regions 200 GeV < m jj < 500 GeV and 2 < |Δy jj | < 2.5, where the difference between the |t| 2 + |u| 2 approximation and the full computation can be above 30%. The inclusion of s-channel contributions cures partly this behaviour by improving the approximation to about 10%. This tends to indicate that interference contributions and/or non-factorisable QCD corrections play a non-negligible role in this phase-space region.
In order to investigate further the jet-pair kinematics, we study the double-differential distribution in the variables m jj and |Δy jj |. In particular, in Fig. 7, we compute in each bin the ratios of the approximated cross sections over the full ones [σ (|t| 2 + |u| 2 )/σ (full) and σ (|s| 2 + |t| 2 + |u| 2 )/σ (full)]. As expected, in the low invariant-mass and low rapidityseparation region of the jet pair (200 GeV < m jj < 500 GeV, 2 < |Δy jj | < 2.5) the VBS approximation fails significantly (by more than 40%). Including the s-channel contributions leads to a difference of less than 10% in this very region. However, in the region of large di-jet invariant mass and low rapidity separation of the jets, the |s| 2 + |t| 2 + |u| 2 approximation overestimates the full computation by more than 40%. 9 Again, this seems to support the fact that interferences and non-factorisable corrections can be non-negligible in this region. On the other hand, in the typical VBS region, the VBS approximation shows a good agreement with the full computation as documented in detail in Sect. 5.2.
In Fig. 8, the distributions in the transverse momentum of the hardest jet and its rapidity are shown. At low transverse momentum, |t| 2 + |u| 2 and |s| 2 + |t| 2 + |u| 2 approximations are lower and higher than the full computation by about 20%, respectively. At high transverse momentum, they have a similar behaviour. They both diverge from the full computation towards larger transverse momentum (about 10% at 1000 GeV). Regarding the rapidity of the hardest jet, the two approximations have opposite behaviours. In the central region, the |t| 2 + |u| 2 approximation differs by 12% with respect to the full computation, while the |s| 2 + |t| 2 + |u| 2 one is good within 5%. In the peripheral region, the |t| 2 +|u| 2 approximation is rather close to the full computation (5%), while the |s| 2 + |t| 2 + |u| 2 one differs by 10%.
Concerning leptonic observables, we show in Fig. 9 the distributions in the di-lepton invariant mass and in the Zeppenfeld variable of the electron, defined as Analogous definitions are later also used for the Zeppenfeld variable of the muon and of the third jet. The |s| 2 +|t| 2 +|u| 2 predictions for m e + μ + agree rather well with the full curve, obtained from MoCaNLO+Recola. The prediction from Bonsay is about 10% lower around 1000 GeV. The Zeppenfeld variable of the positron z e is more strongly affected by the exclusion of s-channel contributions. For increasing z e , the |t| 2 +|u| 2 approximation diverges from the full computation to reach a difference of about 25% at 1.5. On the other hand, including s-channel contributions leads to a better approximation, staying within 10% difference over the whole range.
In conclusion, both the loose minimum di-jet invariantmass cut and the inclusion of QCD radiative corrections render the s-channel contributions less suppressed than at LO, making their inclusion mandatory, in order to provide trustworthy predictions at NLO accuracy. In the inclusive region studied here, neglecting s-channel contributions, nonfactorisable corrections, and EW corrections can lead to discrepancies of up to 30% with respect to the full computation. Nevertheless, the VBS approximation at NLO provides a good approximation of full calculations in the kinematic region where m jj 500 GeV and |Δy jj | 2.5), for both total cross section and differential distributions. This more exclusive region is studied in more detail in the next section.

Comparison in the fiducial region
In Table 5, the cross sections of the various tools at NLO-QCD accuracy are presented. The order considered is again  Fig. 7 Ratios for double-differential distributions in the variables m jj and |Δy jj | at NLO QCD i.e. at order O(α s α 6 ) of the approximated squared amplitudes over the full matrix element. The approximated squared amplitudes are computed as |A| 2 ∼ |t| 2 + |u| 2 (left) and |A| 2 ∼ |s| 2 + |t| 2 + |u| 2 (right). In addition to the cuts of Sect. 3.3, the VBS cuts take the values m jj > 200 GeV and |Δy jj | > 2 the order O(α s α 6 ), and the fiducial volume is the one described in Sect. 3.3. In contrast with Table 3, the NLO predictions differ visibly according to the approximations used.
The first observation is that the predictions featuring two versions of the VBS approximation (Bonsay and the Powheg-Box) are relatively close. 10 Bonsay uses a double-pole approximation for the virtual matrix element, and it is worth noticing that this approximation seems to  approximates the full computation by a per cent. The main contribution due to the s-channel diagrams thereby consists of real-emission contributions, where one of the two leading jets is formed by one quark, or possibly also both quarks, originating from the W-boson decay, and the second one by the extra radiation emitted from the initial state. In such configurations, the hadronically-decaying W boson can become on-shell and hence yield larger contributions than at LO, where the invariant mass cut on the two jets forces the boson into the far off-shell region. However, the agreement between MoCaNLO+Recola and VBFNLO is mostly accidental, as the inclusion of interference effects and some non-factorisable corrections (in the real corrections) in MG5_aMC brings the prediction down and closer to the VBS approximation. Not unexpectedly none of the approximations used here agrees perfectly with the full calculation of MoCaNLO+Recola at NLO. Nevertheless, the disagreement seems never to exceed 2% at the fiducial cross-section level. In Figs. 10, 11 and 12, several differential distributions are shown. All these predictions are performed at NLO accuracy at the order O(α s α 6 ). In the upper panel, the absolute predictions are shown while in the lower panel, the ratio with respect to the full predictions are displayed. The band corresponds to a seven-points variation of the factorisation and renormalisation scales (as defined in Eq. (3.11) of Ref. [19]).
We start with Fig. 10 which displays the invariant mass (left) and the rapidity separation (right) of the two tagging jets. For high invariant mass, all predictions agree rather well. On the other hand, for low invariant mass, the hierarchy present at the level of the cross section is reproduced. The VBS-approximated predictions (Bonsay and Powheg-Box) are lower than the full calculation (MoCaNLO+Recola). The full calculation is rather well approximated by the hybrid VBS approximation implemented in MG5_aMC. Finally, VBFNLO which includes also s-channel contributions provides larger predictions at low invariant mass. For the rapidity difference between the two tagging jets, the hierarchy between the predictions is rather similar. Therefore, depending on the approximation used, it can vary by ± 7% and ± 4% with respect to the full computation at low invariant mass and low rapidity difference for the tagging jets, respectively Concerning the transverse momentum (left) and rapidity (right) of the hardest jet shown in Fig. 11, the situation is rather different. While MG5_aMC is very close to the full prediction for low transverse momentum, it departs from it at larger transverse momentum by about 10%. This is in contrast with the VBS-approximated predictions such as Bonsay, Powheg, and VBFNLO which are lower than the full computation at low transverse momentum and higher for larger transverse momentum. The difference at high transverse momentum between the latter predictions and the full computation can be attributed to EW Sudakov logarithms that become large in this phase-space region. While the predictions of Bonsay and Powheg are rather close over the whole range, the one of VBFNLO is very different at low transverse momentum where it is even higher than the full computation. We note that for the transverse momentum of the second hardest jet, the predictions from MG5_aMC are in good agreement with the other VBS-approximated predictions. Concerning the rapidity of the hardest jet, VBFNLO is in good agreement with MoCaNLO+Recola in the rapidity range |y j 1 | < 3. For larger rapidity, the other codes constitute a better description of the full process at order O(α s α 6 ).
The last set of differential distributions is the invariant mass of the two charged leptons (left) and the Zeppenfeld variable for the anti-muon (right). Concerning the comparison of the predictions, both distributions display a rather similar behaviour. Indeed, the hierarchy mentioned previously is here respected and enhanced towards high invariant mass or high Zeppenfeld variable. The predictions of MoCaNLO+Recola and VBFNLO are in rather good agreement for both distributions for the kinematic range displayed here. The other three VBS approximations are close to each other within few per cent.
In the end, the quality of the VBS approximations is good up to 10% in the fiducial region. These differences are larger than those at LO.
The contributions from the s-channel amplitude can be sizeable especially at low invariant mass for the two tagging jets (comparing the predictions of VBFNLO against the ones of Bonsay and Powheg). This can be explained by the fact that s-channel contributions are less suppressed at NLO. As real radiation, an extra gluon-jet can be radiated from any of the strongly-interacting particles while the two quarks originating from the W-boson decay can be recombined in a single jet. Therefore, the jet requirements (m jj > 500 GeV and |Δy jj | > 2.5) that were suppressing s-channel contributions at LO are partially lifted with the inclusion of a third jet at NLO. Such an effect has also been observed for top-antitop production in the lepton+jet channel at NLO QCD [85].
In phase-space regions where the s-channel contributions are sizeable their interference with the t/u-channel can be of similar size. This can be observed by comparing the predictions of VBFNLO against the ones of MG5_aMC.
Finally, the effect of EW corrections and non-factorisable contributions in the virtual corrections are usually small. But they can be relatively large (about 10%) for large transverse momentum of the hardest jet. These high-energy region of the phase space are where EW Sudakov logarithms become large. Nonetheless these regions are rather suppressed and thus these effects are hardly visible at the level of the cross section.

Matching to parton shower
We now discuss how different predictions compare when the matching to parton shower is included, both at LO (i.e. at order O(α 6 )) and at NLO-QCD (i.e. at order O(α 6 α s )) accuracy. For such a comparison we expect larger discrepancies than what we found at fixed order, as a consequence of the different matching schemes, parton showers employed, and of other details of the matching (such as the choice of the parton shower initial scale). Among the codes capable of providing fixed-order results, presented before, MG5_aMC, the Powheg-Box, and VBFNLO can also provide results at (N)LO+PS accuracy. For VBFNLO matched to Herwig and the Powheg-Box, we restrict ourselves to show results only in the VBS approximation, i.e. the s-channel contributions are neglected here. Besides, also Phantom and Whizard are used for LO+PS results.  1.2). For the latter, the default angular-ordered shower is employed. The same parton showers are employed for the LO results of Phantom. Pythia8 is also employed for the LO results of Whizard. For the Powheg-Box, the namesake matching procedure is employed [53,54], together with Pythia8 (version 8.230). VBFNLO serves as a matrix-element and phase-space provider for the Matchbox module [89] of Herwig7 [87,88], using an extended version of the Binoth Les Houches Accord interface [90][91][92]. The Matchbox module makes it possible to choose between MC@NLOlike and Powheg-like matching. As parton shower, both the default angular-ordered shower as well as the dipole shower can be employed. We use here the subtractive (MC@NLOtype) matching to these parton-shower algorithms. Whenever Pythia8 is used, the Monash tune [93] is selected. Multipleparton interactions are disabled.
Results are presented within the cuts described in Sect. 3.3, applied after shower and hadronisation (this implies that jets are obtained by clustering stable hadrons, and not QCD partons). It follows that at the event-generation level, looser cuts (or even no cuts at all) must be employed in order not to bias the results. This also implies that the tagging jets, whose momenta enter in the renormalisation and factorisation scales, Eq. (7), are now defined without imposing the ΔR j cut. The effect of this change is below one per cent at the level of the fiducial cross sections at NLO.
Compared to the fixed-order computations, a slightly different set-up has been employed for MG5_aMC in order to simplify the calculation: instead of generating the full pp → μ + ν μ e + ν e jj process, since it is dominated by doublyresonant contribution, the events are produced for the process with two stable W + -bosons (pp → W + W + jj), and the decay of these W + -bosons is simulated with MadSpin [94] (ensuring spin correlations) before the parton shower. Since MadSpin computes the partial and total decay widths of the W bosons at LO accuracy only, while in Sect. 3.3 the NLO width is employed, an effect (6%) on the normalisation is induced.
We now present the results of predictions matched to parton showers. The total rates within VBS cuts are displayed in Tables 6 and 7, at LO and NLO accuracy respectively. For MG5_aMC, the numbers with Γ resc are rescaled to take into account the width effects described in the above paragraph. At NLO accuracy, for MG5_aMC + Pythia8 and VBFNLO+Herwig7-Dipole, we also quote theoretical uncertainties. For the former, we show both PDF and   11 obtained via exact reweighting [72] by varying independently the renormalisation and factorisation scales by a factor of two around the central value, Eq. (7) (nine-point variations). For the latter, we show the threepoint scale uncertainties, obtained by considering correlated variations of the renormalisation, factorisation, and shower starting scale. Theory uncertainties should have very little dependence on the tool employed. We observe that, once the width effect is taken into account, total rates from different tools agree within some per cents (≤ 7%), both at LO and NLO. Larger discrepancies, however, appear for differential observables, which we discuss in the following. Theory uncertainties on the total rates are very small, regardless of whether scale variations are estimated with independent or correlated variations of the renormalisation and factorisation scales. Concerning differential distributions, for each observable we display results in two plots, shown side-by-side. In the plot on the left (right), (N)LO+PS predictions are shown with different colours in the main frame. In the inset, these predictions are compared in both cases with a fixed-order prediction at NLO accuracy (obtained with VBFNLO, i.e. the VBS approximation with s-channel). For the differential observables, the MG5_aMC predictions are not rescaled to compensate for the width effect mentioned above. As for the table, we show theoretical uncertainties for the NLO+PS samples obtained with VBFNLO and MG5_aMC: again, for the first the band corresponds to three-point variations, while for the second the darker (lighter) band corresponds to nine-point scale variations (plus PDF uncertainties, linearly added). The first observable we investigate is the exclusive jet multiplicity, shown in Fig. 13. Looking at the LO+PS predic- tions, one can appreciate that the main effects are driven by the parton shower that is employed (Herwig7 or Pythia8), with the clear tendency of producing more radiation for the latter, leading to higher jet multiplicities. Differences among tools that employ the same parton shower are typically smaller, and can be traced back to different values of the initial scale of the parton shower (the scalup entry of the Les Houches Event (LHE) file [95,96]). This event-byevent number corresponds to the maximum hardness (translated into the shower-evolution variable) of the radiation that can be generated by the shower. 12 The main effect of NLO corrections for this (rather inclusive) observable is to stabilise the predictions for the two-jet bin, where discrepancies among tools are reduced to about 10%. For the three-jet bin, which is described only at LO accuracy, differences among tools remain large, and are possibly related to the underlying approximation performed (MG5_aMC is the only tool to use the full matrix element for the real radiation), in particular the inclusion of the s-channel contributions: the largest rate is predicted by MG5_aMC, while the smallest rate is predicted by the Powheg-Box, both matched to Pythia8. Despite the fact that the same parton shower is employed, the way emissions are treated is different among the two tools. In particular, for the Powheg-Box, the first emission is generated with an internal Sudakov form factor (the prediction dubbed Powheg-LHE corresponds to stopping after the first emission), while for MG5_aMC there is an interplay between the real-emission matrix element and the shower emission. For this observable we also show the prediction obtained with MG5_aMC+Pythia8 by reducing the shower starting scale by a factor 2 with respect to the default value, dubbed MG5_aMC+Py8, Q sh /2. 13 The main effect of reducing the shower scale is that events migrate from the three-jet bin into the two-jet bin, i.e. less radiation is generated. The size of this effect on the jet rate is + 4% (− 8%) on the two (three) jet bin, while the total rate within cuts is left unchanged.
The next observable that we study is the invariant mass of the two tagging jets, shown in Fig. 14. For this observable, both at LO+PS and NLO+PS, the spread of predictions matched with parton shower is rather small ( 10%, if one compensates for the 6% width effect for MG5_aMC). The LO+PS predictions tend to be significantly softer than the fixed NLO one, with an effect of about − 30% at the end of the displayed range. At NLO+PS, this effect is mitigated, owing to the better description of the first QCD emission which is now driven by the real-emission matrix element. For this observable (and all the others which are NLO accurate) the effect of reducing the shower scale is negligible, hence it is not shown.
The rapidity difference between the two tagging jets, shown in Fig. 15, has some similarities with the invariantmass distribution. At LO+PS all predictions show the tendency to deplete the large-separation region with respect to the fixed-order prediction, in a quantitatively similar way, except for VBFNLO+Herwig7 where the effect is mitigated. At NLO+PS, when the extra radiation is described by the real matrix element, such an effect is greatly reduced. A notable exception is the Powheg-Box prediction, which still shows a suppression at large separations. Since such a suppression is already there for the Powheg-LHE sample, it is very likely that it is driven by the way the first emission is generated. A minor effect in the same direction is visible in the last two bins of the MG5_aMC+Herwig7 prediction (although with rather large statistical uncertainties).
The transverse momentum of the hardest and secondhardest jets are shown in Figs. 16 and 17, respectively. In general, for both observables, predictions from different tools agree rather well with each other, with a spread at most at the 10% level. At LO+PS, typically the transverse-momentum spectra are softer than the fixed-NLO ones, and this effect is more marked for the second-hardest jet which, as expected, is more sensitive to the description of the extra radiation.
Again, this effect is mitigated by NLO corrections. The only feature that may be worth noticing among the NLO+PS predictions is the tendency of the Powheg-Box to suppress the hardest-jet spectrum at low transverse momentum ( p T,j 1 < 100 GeV).
If we consider the rapidity of the second jet, Fig. 18, we observe again rather small differences among tools, with the tendency towards a general stabilisation at NLO+PS. However, some (small) differences in the shape remain at NLO+PS, which are worth to be briefly discussed: predictions obtained with MG5_aMC are very close to the fixedorder prediction; the Powheg-Box displays an enhancement of the central region, and a consequent suppression in the peripheral region, while VBFNLO shows an opposite behaviour. However, the effect is rather small, with the largest departure from the fixed-order prediction being at most 10%. 14 Finally, focusing on the third jet, we conclude the list of differential observables by showing the Zeppenfeld variable defined in Eq. (12), Fig. 19. This variable is closely related to the third jet rapidity, and small (large) values of z correspond to central (peripheral) rapidities. In general, for observables which involve the third jet, one can clearly see a degradation of the agreement among the various 14 If the setting SpaceShower:dipoleRecoil = on (discussed in the following) is used when Pythia8 is employed together with the Powheg-Box, the enhancement at central rapidities and the depletion at small value of transverse momentum are partially compensated. tools, because of the poorer perturbative description of these observables. The Zeppenfeld variable is a striking example: both at LO and NLO, the tendency of Pythia8 to generate more hard and central radiation, corresponding to low val-ues of z, is clearly visible. Such an effect, which is related to the way Pythia8 deals with the recoil of the radiation in VBF(VBS)-type processes, can be mitigated by setting SpaceShower:dipoleRecoil = on in the Pythia8 input file. 15 It is interesting to notice that the effect survives beyond the first emission, as it can be observed by prefer the description by Herwig++ [98,99]. The central enhancement is a bit mitigated if NLO+PS tools are used (compare LO+PS and NLO+PS from MG5_aMC+Pythia8 with the fixed-NLO prediction), however even at NLO+PS the central region (z j 3 < 0.5) is cursed by huge differences between tools. Large differences, reaching a factor 2, persist also away from the central region. These findings are consistent with behaviour displayed in Refs. [11,[100][101][102][103] where the behaviour of NLO matching in VBS processes has been reported.
In conclusion, the comparison of tools including matching to parton shower clearly shows the benefits of the inclusion of NLO corrections: for most observables described effectively at NLO accuracy differences between tools are at (or below) the 10% level. Some exceptions exist, e.g. the rapidity separation of the two tagging jets, which on the one hand clearly suggest not to rely on a single tool/parton shower, and on the other make it worth investigating more in detail the way QCD radiation is generated, e.g. when fully-differential computations at NNLO will become available (for VBF Higgs production, see Refs. [104,105]). It is a remarkable fact that, even for those observables that display small discrepancies, the theoretical uncertainty obtained via scale variations (renormalisation, factorisation, and shower scale) systematically underestimates the spread of predictions. We note that in the only VBF process where NNLO corrections are known, i.e. VBF Higgs production [104,105], the NLO scale-uncertainty band does not include the NNLO prediction. This suggests that the NLO scale variation underestimates the size of the perturbative uncertainty. Again, this stresses the need to employ at least two different tools in order to obtain a more realistic estimate of theoretical uncertainties. Finally, the size of discrepancies for observables that are described at a lower perturbative accuracy, notably those related to the third jet, suggests that experimental analyses should rely as little as possible on those observables and, in any case, use conservative estimates of the theory uncertainties. On the one hand, in order to improve the description of these observables, a simulation of VBS+j at NLO accuracy, currently unavailable but within the reach of modern automated tools, is certainly desirable. On the other hand, measurements of processes with similar colour flow (EW production of a single vector boson plus jets, VBF, . . .) can certainly help in order to discriminate which tools perform better in the comparison with data [97,106].

Conclusions and recommendations
In the present article, a detailed study of the process pp → μ + ν μ e + ν e jj + X at the LHC has been presented, mainly focused on the EW production mechanism which involves the scattering of massive vector bosons. Until very recently, when the complete calculation became available for the NLO QCD corrections (order O(α s α 6 )), the so-called VBS approximation was the standard for this kind of simulations. For this reason, various theoretical predictions have been compared to the full computation, both in a typical VBS fiducial region and also in more inclusive phase space. We have precisely quantified the differences that arise for several physical observables, in particular for the di-jet invariant mass and the rapidity separation of the leading two jets. This is the first time that such an in-depth study is performed. Besides the study of fixed-order predictions, we have also investigated the impact of parton showering. To that end, several LO and NLO event generators which are able to perform matching to parton showers have been employed, and various observables have been thoroughly compared. While in general observables which are described at NLO accuracy show reasonable agreement among the tools, larger differences can appear for those observables described at a lower accuracy, such as those that involve the third jet. In particular such differences are quite prominent in the centralrapidity region, and are the largest for those simulations which employ Pythia8. The effect has been understood, and it can be partially mitigated by changing the recoil scheme of Pythia8 to distribute momenta within initial-final colour connections. These findings make it worth to further investigate these issues not only in the theoretical community, but also by experimental collaborations, for example by measuring related observables for similar processes.
The last part of our work is devoted to remarks and recommendations concerning the usage of theoretical predictions by experimental collaborations.
-As found in Ref. [19], the NLO EW corrections of order O α 7 are the dominant NLO contribution to the process pp → μ + ν μ e + ν e jj + X. It is thus highly desirable to combine them with NLO-QCD predictions matched with parton shower, or at least to include them into experimental analyses. Since, as shown in Ref. [77], these large EW corrections originate from the Sudakov logarithms which factorise, we recommend to combine them with QCD corrections in a multiplicative way. The estimate of missing higher-order EW corrections can be obtained, in a first approximation, by considering ±δ 2 NLOEW , 16 while the missing higher-order mixed QCD-EW corrections can be estimated by taking the difference between the multiplicative and additive prescriptions. For more detailed studies of the combination of QCD and EW higher-order corrections, see e.g. Ref. [107] in the context of top-pair production, or Ref. [108] for SM backgrounds for dark matter searches at the LHC.
-For the typical fiducial region used by experimental collaborations for their measurements, the agreement between the approximations and the full calculation is satisfactory given the current experimental precision, as 16 The quantity δ NLOEW is defined through the relation σ NLOEW = σ LO (1 + δ NLOEW ).
well as the one foreseen for the near future [5,6,109]. Nonetheless, care has to be taken when using such approximations, in particular if more inclusive phasespace cuts are used. -In addition to the standard interpretation of EW signal versus QCD background, combined measurements should also be presented as they are better defined theoretically. In fact, while at LO the interference term can be included in the background component, at NLO the separation of EW and QCD components becomes more blurred, as, e.g. at the order O α s α 6 both types of amplitudes contribute. Therefore, a combined measurement including the EW, QCD, and interference contributions is desirable. Note that with such a measurement a comparison to the SM would be straightforward and still be sensitive to the EW component. In addition, the QCD component could be subtracted based on a well-defined Monte Carlo prediction. -Since the inclusion of NLO QCD corrections gives a better control of extra QCD radiation and reduces the ambiguities related to the matching details and/or the parton shower employed, we encourage the use of NLOaccurate event generators in experimental analyses. In doing so, special care should be employed in order to estimate the theoretical uncertainties, as the standard prescription based on renormalisation and factorisationscale variation is clearly inadequate. Rather, different combinations of generators and parton showers should be employed. -The present study has focused on the orders O α 6 at LO and O α s α 6 at NLO. NLO computations and publicly-available tools also exist for the QCD-induced process [13][14][15][16][17]19,58]