Full NLO predictions for vector-boson scattering into Z bosons and its irreducible background at the LHC

Vector-boson scattering into two Z bosons at the LHC is a key channel for the exploration of the electroweak sector of the Standard Model. It allows for the full reconstruction of the scattering process but at the price of a huge irreducible background. For the first time, we present full next-to-leading-order predictions for $pp \to e^+e^-\mu^+\mu^-jj+X$ including all electroweak and QCD contributions for vector-boson scattering signal and irreducible background. The results are presented in the form of cross sections and differential distributions. A particular emphasis is put on the newly computed $O(\alpha_s^2 \alpha^5)$ corrections.


Introduction
An important part of the physics programme conducted at the Large Hadron Collider (LHC) aims at extracting Standard Model parameters by measuring precisely particular processes. This requires the definition of tailored event selections and the subtraction of undesirable background processes. In the presence of irreducible backgrounds, this picture gets complicated through non-vanishing interferences between the signal and background processes. For vectorboson scattering (VBS), the irreducible background can be overwhelming which warrants a detailed analysis of both the signal and background in order to fully exploit its physics potential [1,2].
The present article continues a series of studies [3][4][5][6] aiming at describing not only VBS processes but also their irreducible backgrounds with next-to-leading-order (NLO) accuracy. Given that the signal and background are connected through interferences, a physical definition of a signal and background separately is not possible (even when using exclusive cuts to enhance the electroweak (EW) component). At NLO, the situation becomes even more complicated as a separation of signal and background in an IR-finite way must necessarily rely on approximations [7].
On top of these considerations, there are more pragmatic reasons to investigate in detail both the EW and QCD components of VBS at full NLO accuracy. While for the same-sign W channel the signal-to-background ratio is about 8/1, it is only of the order of 1/2 for the ZZ case. In addition, given the smallness of the cross sections it is critical to know all contributions, including the loop-induced contribution.
Important steps towards describing pp → e + e − µ + µ − jj + X at full NLO accuracy have already been taken, especially regarding NLO QCD corrections. These are known for the EW component [8] in the VBS approximation and have been implemented in the POWHEG BOX framework [9] to also include parton-shower (PS) corrections [10]. In addition, NLO QCD corrections to the QCD component are known [11] and can in principle be obtained at NLO QCD+PS accuracy from public multi-purpose Monte Carlo generators such as MadGraph5_aMC@NLO [12] or Sherpa [13]. The full NLO corrections of orders O α 7 and O α s α 6 have been computed in Ref. [6] along with the loop-induced contribution of order O α 4 s α 4 . Finally, in Ref. [14], results for loop-induced ZZ production with up to 2 jets merged and matched to parton showers have been presented. Therefore, at NLO level the only missing piece is the order O α 2 s α 5 , which is for simplicity sometimes referred to as EW corrections to the QCD component in the following. With the present publication we fill this gap by presenting for the first time the full NLO corrections to pp → e + e − µ + µ − jj + X at the LHC in a unique setup.
Previous works [3][4][5][6] have put a strong emphasis on EW corrections of order O α 7 . These are exceptionally large for EW corrections and are an intrinsic feature of VBS at the LHC [3]. They are even the largest NLO contribution for the same-sign WW (ss-WW) channel which has thus motivated their implementation in the POWHEG BOX [15]. Such a hierarchy of the NLO contributions is not only due to the general characteristics of VBS but also to the large signal-over-background ratio in the ss-WW channel. In other cases, like for vector-boson scattering into ZZ where the background is large, such a hierarchy does not necessarily hold.
In this article we study these implications and give new insights in VBS and its irreducible background at the LHC.
The results presented in this article should be of great interest for experimental collaborations. In particular, the ATLAS and CMS collaborations have already observed the EW ZZjj production [16][17][18]. Along with state-of-the-art theoretical predictions, upcoming data will allow for deeper analyses of VBS at the LHC.
This article is structured as follows: In Section 2 the process under consideration is introduced. In Section 3 the results are presented and analysed in detail. Finally, Section 4 contains a summary and concluding remarks.

Description of the calculation
The present article complements Ref. [6] by providing the two missing NLO contributions of orders O α 2 s α 5 and O α 3 s α 4 . Therefore, in the following, we often refer to this reference, and the interested reader is invited to look into it. For the sake of simplicity some information is not repeated here as we focus on the salient features of the newly computed contributions.

The process
The physical process under investigation is at the LHC. In the same way as all processes containing VBS contributions, its cross section possesses three leading-order (LO) components: an EW one of order O α 6 including VBS, an interference of order O α s α 5 , and a QCD contribution of order O α 2 s α 4 . Exemplary Feynman diagrams of order O g 6 and O g 2 s g 4 are shown in Fig. 1. At order O g 2 s g 4 , besides partonic processes involving four quarks, also processes with two quarks and two gluons appear (see Fig. 1c) with all possible crossings of gluons and quarks.
At NLO, there are four orders contributing: O α 7 , O α s α 6 , O α 2 s α 5 , and O α 3 s α 4 . The first two have been computed without approximations in Ref. [6]. The latter has been Figure 2: Examples of real tree-level Feynman diagrams: photon emission (left) and gluon emission off LO EW (middle) and QCD diagrams (right). evaluated in Ref. [11], while the order O α 2 s α 5 is calculated here for the first time. Given that the first two orders are discussed in detail in Ref. [6], the discussion in the following is focused on the two remaining orders, which receive contributions from partonic processes involving four quarks as well as from those with two quarks and gluons.
The order O α 2 s α 5 is made of corrections of both QCD and EW types [as the order O α s α 6 ]. In particular, it features both types of real corrections: photon emissions and emissions of QCD partons. The first contribution entails squared matrix elements of photonemission diagrams of order O g 2 s g 5 as shown in Fig. 2a. The second one is made of interferences of QCD real matrix elements of orders O g s g 6 (shown in Fig. 2b) and O g 3 s g 4 (shown in Fig. 2c). Their infrared (IR) singularities are compensated by the corresponding virtual corrections, which also receive two types of contributions. The first one consists of one-loop amplitudes of order O g 2 s g 6 (as shown in Fig. 3a) interfered with tree-level amplitudes of order O g 2 s g 4 (Fig. 1b). The second type is furnished by one-loop amplitudes of order O g 4 s g 4 (as shown in Figs. 3b) interfered with tree-level amplitudes of order O g 6 (Fig. 1a). We note that the two types of NLO corrections cannot be defined in an IR-finite way on the basis of Feynman diagrams in a full computation [4], i.e. without relying on any approximations [7]. For instance, the diagram in Fig. 3a can be viewed as an EW correction to a LO diagram of order O g 2 s g 4 or as a QCD correction to a LO diagram of order O g 6 .
As explained above, at the order O α 2 s α 5 , some of the real corrections are made of a photon emission of a LO QCD amplitude. It also means that in the final state one can have a photon as well as one (or two) gluon(s). In the case where a hard photon is recombined with a soft gluon in a single jet, the QCD singularity associated to the soft gluon is not accounted for by the QED subtraction term. To avoid such configurations, a veto is typically applied on jets that have a too large photon-jet energy fraction where E γ and E a are the energy of the photon and the QCD parton a, respectively. In the present case, the numerical limit has been chosen to be z γ < z γ,cut = 0.7. By cutting into the collinear region, IR-safety is lost but can be restored upon including the nonperturbative fragmentation of quarks into photons via the fragmentation function [19][20][21][22]. Schematically, the NLO cross section in the relevant order can be written as where the index at the integrals denotes the number of particles in the corresponding phase spaces. The modification of the dipoles due to the rejection of hard photons is thus encoded in In addition, the contributions of final-state quarks fragmenting into a hard photon are accounted for by where the sum runs over all the final state quarks and D q i →γ is the quark-to-photon fragmentation function. For the present computation, the fit parameters entering the fragmentation function have been taken from Ref. [23] and read µ 0 = 0.14 GeV, C = −13.26.  (Fig. 1b). Note that we do not consider contributions with external bottom quarks in any of the LO or NLO contributions. These are small or can be removed in experimental analyses with the help of bottom-jet vetoes [5,6].

Details and validation
The results presented here have been produced using the Monte Carlo program MoCaNLO and the matrix-element generator Recola. The program MoCaNLO is able to compute arbitrary processes within the SM at NLO QCD and EW accuracy. In order to obtain a fast integration for high-multiplicity processes with many resonances, it takes advantage of phasespace mappings similar to the ones of Refs. [24][25][26]. On the other hand, Recola [27,28] is a general tree and one-loop matrix-element provider. It relies on the Collier library [29,30] that provides numerically the one-loop scalar [31][32][33][34] and tensor integrals [35][36][37]. This combination has already shown to work well for high-multiplicity processes (2 → 6 and beyond). Many of these applications concerned VBS processes [3][4][5][6][7]38]. During the course of these studies, the set of tools was verified against BONSAY+OPENLOOPS [5] for WZ scattering, the Monte Carlo BBMC for same-sign W scattering [3,4], and against Powheg [15] for same-sign W scattering at order O α 7 . For non-VBS processes, it was checked against Sherpa in Refs. [39,40] and a multitude of tools in Ref. [41] for di-boson production at NLO EW.
Finally, in Ref. [6], the partonic processes ud → e + e − µ + µ − ud, us → e + e − µ + µ − dc, uu → e + e − µ + µ − uu, and uc → e + e − µ + µ − ds at order O α 7 were successfully compared between MoCaNLO+Recola and BBMC+Recola. Also, selected contributions of order O α s α 6 were found in agreement within integration errors. For the present article, about 10 representative partonic channels were confirmed to agree between MoCaNLO and BBMC at the NLO orders O α 2 s α 5 and O α 3 s α 4 within statistical errors at the level of a few per cent. In MoCaNLO and BBMC, the subtraction of IR divergences in the real radiation is realised with the help of the Catani-Seymour dipole formalism [42][43][44]. Within this method, it is possible to vary the α dipole parameter [45], which restricts the subtraction to IR-singular regions of the phase space. Obtaining consistent results with different values α dipole ensures the correctness of the subtraction mechanism. We have therefore performed two full computations of all NLO orders using MoCaNLO+Recola with α dipole = 1 and α dipole = 10 −2 . Full statistical agreement has been found, and the results presented here are the ones obtained with α dipole = 10 −2 . Also, for the computation of order O α 2 s α 5 two different modes of Collier have been successfully compared. The results presented here have been obtained with the COLI mode throughout. In all computations, the complex-mass scheme [25,[46][47][48] is utilised.
Finally, the implementation of the fragmentation function closely follows Refs. [19][20][21][22]. Its implementation in MoCaNLO has been validated against BBMC which was used in Ref. [22] for the computation of the EW corrections to pp → + − jj + X. Also, this implementation was employed in Ref. [40] for pp → e + ν e µ −ν µ j + X where it was compared against an approximate treatment.

Input parameters and event selection
The input parameters and event selection used in the present article are exactly the same as those in Ref. [6]. For completeness we reproduce them here.

Input parameters
The calculation is done for LHC characteristics and a centre-of-mass energy of 13 TeV. The set of parton distribution functions (PDF) NLO NNPDF-3.1 Lux QED with α s (M Z ) = 0.118 [49,50] and with fixed N F = 5 flavour scheme is used throughout. All collinear initial-state splittings are treated by the MS redefinition of the PDF. The PDF set is interfaced to our Monte Carlo programs using LHAPDF [51,52].
The renormalisation and factorisation scales, µ ren and µ fac , are chosen to be for all contributions, where j 1 and j 2 are the two hardest (in p T ) identified jets (tagging jets). This is the scale choice that is usually used for VBS processes in the literature and that has been also employed in Ref. [6]. In cross sections and differential distributions, the scale uncertainty is obtained with the standard 7-point scale variation, i.e. observables are evaluated for the 7 pairs of renormalisation and factorisation scales and the scale variation is determined from the resulting envelope. The electromagnetic coupling is obtained through the G µ scheme [53] via where G µ is the Fermi constant. The numerical values of the masses and widths used as input in the numerical simulation read [54] m t = 173.0 GeV, Γ t = 0 GeV, Everywhere, it is assumed that m b = 0 GeV and as mentioned above no partonic channels with external bottom quarks are considered. For the massive gauge bosons (W and Z), the pole masses and widths utilised in the calculation are obtained from the on-shell (OS) values [55] using (3.5)

Event selection
The event selection of the present computation is largely inspired from the CMS analyses [17,18]. The process features four charged leptons and two jets in the final state. Quarks and gluons with pseudorapidity |η| < 5 are clustered into jets with the anti-k T algorithm [56] using R = 0.4. Following the same method and radius parameter, the photons are recombined with the final-state jets or leptons. The four leptons must fulfil where and are any leptons while + and − are oppositely charged leptons regardless of their flavour. On top of these cuts, the invariant masses of the leptonic decay products of the Z bosons are restricted to Once the jet clustering is performed, at least two jets still have to fulfil the criteria Such jets are then denoted as identified jets. Out of these, the two with the highest transverse momenta are called tagging jets. In our default setup (inclusive setup for short) the tagging jets have to fulfil the constraint We consider in addition a setup (VBS setup) where the last cut is replaced by a stronger one M j 1 j 2 > 500 GeV. (3.10)

Cross sections
In this section, various cross sections at LO and NLO accuracy are discussed. We start by recalling the LO cross sections computed in Ref. [6] in Table 1. As opposed to Ref. [6], here the sum comprises only the contributions of orders O α 6 , O α s α 5 , and O α 2 s α 4 , while the loop-induced contribution of order O α 4 s α 4 from gg → e + e − µ + µ − gg does not enter the sum and is only shown for completeness. Contributions of partonic channels with four external quarks, σ 4q LO , and of those involving gluons, σ 2q2g LO , which enter at order O α 2 s α 4 , are shown separately as well. When referring to relative corrections in the following, they are always normalised to the above LO sum, and the loop-induced contribution is not included in the NLO predictions. Table 1 clearly shows that the QCD contribution of order O α 2 s α 4 is dominating and reaches 91% of the LO sum for M j 1 j 2 > 100 GeV. When applying the tighter VBS event selection, M j 1 j 2 > 500 GeV, it is reduced to 63%, demonstrating hence the impact of such cuts in enhancing the EW contribution. While the partonic channels involving quarks and 0.097683(2) 0.008628 (1) (2)  for M j 1 j 2 > 100 GeV, for M j 1 j 2 > 500 GeV, (3.12) i.e. the scale dependence is reduced by more than a factor 3 when including NLO corrections.
Order  In Table 2 Table 1, two setups with M j 1 j 2 > 100 GeV and M j 1 j 2 > 500 GeV are displayed. For the case M j 1 j 2 > 100 GeV, the largest correction is the one of order O α 2 s α 5 with −7.1%, dominated by the EW corrections to the LO QCD-induced contribution, followed by the one of order O α 3 s α 4 with −4.5%, the corresponding QCD corrections. The corrections of order O α 7 and O α s α 6 are on the other hand between 1% and 2%. The picture changes when imposing the restriction M j 1 j 2 > 500 GeV. Then, the two largest NLO contributions are the EW corrections of orders O α 7 and O α 2 s α 5 being both negative and between −6% and −6.5%. The contributions of order O α 2 s α 5 are dominated by partonic channels involving two quarks, which amount to 80% for the inclusive setup and 66% for the VBS setup. At O α 3 s α 4 , the four-quark channels dominate for M j 1 j 2 > 100 GeV, while the gluon-quark channels take over for M j 1 j 2 > 500 GeV. In the sum of all NLO corrections contributions from four-quark and gluon-quark channels are of similar size.
The hierarchy of NLO contributions is quite different from the one in pp → µ + ν µ e + ν e jj+X, which was computed in Ref. [4]. There, the largest NLO corrections were those of order O α 7 with −13.2% followed by those of order O α s α 6 with −3.5%, while the contributions of orders O α 2 s α 5 and O α 3 s α 4 were subleading and below 1%. This different hierarchy between the NLO corrections for the two signatures is linked to the rather different LO hierarchy between the EW and QCD components, which originates to a considerable extent from the appearance of partonic channels involving gluons, which are not present for ss-WW. While the EW component of order O α 6 represents 86.5% of the total LO in the case of ss-WW, it contributes only 34.2% in the case of ZZ for the most exclusive setup. The LO hierarchy strongly determines the one of the NLO corrections. For example, when normalising the corrections of order O α 7 for VBS into ZZ to the LO EW contribution of order O α 6 , one obtains −17.6% in line with the results of Refs. [3][4][5].
Interestingly s α 4 only slightly from −8.2% to −7.9% for M j 1 j 2 > 100 GeV and from −11.6% to −10.3% or M j 1 j 2 > 500 GeV. The situation is different for some partonic channels (see below). Note also that for ZZ the behaviour of the O α s α 5 corrections is largely determined by the partonic channels involving gluons, which are absent for ss-WW.
Some qualitative understanding of the NLO EW corrections of order O α 2 s α 5 can be obtained from an high-energy approximation of the virtual corrections in the pole approximation, i.e. for the process pp → ZZjj → e + e − µ + µ − jj + X. To keep it simple, we restrict the discussion to the dominant contributions of left-handed quarks and transverse Z bosons and consider only the double EW logarithms, the collinear single EW logarithms, and the single EW logarithms resulting from parameter renormalisation. We do not include the angular-dependent leading logarithms, which have a much more complicated structure. Based on Ref. [57], we find for the EW correction factor to the cross section of order O α 2 s α 4 for qq → ZZqq (n q = 2) and qg → ZZqg (n q = 1) in the leading-logarithmic approximation: The constants are given by and where I 3 w,q and Q q denote the weak isospin and relative charge of the quark q, respectively, while s w and c w stand for the sine and cosine of the EW mixing angle. The correction factor (3.13) also applies unchanged to processes resulting from crossing symmetry and/or if one of the quark pairs is replaced by a different quark pair with the same quantum numbers, e.g. u → c. For processes involving quarks with different quantum numbers, the simple formula (3.13) does not hold, however, in the limit of vanishing hypercharge it is valid for arbitrary quark combinations with F q = 1. Since the correction factor (3.13) only describes EW corrections to the LO cross section of order O α 2 s α 4 , this approximation does not apply if QCD corrections to the LO interference of order O α s α 5 are sizeable, which is usually the case if the LO EW cross section of order O α 6 is comparable or larger than the QCD one of order O α 2 s α 4 . This happens, in particular, for the partonic processes dd → e + e − µ + µ − uū, uū → e + e − µ + µ − dd, du → e + e − µ + µ − du, anddū → e + e − µ + µ −dū , where the contributions of orders O α 2 s α 4 and O α 6 are comparable. The only free parameter of (3.13) is the scale Q. While for a 2 → 2 process one usually picks the centre-of-mass energy, a 2 → 4 process involves many more scales. As a typical scale, we choose the invariant mass of the ZZ pair, Q = M ZZ = M 4 , which is the scale appearing in the leading double logarithms multiplied with the largest coupling factor. Determining the average of M 4 from LO distributions, yields values for Q in the range 330-470 GeV for all types of processes in the inclusive setup. Using Q = M 4 in (3.13) event by event results in correction factors that agree with the full calculation within about 2% for all partonic processes with external gluons, which dominate the cross section at O α 2 s α 4 , and within 3% for most four-quark processes. In the VBS setup, we find average values for Q within 390-550 GeV for four-quark processes and 370-530 GeV for gluon-quark processes as well as correction factors δ LL that agree with the full relative corrections within about 3% for the gluonic processes and within 4% for most four-quark processes. The approximation fails completely for the four-quark processes with sizeable O α 6 contributions mentioned at the end of the preceding paragraph. Besides those, the approximation is worse than the above-stated figures only for dd → e + e − µ + µ − cc and uū → e + e − µ + µ − ss, owing to somewhat larger QCD corrections to the EW LO interference, as well as for uu → e + e − µ + µ − uu, which has the highest scales M 4 among all partonic processes. In general, the logarithmic approximation predicts larger corrections for processes involving up-type quarks as compared to those with down-type quarks, which is also observed in the full calculation. In summary, the simple approximation (3.13) with the scale choice Q = M 4 reproduces the O α 2 s α 5 corrections to the fiducial cross section within a few per cent for all partonic processes that are dominated by the LO O α 2 s α 4 contributions and for the sum of all partonic channels.

Differential distributions
In this section, results for differential distributions at 13 TeV are discussed for the inclusive setup defined in Section 3.1. In the first part, we display the contributions of the different orders separately. In the upper panels of Figs The first set of differential observables in Fig. 4 relates to the two tagging jets as defined in Eqs. (3.8) and (3.9). The distributions in the invariant mass (Fig. 4a) and the rapidity separation (Fig. 4b) of the two tagging jets are typically used to improve the ratio of the EW component over the QCD one. These two distributions differ substantially from all others shown in this paper, since they receive sizeable contributions from the order O α 6 for large M j 1 j 2 or large ∆y j 1 j 2 , as can be seen in the upper part of both plots, while all other distributions are dominated by the order O α 2 s α 4 throughout. Thus, the normalisation of the relative corrections is dominated by the O α 2 s α 4 contributions for small M j 1 j 2 and ∆y j 1 j 2 , but by the O α 6 ones for large variables. Owing to this varying normalisation, the EW corrections of order O α 7 are large for large M j 1 j 2 or large ∆y j 1 j 2 (reaching −18% at M j 1 j 2 = 2 TeV) and small otherwise. The normalisation also explains the opposite behaviour of the (EW) corrections of order O α 2 s α 5 , which reach about −15% at M j 1 j 2 = 400 GeV but are reduced to about −5% to −7% at 2 TeV in the invariant-mass distribution. Despite the fact that these large EW corrections can be traced back to Sudakov logarithms, they become relatively smaller at high energies as the LO contribution of order O α 2 s α 4 (to which these corrections act on) is suppressed there. The corrections of QCD type, O α s α 6 and O α 3 s α 4 , stay within ±10% apart from the region of large rapidity separations. In particular, the QCD corrections of order O(α 3 s α 4 ) turn very large (over +40%) there, but this part of the phase space is rather suppressed. These QCD corrections are positive for invariant masses above M j 1 j 2 = 500 GeV and tend to increase towards higher invariant masses.
The corrections to the distributions in the azimuthal angle (Fig. 4c) and the cosine of the angle between the two tagging jets (Fig. 4d) show rather mild variations that do not exceed 10% over the full phase-space range. The EW corrections of order O(α 7 ) vary by less than 2% with the exception of cos θ j 1 j 2 ≈ −1. The EW corrections of order O(α 2 s α 5 ) are almost −10% for small angles between the jets, while they decrease in magnitude to −5% for large angles. The pure QCD corrections of order O(α 3 s α 4 ) grow in the azimuthal-angle distribution from −10% up to almost 0%, while they are minimal for small values of cos θ j 1 j 2 and maximal for cos θ j 1 j 2 = ±1.
In Fig. 5, several leptonic observables are shown, including distributions in the invariant mass (Fig. 5a) and transverse momentum (Fig. 5b) of the four leptons as well as in the invariant mass (Fig. 5c) and the transverse momentum (Fig. 5d)  negative QCD corrections. Note that the O(α s α 6 ) corrections, i.e. the QCD corrections to the VBS contributions, stay within ±10% for all distributions in Fig. 5. The distribution in the invariant mass of the electron-positron pair (Fig. 5c) is characterised by the Z-boson resonance at ∼ 91 GeV in all LO contributions. While the QCD corrections are rather flat across the resonance, the EW corrections of order O(α 2 s α 5 ) exhibit a pronounced radiative tail below the resonance. The corrections due to final-state photon radiation exceed +30% around 70 GeV. The corresponding radiative tail is also present in the corrections of order O(α 7 ) but suppressed by one order of magnitude by the LO QCD contributions.
In Fig. 6 we present distributions related to the leading jet or a single lepton. In the distribution in the transverse momentum of the hardest jet (Fig. 6a), the corrections of order O(α 3 s α 4 ) show the most interesting behaviour. They start around 6% at 30 GeV, reach a minimum of about −12% around 350 GeV to finally grow to 1% at 800 GeV. The moderate rise of these corrections towards high transverse momenta can be attributed to the choice 3.1 of renormalisation scale. The increase of the QCD corrections towards small p T,j 1 has already been observed in other VBS/VBF processes [4][5][6]58] and is due to the real QCD radiation. The corrections of order O(α 2 s α 5 ) decrease steadily towards high energy by about 10%. Regarding the rapidity distribution of the hardest jet (Fig. 6b), the corrections of order O(α 2 s α 5 ) vary only by few per cent in the whole kinematic range. The O(α 3 s α 4 ) corrections, on the other hand, go from +8% in the peripheral region down to −7% in the central one. The differential distributions in the transverse momentum of the muon (Fig. 6c) and the rapidity of the positron (Fig. 6d) are rather standard. The p T,µ − distribution behaves like other transverse-momentum distributions of leptons shown above, while the variation of all corrections to the y e + distribution is below about two per cent.
Next we consider leptonic angular distributions in Fig. 7. The fraction of the LO EW contribution is somewhat enhanced for small and moderate cos θ e + µ − . In all distributions of  4 which have already been shown in Ref. [6]. The total corrections to the distribution in the invariant mass of the two tagging jets (Fig. 8a) vary between −10% and −30% and are smallest at about 800 GeV. These corrections are not dominated by one particular contribution but result from the interplay of the four NLO contributions (see Fig. 4a). The scale uncertainty is roughly of the same size as for the fiducial cross section. On the other hand, the corrections to the distribution in the rapidity difference of the two tagging jets (Fig. 8b)  jet (Fig. 8c) range between 0% and −25%, where the overall behaviour is directed by the one of the O(α 3 s α 4 ) corrections. It is worth noting that the NLO prediction is beyond the edge of the LO scale-uncertainty band for p T,j 1 > 200 GeV and that the scale uncertainty for this distribution is rather large. As for all other distributions, the central value is close to the upper edge of the scale envelope. The behaviour of the NLO corrections to the distribution in the transverse momentum of the second hardest jet (Fig. 8d) is quite different from the one of the hardest jet. The NLO corrections are at the level of −20% at 30 GeV and reach a minimum at 100 GeV close to 0% to steadily decrease to −20% at 800 GeV. The NLO scale uncertainty stays small over the whole considered kinematic range.  O α 2 s α 4 contributions at high leptonic energies, where the scales are underestimated. The azimuthal angle between electron and positron is correlated to the transverse momentum of the electron-positron pair. Accordingly, the corrections are small for ∆φ e + e − ≈ 180 • and increase smoothly to −30% for small ∆φ e + e − , where the corresponding rate is minimal and the scale uncertainty is about 20%. Since the correlation of the angle between the positron and the muon to the transverse momentum of the electron-positron pair is smaller, the corrections vary only from −5% to −20% with increasing angle. Finally, we show the distribution in the total transverse energy, defined as the sum of transverse energies of the tagged final-state objects, H T = E T,e + + E T,e − + E T,µ + + E T,µ − + E T,j 1 + E T,j 2 , (3.16) which are defined from the transverse momenta p T,i as E T,i = p 2 T,i + M 2 i , (3.17) where the invariant mass M i is nonzero after recombination. The large negative corrections below 400 GeV result exclusively from the order O α 3 s α 4 corrections, more precisely from a suppression of the corresponding real contributions resulting from the restriction of the related phase space for small H T . The smooth increase of the corrections to 20% at 2 TeV is basically driven by the positive O α 3 s α 4 corrections combined with negative O α 2 s α 5 corrections. The scale uncertainty remains small apart from the region of very small H T .

Conclusion
Given the current and upcoming expected precision of the LHC experiments, VBS processes offer great opportunities to further probe the EW sector of the Standard Model. In this context, precision does not only refer to the EW signal that contains VBS contributions but also to the irreducible background that can be overwhelming. In addition of being sometimes large, the background can not be trivially separated from the signal making it a crucial component of any VBS studies. This calls for a full NLO description of VBS signatures including all possible contributions.
In this article we have followed this avenue by presenting full NLO predictions for pp → e + e − µ + µ − jj + X. While some of the NLO contributions were already known in the literature, they are shown together in a unique setup for the first time here. This allows to single out salient features of VBS into ZZ pairs. It is worth noting that the hierarchy of NLO corrections is rather different from the one in the ss-WW case (which was the only channel known at full NLO accuracy up to now), in particular, owing to the appearance of partonic channels with gluons. Comparing two different setups, with different invariant-mass cuts on the two tagging jets, we have demonstrated that the LO hierarchy strongly influences the size of the NLO corrections. The newly computed corrections of order O α 2 s α 5 turn out to be much more sizeable here than for the ss-WW signature. While for the ss-WW cross section, they were found to be below a per cent, they are between −6% and −8% (depending on the setup) for the ZZ case discussed here. We have traced back these corrections to typical EW Sudakov logarithms which appear for both signatures. While these corrections contribute for VBS into ZZ with their natural order of magnitude, they were accidentally compensated for ss-WW by QCD corrections to the LO interference appearing at the same order. This cancellation does not happen here due to the larger LO QCD contribution upon which these EW corrections act. This has an important implication: in arbitrary experimental setups, none of the four NLO corrections can safely be neglected at the 10% level. Indeed, the impact of the various NLO corrections strongly depends on the hierarchy between the LO contributions which eventually originates from the experimental event selections.
In high-energy tails of distributions, the O α 2 s α 5 corrections reach 30%, while the O α 3 s α 4 corrections can amount to 50%. Angular distributions are distorted by up to 10% and 25% by the corrections of orders O α 2 s α 5 and O α 3 s α 4 , respectively. Owing to the enhancement from Sudakov logarithms and the scale choice adapted to the VBS contributions, the total NLO corrections exceed the LO scale uncertainty band in various regions of phase space.
The results presented here should prove particularly useful for current and upcoming analyses for ZZ production in association with two jets at the LHC. We hope that experimental collaborations will take into account such NLO corrections, paving the way to precise comparisons between experimental measurements and theoretical predictions. We would like to emphasise that tailored corrections can be provided upon request.