NLO QCD and EW corrections to vector-boson scattering into ZZ at the LHC

We present the first calculation of the full next-to-leading-order electroweak and QCD corrections for vector-boson scattering (VBS) into a pair of Z bosons at the LHC. We consider specifically the process pp → e+e−μ+μ−jj + X at orders O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(α7) and O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(αsα6) and take all off-shell and interference contributions into account. Owing to the presence of enhanced Sudakov logarithms, the electroweak corrections amount to −16% of the leading-order electroweak fiducial cross section and induce significant shape distortions of differential distributions. The QCD corrections on the other hand are larger (+24%) than typical QCD corrections in VBS. This originates from considering the full computation including tri-boson contributions in a rather inclusive phase space. We also provide a leading-order analysis of all contributions to the cross section for pp → e+e−μ+μ−jj + X in a realistic setup.


Introduction
The non-observation of physics not described by the Standard Model (SM) at the Large Hadron Collider (LHC) shifts into the focus of interest tests of the SM of particle physics and in particular of the mechanism of electroweak (EW) symmetry breaking. A particularly important class of processes in this respect is vector-boson scattering (VBS). It is not only sensitive to the scalar sector of the EW theory but also to non-standard triple and quartic gauge-boson couplings.
The scattering of like-sign W bosons as well as WZ has been observed [1][2][3][4][5][6][7][8] by both ATLAS and CMS in leptonic final states. Quite recently, measurements of VBS into a pair of leptonically decaying Z bosons, which allows for a cleaner experimental detection albeit with lower hadronic cross section, have been published by both experiments [9][10][11]. Moreover, EW di-boson production in association with a high-mass di-jet system in semileptonic final states has been searched for [12,13].
Theoretically, VBS has been studied since a long time. NLO QCD corrections to all VBS processes exist for more than ten years (see refs. [14,15] and references therein). On the other hand, EW corrections to these processes have become available only recently [16][17][18]. Thereby it has been found that EW corrections to fiducial cross sections of VBS are generically at the level of −15% [16]. This has been confirmed by complete calculations for the scattering of W + W + [17] and WZ pairs [18]. Furthermore, for like-sign W scattering an event generator including EW and QCD corrections is available [19].
In this article, we present for the first time NLO EW corrections to VBS into a pair of Z bosons at the LHC, considering specifically the fully leptonic final state e + e − µ + µ − jj.

JHEP11(2020)110
The prospects of this channel for the high-luminosity and high-energy upgrade of the LHC have been studied in ref. [20]. NLO QCD corrections to the purely EW VBS process in the so-called VBS approximation, which neglects s-channel diagrams and interferences between u-channel and t-channel diagrams, have been presented in ref. [21]. These corrections have been matched [22] to a QCD parton shower via the POWHEG BOX framework [23]. NLO QCD corrections to the corresponding QCD-induced process have been provided in ref. [24]. Very recently, loop-induced ZZ production with up to 2 jets merged and matched to parton showers has been studied [25].
The leading-order (LO) cross section for pp → e + e − µ + µ − jj + X receives contributions of orders O α 6 (EW contributions), O α s α 5 (interference), and O α 2 s α 4 (QCD-induced contributions). The VBS subprocess is part of the gauge-invariant EW contributions. In the experimental analysis, the three contributions are often referred to as signal, interference, and QCD background, respectively. While the QCD-induced contributions are larger than the EW ones for typical VBS cuts, the interference is smaller than the latter. Since gluon-induced contributions at the order O α 4 s α 4 contribute sizeably, we include them in our analysis. In this paper we focus on the NLO EW and QCD corrections to the LO EW contributions. More precisely, we consider the complete gauge-invariant set of contributions at orders O α 7 and O α s α 6 . While the former are pure EW corrections to the LO EW contributions, the latter receive contributions from the QCD corrections to the LO EW contributions as well as from EW corrections to the LO interferences. In our NLO calculation we do not rely on approximations but take into account the full set of diagrams relevant at the corresponding perturbative order including all interference contributions. We also compare the full NLO EW corrections with the ones obtained within a Sudakov approximation.
The final state e + e − µ + µ − jj leads to a wide-spread variety of contributing partonic channels. This renders the full calculation at NLO EW and especially at NLO QCD technically very demanding and CPU-time intensive.
This article is structured as follows: in section 2 the considered process and the different contributions are described. In addition, details of our calculation including the checks to validate our results are presented. Section 3 contains numerical results and their description. Lastly, in section 4 a summary and concluding remarks are given.

Leading-order contributions
We are studying the process pp → e + e − µ + µ − jj + X. (2.1) At LO, VBS appears in quark-induced partonic channels qq → e + e − µ + µ − qq (q generically stands for a quark or anti-quark). The amplitudes for these processes receive contributions of order O g 6 as well as of order O g 2 s g 4 , where g and g s denote the EW and strong coupling constant, respectively. JHEP11(2020)110 We neglect quark mixing and use a unit quark-mixing matrix. Since a non-trivial quark-mixing matrix would only affect the s-channel contributions and the NLO corrections, its effects are suppressed.
Some sample diagrams of order O g 6 are shown in figures 1(a)-1(h). The diagrams in figures 1(a) and 1(b) represent the characteristic t-channel VBS topology. While the first diagram illustrates the scattering of W + W − into ZZ, the second one exemplifies ZZ scattering, which takes place only via Higgs exchange. Higgs-exchange s-channel diagrams arise in VBS into ZZ as a new feature with respect to same-sign W and WZ scattering. 1 However, partonic channels involving only ZZ scattering are suppressed with respect to those involving the subprocess W + W − → ZZ. Doubly-resonant, singly-resonant, and nonresonant diagrams of non-VBS type contribute to the partonic processes at order O g 6 , as illustrated in figures 1(c), 1(d), and 1(e), respectively. Besides t-and u-channel diagrams also s-channel diagrams show up for some partonic processes (figure 1(f)). In particular, s-channel diagrams corresponding to triple gauge-boson production, both for WZZ and ZZZ production appear as depicted in figures 1(g) and 1(h). Lastly, diagrams at the order O g 4 g 2 s are characterised by t-channel gluon exchange between the two quark lines 1 For our set of cuts, the Higgs resonance does not appear in the fiducial phase-space region. JHEP11(2020)110 (see figure 1(i)). While in the diagrams we specify only Z bosons in s-channel propagators to stress the contributions to the intermediate ZZ state, corresponding diagrams with photons are also included in our computation. In fact, we take into account the complete set of diagrams for the given six-fermion final state. Consequently, at the level of squared amplitudes three kinds of gauge-invariant contributions exist, purely EW contributions at O α 6 , QCD-induced contributions at order O α 2 s α 4 and contributions at order O α s α 5 , which result from the interference of diagrams at order O g 6 and O g 2 s g 4 . Owing to the colour structure, the latter contributions are only non-vanishing, if diagrams of two different kinematic channels (s, t, u) with all 4 external quarks in the same generation are interfered, such as the diagrams in figure 1(f) and figure 1(i).
We do not consider contributions with bottom quarks in the initial state, which are PDF suppressed, and also do not include final states with bottom quarks. 2 Then 60 partonic quark-induced channels contribute compared to 40 for WZ and 12 for W ± W ± scattering (not counting qq and q q initial states separately). Out of these 60 channels, 24 receive non-vanishing interference contributions between different coupling orders that make up the contribution of order O α s α 5 . At order O α 2 s α 4 , in addition to the 60 quark-induced channels, channels with one or two gluons in the initial state contribute. Given the large gluon luminosity at the LHC, the latter are one of the reasons for the enhancement of the QCD-induced contributions over the EW ones.
Further contributions at orders O α 6 and O α s α 5 result from photon-induced processes with γγ, γg and γq initial states. Such contributions were found to be below 0.5% for WZ scattering [18], which is also expected for VBS into ZZ. These contributions are neglected in this work.
In contrast to final states corresponding to charged W ± W ± and WZ scattering, the e + e − µ + µ − jj final state receives contributions from the loop-induced partonic process gg → e + e − µ + µ − gg at order O α 4 s α 4 (see figure 2 for sample diagrams). We include these contributions in our leading-order analysis.

Virtual corrections
We compute NLO corrections of orders O α 7 and O α s α 6 to the process (2.1).
Virtual corrections of order O α 7 result from interfering purely EW loop diagrams of order O g 8 with tree diagrams of order O g 6 and are thus EW corrections to the JHEP11(2020)110 s . As a consequence EW and QCD-induced contributions cannot be separated at order O α s α 6 on the basis of Feynman diagrams. While diagrams with gluons attached to a single quark line (see figure 3(c)) contribute for all partonic channels, diagrams with gluon exchange between different quark lines only contribute for partonic processes that receive contributions of different kinematic channels such as those in figures 3(e) and 3(f), which interfere with the LO diagram figure 1(a). In the VBS approximation only t-and u-channel diagrams are taken into account and no interferences between different kinematic channels. Then all corrections of O α s α 6 can be interpreted as QCD corrections to the EW LO diagrams.
We do not include contributions of partonic channels with external bottom quarks or photons in the initial state in the virtual corrections.

Real corrections
In addition to the virtual corrections also real photon and gluon emission needs to be considered. Some related Feynman diagrams are shown in figure 4.
In the quark-induced channels, emission of a real photon from a quark, a W boson, or a charged lepton gives rise to the partonic processes qq → e + e − µ + µ − qqγ. These furnish the real EW corrections to the LO EW diagrams with diagrams of order O g 7 . As exam-JHEP11(2020)110 As the virtual corrections, also the real contributions at order O α s α 6 emerge from different sources. First, as part of the NLO QCD corrections, emission of a real gluon from quarks of the LO EW diagrams exemplarily depicted in figure 4(c), results in diagrams of order O g 6 g s . These contribute to the partonic process qq → e + e − µ + µ − qqg and upon squaring yield contributions of order O α s α 6 . Contributions of the same order result from interferences of real-photon-emission diagrams of orders O g 7 (figure 4(b)) and O g 5 g 2 s (figure 4(d)) corresponding to different kinematic channels. Besides quarkinduced channels also gq channels appear with diagrams resulting from crossing the gluon and a quark in qq → e + e − µ + µ − qqg, as for example shown in figures 4(e) and 4(f), and contribute at order O α s α 6 .
The number of partonic channels for the quark-induced processes is the same as for the corresponding LO processes. For the gluon-induced gq processes 40 partonic channels contribute compared to 28 for WZ and none for W ± W ± scattering. Note that we do not take into account partonic channels with external bottom quarks.
Crossing the photon and a quark in qq → e + e − µ + µ − qqγ gives rise to photon-induced partonic processes that contribute both at orders O α 7 and O α s α 6 . Such contributions have been found to be below 2% for same-sign W scattering [17] and are neglected in this calculation.
For the treatment of the infrared (IR) singularities arising in the phase-space integration of the real corrections we use Catani-Seymour dipole subtraction for QCD [26] and JHEP11(2020)110 its variant for QED [27,28]. In this regard, the strategy for WZ scattering of ref. [18] can be taken over, and no new features appear. At the order O α s α 6 both QCD and QED IR singularities due to soft and/or collinear gluon or photon emission or via forward branchings of QCD partons in the initial state appear. Another type of singularity arises in diagrams with photons converting into a quark-anti-quark pair in the final state. The low virtuality of the photon leads to a collinear singularity, which would cancel against the virtual corrections to the e + e − µ + µ − jγ final state. While this final state is considered separately from VBS, this singularity can be absorbed into a photon-to-jet conversion function [29], which can be related to the non-perturbative hadronic vacuum polarisation. In general, different kinds of singularities appear in the same diagrams and have to be dealt with simultaneously. More details on the treatment of IR singularities can be found in refs. [17,18].

Details of the computation and validation
The results presented here have been obtained from the combination MoCaNLO+Recola of two computer programs. MoCaNLO is a generic Monte Carlo integration code that can compute arbitrary processes at NLO QCD and EW accuracy in the SM. The use of phasespace mappings similar to the ones of refs. [30][31][32] ensures an efficient integration even for high-multiplicity processes. Recola [33,34], on the other hand, is a tree and one-loop matrix-element provider. It uses the Collier library [35,36] to obtain numerically the oneloop scalar [37][38][39][40] and tensor integrals [41][42][43]. The tandem MoCaNLO+Recola has already been used for many complex computations, in particular for VBS processes [16-18, 44, 45]. It has also been successfully compared for WZ scattering with results of BONSAY+OPENLOOPS [18] and for same-sign W scattering with the combination of Recola and the Monte Carlo integration code BBMC [17]. MoCaNLO+Recola has also served as baseline for the validation of the implementation of O α 7 corrections in Powheg+Recola [19] for same-sign W scattering and of O α s α 5 corrections for pp → WWj in Sherpa+Recola [46,47]. For VBS into ZZ, we compared the dominant partonic process ud → e + e − µ + µ − ud as well as the channels us → e + e − µ + µ − dc, uu → e + e − µ + µ − uu, and uc → e + e − µ + µ − ds at order O α 7 and selected contributions at order O α s α 6 between MoCaNLO+Recola and BBMC+Recola and found agreement within integration errors.
In MoCaNLO, the IR divergences in the real radiation are subtracted using the dipole method for QCD [26] and its extension to QED [27,28]. Within the dipole-subtraction scheme, the α dipole parameter [48] allows to restrict the phase space to the singular regions. The value α dipole = 1 corresponds to the full phase space (within the acceptance defined by selection cuts) without additional restrictions. In the present case, two full computations both at order O α 7 and O α s α 6 have been performed, one with α dipole = 10 −2 and one with α dipole = 1. Both computations agree within statistical errors at each order and provide thus a robust check of the subtraction procedure. The results presented here have been obtained from the computation with α dipole = 10 −2 . In addition, representative contributions have also been computed with two different numerical values of the IRregulator parameter, proving IR finiteness of the results. Finally, for the computation JHEP11(2020)110 of the one-loop amplitudes of order O g 6 g 2 s , the two different modes of Collier have been used: the DD mode in the α dipole = 1 computation and the COLI mode in the α dipole = 10 −2 computation. Throughout, the massive resonant particles are described in the complex-mass scheme [31,49,50] at both tree and one-loop level.

Input parameters and event selection
Input parameters. The numerical simulations have been carried out for the LHC at a centre-of-mass energy of 13 TeV. As parton distribution functions (PDF) we employ the NLO NNPDF-3.1 Lux QED set with α s (M Z ) = 0.118 [51,52]. This is incorporated in the calculation through LHAPDF [53,54]. We are working with the fixed N F = 5 flavour scheme throughout. Both the EW and QCD collinear initial-state splittings are treated by MS redefinition of the PDF. Note that the same PDF set is used for both the LO and NLO predictions.
The renormalisation and factorisation scales are set to the geometric average of the transverse momentum of the jets where j 1 and j 2 are the two hardest identified jets (see definition below) ordered according to transverse momentum. In the following, we perform the usual 7-point scale variation of both the renormalisation and factorisation scale, i.e. we calculate the observables for the pairs of renormalisation and factorisation scales with the central scale defined in (3.1) and use the resulting envelope.
To fix the electromagnetic coupling, the G µ scheme [55] is used, which defines the coupling from the Fermi constant as The masses and widths of the massive particles are chosen as [56] m t = 173.0 GeV, The bottom quark is taken to be massless, and no partonic channel with initial-state and/or final-state bottom quarks is included. Since there are no resonant top quarks in the considered processes, we set the top-quark width to zero. The values of Higgs-boson mass and width are taken from ref. [57]. The pole masses and widths of the W and Z bosons used in the calculation are determined from the measured on-shell (OS) values [58] via (3.5) Event selection. The event selection used here is inspired by the CMS measurement [9]. Experimentally, the final state of the process is described by four charged leptons and at least two QCD jets. The QCD partons are combined into jets with the anti-k T algorithm [59] using R = 0.4 as jet-resolution parameter. In the same way, the real photons are recombined with the final-state quarks into jets or with the charged leptons into dressed leptons. In both cases the anti-k T algorithm and a resolution parameter R = 0.4 is utilised.
The four charged leptons are required to fulfil where and can be any type of leptons and + and − have to be oppositely charged leptons regardless of flavour. In addition, a cut on the invariant mass of the leptonic decay products of the Z bosons is applied Note that this cut also removes the phase-space region containing the Higgs resonance. After jet clustering, jets that fulfil the conditions are called identified jets, and at least two of them are required. Note that the last condition demands a minimal distance between the jet and all four charged leptons. The two identified jets with highest transverse momenta, called hardest or leading jets, should obey

Cross sections
In this section, numerical results are discussed for the fiducial cross section in the setup defined above (inclusive setup for short) as well as corresponding results with a stronger cut on the invariant mass of the two leading jets (VBS setup) In the experimental analysis, the first three contributions are often referred to as signal, interference, and QCD background, respectively. We note that the contributions of bottom quarks, which are not included in the results of table 1, amount to 0.2%, −1.1% and 2.8% of the three contributions, respectively. The most striking fact is the size of the QCD contributions with respect to the EW ones. In our default setup, the QCD component is an order of magnitude larger, and the EW component containing VBS contributions is only about 0.1 fb. This is in contrast to same-sign-WW and WZ scattering which have larger cross sections as well as larger signal-to-background ratios [17,18]. Imposing the JHEP11(2020)110 (1)  stronger cut M j 1 j 2 > 500 GeV, the QCD background is suppressed and is only twice as large as the EW contribution, which is only moderately decreased. 3 The relevance of the interference grows from 1% to 2.5% of the cross section, while the relative contribution of the loop-induced process decreases from 10% to 6%.
Including 7-point scale variations as defined above, the EW contribution to the LO cross section of order O α 6 reads for M j 1 j 2 > 500 GeV. (3.11) In table 2 we present results for the NLO cross section including EW corrections of order O α 7 , QCD corrections of order O α s α 6 , or both as well as the corresponding corrections normalised to the LO cross section of order O α 6 in percent. We also list the numbers corresponding to the maxima and minima of the 7-point scale variation in absolute terms and relative to the central values in percent. The scale dependence is reduced by a factor of 2 for the inclusive setup and even more in the setup with the additional VBS cut when including the QCD corrections of order O α s α 6 . We find negative EW corrections at the level of 16%-17%, i.e. in the same range as for other VBS processes [16,18]. While the relative EW corrections are only marginally affected by the additional cut, the QCD corrections are drastically reduced. The reason for this sizeable effect is discussed further below.
We turn to the discussion of different partonic channels. To this end, we split all partonic channels into 4 subsets: VBS-WW encompasses the 16 partonic channels that contain WW → ZZ as subprocess, VBS-ZZ is made up of the remaining 32 channels that include the ZZ → ZZ subprocess. The left-over channels are further separated into 4 that contain pp → WZZ as subprocess (WZZ) and 8 that then always include the pp → ZZZ subprocess (ZZZ). We note that in total 36 partonic channels involve ZZ → ZZ, 8 involve WZZ, and 16 involve ZZZ. None of the channels involves both WW → ZZ and WZZ.

JHEP11(2020)110
Order The contributions of these different partonic processes are compiled in table 3, where we show the corresponding contributions of orders O α 6 , O α 7 , and O α s α 6 in ab, as well as the NLO corrections in percent. The LO O α 6 cross section is dominated by the 16 partonic channels containing WW → ZZ as subprocess. The remaining partonic channels contribute about 2.5% and 1.0% in the inclusive and VBS setup, respectively, at LO and similarly at the order O α 7 . The relative EW corrections are smaller for the non-VBS-WW channels than for the VBS-WW channels apart from ZZZ in the VBS setup, which is however very small. The O α s α 6 contributions, on the other hand, are dominated by channels involving triple-vector-boson production in the inclusive setup. In the inclusive setup more than 70% of the VBS-ZZ contribution in the fifth column results from partonic channels that also involve WZZ. Note that at this order also gq channels contribute at the same level as the qq channels and are included in columns 5 and 6 of table 3. In the VBS-setup, the VBS channels and the non-VBS channels practically cancel at order O α s α 6 . The cut M j 1 j 2 > 500 GeV reduces the O α s α 6 contributions of the WZZ/ZZZ channels by almost an order of magnitude. Note that the QCD corrections are small for the dominating VBS-WW channels, but huge for the WZZ/ZZZ channels. The huge QCD corrections result from contributions with three resonant vector bosons that are present at NLO QCD but not at LO (see below). 4 4 The huge QCD corrections raise the question about the relevance of the non-trivial quark-mixing matrix.
Out of the 24% QCD corrections, 14% result from partonic channels with (anti-)quark-gluon in the initial states. For these channels, owing to the unitarity of the CKM matrix, the effects of a non-trivial quarkmixing matrix only result from mixing of the first two generations with the top quark, which is very small. The leading effect of quark-mixing results from the quark-antiquark s-channel contributions and is of order JHEP11(2020)110  In table 4 we show the relative EW corrections of order O α 7 for selected partonic processes. In the LO cross sections σ α 6 the contributions of qq and q q are combined and the channels resulting from interchanging all quarks of the first and second generations are included. The last column shows the subprocesses present in the channels. We have selected those channels with the relevant subprocesses that give the largest contributions. While the relative EW corrections for the channels involving VBS-WW are in the range 15%-18%, they are between 7% and 19% for the other channels.
The large EW corrections for the dominant VBS-WW channels can be explained based on a Sudakov approximation applied to the WW → ZZ subprocesses, as was already shown for like-sign WW scattering [16] and WZ scattering [18]. Following refs. [60,61] one can show that the leading logarithmic corrections to the scattering of transverse vector bosons, which is the dominant contribution, yields the simple correction factor It includes all logarithmically enhanced EW corrections apart from the angular-dependent subleading soft-collinear logarithms and applies to all VBS processes that are not mass |Vus/V ud | 2 ∼ |V cd /Vcs| 2 ∼ 5%. Since these channels cause 10% QCD corrections, the effect of neglecting quark mixing is at the level of 10% × 5% = 0.5%.
suppressed, such as WW → ZZ, owing to the fact that these scattering processes result from the same SU(2) w coupling. Equation (3.12) is, however, not valid for mass-suppressed processes like ZZ → ZZ. The constants are given by C EW W = 2/s 2 w and b EW W = 19/(6s 2 w ), where s w represents the sine of the weak mixing angle. Further, Q is a representative scale of the V V → V V scattering process, which can conveniently be chosen as the four-lepton invariant mass M 4 . Using Q = M 4 event by event, results in the numbers for δ LL shown in the 4th column of table 4, which agree within 2% with the exact NLO results. It should be clear that the formula (3.12) is not applicable to the non-VBS-WW processes. For the angular-dependent subleading soft-collinear logarithms, correction factors can be derived as well in the Sudakov limit based on the results of refs. [60,61]. These depend on the specific VBS process, and the correction factor for WW → ZZ reads where s 12 , s 13 , and s 23 are the Mandelstam variables of the VBS process. 5 Applying the correction factor (3.13) combined with (3.12) event by event yields the approximations for the EW corrections in column 5 of table 4. While the inclusion of these angular-dependent logarithmic terms somewhat deteriorates the agreement with the full correction factors, it shows that the effect of the angular-dependent logarithmic corrections is only at the level of 2% and thus well within the accuracy of the approximation, which is expected to be a few percent. 5 To determine the Mandelstam variables, we need the momenta of the scattering vector bosons. The momentum of one of them is fixed by combining the momentum p1 of one of the incoming partons with the momentum pparton,out of the outgoing parton that maximises (p1 − pparton,out) 2 . The momentum of the second scattering vector boson is obtained from the momentum of the second incoming parton and the other outgoing parton.

Differential distributions
In this section we discuss distributions for the inclusive setup (3.9) only.
LO distributions. First, we display distributions at LO in figure 5 including theoretical predictions at orders O α 6 , O α s α 5 , and O α 2 s α 4 . In addition, we include predictions for the loop-induced contributions at the order O α 4 s α 4 . While the upper panels show the absolute predictions, the lower ones display each contribution relative to the sum of all of them.
The rapidity difference between the two hardest jets, shown in figure 5(a), is a typical handle to enhance the EW contribution. As expected, the EW contribution becomes dominant only in a rather extreme part of the phase space, typically for |∆y j 1 j 2 | > 5. The central region is largely dominated by the QCD contributions which amount to more than 90% of the total process. The interference contribution is below 1% and hardly visible in the figure. The loop-induced process is particularly interesting. It is minimal for low rapidity difference but relatively increasing towards large rapidity difference like the EW contribution.
A further variable that is used to enhance the EW contribution is the di-jet invariant mass, whose distribution is shown in figure 5(b). The EW contribution exceeds the QCD contribution only for M j 1 j 2 > 1200 GeV. The fraction of the loop-induced contribution decreases with increasing invariant mass, while the interference increases but stays below 5% for M j 1 j 2 < 2000 GeV.
In figure 5(c), the distribution in the transverse momentum of the hardest jet is shown. The relative EW contribution is slowly increasing with p T,j 1 to reach almost 20% at 400 GeV. The interference contribution becomes non negligible in this part of phase space and amounts to about 10% for p T,j 1 = 800 GeV. Finally, the loop-induced contribution drops quickly and is below 3% above 200 GeV.
The distribution in the transverse momentum of the four leptons, shown in figure 5(d), behaves qualitatively similar to the distribution in the transverse momentum of the hardest jet. This is expected as these two observables are correlated. It is worth noticing that the interference contribution does not increase towards high transverse momentum as in the previous case and is thus almost imperceptible over the whole range. Also, the loop-induced contribution is dropping less quickly than in the previous case. In particular, in the first bin of the distribution in the transverse momentum of the hardest jet, the loop-induced process represents about 25% of the total predictions, while here it is about 15%. The observables shown in figure 6 are related to the two hardest jets. We start with the two observables that are typically used to enhance EW contributions over its irreducible JHEP11(2020)110 QCD background: the invariant mass ( figure 6(a)) and the rapidity difference ( figure 6(b)) of the two leading jets. The most interesting feature is that the bulk of the positive QCD corrections is located at low di-jet invariant mass. The O α s α 6 corrections are well above 100% in the first bins but become of the order of a few percent above 600 GeV. This effect has been already observed and partially discussed in ref. [44] for like-sign W-boson scattering and is due the appearance of W or Z bosons becoming resonant and decaying to a pair of jets when including real gluon radiation. As shown in figures 1(g) and 1(h), JHEP11(2020)110 pp → e + e − µ + µ − jj also includes contributions from triple vector-boson production (WZZ and ZZZ) at order O α 6 . At LO, the massive vector boson decaying hadronically cannot become resonant due to the cut on the invariant mass of the two jets at 100 GeV [see eq. (3.9)]. When including real radiation at NLO, the cut (3.9) does not necessarily apply to the two quarks coming from the vector-boson decay. An extra gluon jet can make up the M j 1 j 2 > 100 GeV and allows the two quarks to originate from a resonant W or Z boson. Thus, the relatively large QCD corrections found here are a combined effect of the JHEP11(2020)110 real radiation, the event selection, and the inclusion of tri-boson contributions. We would like to emphasise that the phase-space region of our default setup is rather inclusive and should be avoided if one does not include tri-boson contributions in the theory predictions. Otherwise, the fiducial corrections at NLO QCD will most likely be off by about 20%, and this offset is by no means accounted for by scale-variation uncertainties. This consideration is particularly important if Monte Carlo programs using the VBS approximation [44] are used to extrapolate measurements from the inclusive region to more VBS-enriched regions. An alternative is to subtract the on-shell tri-boson contributions from the results for the full process. Such a strategy is often used in experimental analysis to extract VBS contributions. While it allows one to claim a VBS measurement, care has to be taken not to violate gauge invariance. Furthermore, it has the disadvantage to make the measurement even more theory dependent. In our opinion, the most straightforward and physical measurement would include both the QCD and the EW processes and in the latter all possible contributions including tri-boson production, i.e. all contributions to a given physical final state. The EW corrections, on the other hand, do not display an unexpected behaviour but confirm the results known from other VBS signatures [16][17][18]. They become negatively large for large invariant masses owing to enhanced EW logarithms to reach −20% at 2 TeV. Turning to the distribution in the rapidity difference shown in figure 6(b), the QCD corrections reach almost 300% in the central rapidity region. The rapidity separation of the two hardest jets is strongly correlated to their invariant mass (see, for instance, figure 3 of ref. [44]). Thus, the arguments given for the distribution in M j 1 j 2 can be transfered to the distribution in ∆y j 1 j 2 . Events with small ∆y j 1 j 2 are depleted at LO owing to the cut (3.9), while this is not the case at NLO QCD where extra gluons can provide a leading jet. The distribution also shows that a cut on the rapidity difference would be very effective in removing the sizeable QCD corrections linked to triple-vector-boson production in a similar way as a stronger cut on M j 1 j 2 . Thus, the large QCD corrections could be reduced by either a cut on ∆y j 1 j 2 or a stronger one on M j 1 j 2 , which are usually imposed in VBS studies. The EW corrections are moderate and vary between −10% for zero rapidity difference and −20% for large rapidity differences.
Figures 6(c) and 6(d) show distributions in the azimuthal-angle difference and the cosine of the angle between the two leading jets, respectively, which provide information on the correlation between the two jets. The EW corrections are rather stable throughout the kinematic range and vary by less than 7%. For the azimuthal-angle difference, the QCD corrections are maximal near 30 • where they are about 55%. When the two jets have maximal azimuthal-angle difference, the LO contribution is maximal and receives QCD corrections at the level of 10%. The distribution in cos θ j 1 j 2 peaks at −1, i.e. when the two jets are back-to-back. The QCD corrections are minimal there but exceed 200% when the two jets are close to each other.
In figure 7 we display distributions related to the 4-lepton system, i.e. the Z-boson pair, and the electron-positron pair, i.e. one of the Z bosons. The distribution in the invariant mass of the 4-lepton system receives the typical EW Sudakov corrections growing to −40% at M 4 = 2 TeV. The QCD corrections increase from 20% to 40% in the considered range. As a consequence, the overall NLO corrections experience a cancellation for larger invariant masses and decrease from 20% to almost zero towards large

JHEP11(2020)110
While the QCD corrections hardly modify the shape of this distribution, the EW ones exhibit a radiative tail below the resonance reaching more than +20%. This tail is due to real photon radiation that takes away part of the energy of the electron-positron system and thus shifts events from the peak towards lower invariant masses. Such an effect is well known and has been observed already for Drell-Yan, di-boson, or top-anti-top production processes. The distribution in the transverse momentum of the electron-positron system, i.e. the transverse momentum of one of the Z bosons, is presented in figure 7(d). The EW corrections increase from −10% to −40% and the QCD corrections from 25% to 40%, leading to an overall NLO correction decreasing from 15% to almost zero. Next we study distributions in transverse momenta and rapidities of the leading jet and the positron in figure 8. As observed in previous computations of VBS processes [17,18], the distribution in the transverse momentum of the leading jet ( figure 8(a)) is suppressed for small transverse momenta and receives large QCD corrections in this region. This can be attributed to the presence of an extra jet and reshuffling of energy between the jets. In the rest of the spectrum, the QCD corrections are at the level of +20% as for the fiducial cross section. The EW corrections show the typical high-energy behaviour of other transverse-momentum distributions growing negatively large with increasing transverse momentum and reaching about −40% at 800 GeV. The distribution in the rapidity of the hardest jet ( figure 8(b)) displays a similar behaviour as the rapidity difference of the jet pair. Finally, we show distributions in angular variables related to pairs of leptons in figure 9. Such observables are particularly useful in ZZ scattering as they give insight in details of the scattering process which is not possible for other VBS signatures due to the presence of neutrino(s). We start with distributions in azimuthal-angle differences for the electron-positron system (figure 9(a)) and the positron-muon system ( figure 9(b)). Given the flavour, the first observable relates two opposite-sign leptons originating from the same Z-boson decay, while the second one links leptons of two different Z-boson decays. The shapes of the two distributions are therefore rather different. The distribution for the electron-positron system peaks near 50 • , while the one for the positron-muon case is maximal at 180 • . The corrections behave also rather differently. In the electron-positron case, the QCD (EW) corrections are maximally positive (negative) at low angle, leading to an overall correction which is rather stable around +10% over the whole spectrum. For the positron-muon azimuthal separation, the QCD corrections slightly increase by almost ten percent from low to large angles, while the EW ones have the opposite trend. Thus, the overall corrections are again rather steady at about 10%. We emphasise that such cancellations are, a priori, largely accidental and should not be taken for granted in other VBS signatures or setups. The last two distributions concern the cosine of the angle between the two leptons for the same two leptonic systems (figures 9(c) and 9(d)). It is interesting to notice that for these distributions, again both types of corrections do not induce large shape distortions (less than 10%). This makes such observables particularly attractive for correlation analysis or to test models of new physics.

Conclusion
In this article we have presented a calculation of the NLO EW and QCD corrections of orders O α 7 and O α s α 6 , respectively, for the process pp → e + e − µ + µ − jj + X. Our results include, in particular, the NLO EW and QCD corrections to the LO EW contribution of order O α 6 , which is dominated by vector-boson scattering (VBS) into a pair of Z bosons. As the full matrix elements are used at the corresponding orders, our computation accounts for all off-shell, non-resonant, and interference effects. In particular,

JHEP11(2020)110
triple-boson-production processes are part of the EW contributions. We have included all partonic channels apart from those involving initial-state bottom quarks or photons, which are suppressed. The EW corrections of order O α 7 are found to be relatively large, in agreement with similar results obtained previously for W ± W ± and WZ scattering and the expectation that large EW corrections are an intrinsic feature of VBS at the LHC. For the chosen fiducial cross section, the corrections are −16% and can be well reproduced by a simple logarithmic approximation. In the high-energy tails of distributions the EW corrections can reach −40%. The QCD corrections of order O α s α 6 exceed 20% for the fiducial cross section. While such a magnitude for QCD corrections at the LHC is perfectly normal, it is rather large for VBS processes. In particular, previous computations relying on the VBS approximation have found QCD corrections at the percent level, in strong contrast to our findings. These differences are due to the fact that our computation is done for a rather inclusive phase-space region (M j 1 j 2 > 100 GeV as opposed to M j 1 j 2 500 GeV and/or ∆y j 1 j 2 3 usually) and includes tri-boson contributions. Hence, when including real QCD radiation at NLO, Z or W bosons decaying hadronically can become resonant thus leading to very large corrections that are not captured by scale variation. We emphasise that there is nothing wrong in using such inclusive fiducial regions as long as the theoretical simulations used include tri-boson contributions at NLO QCD. In particular, Monte Carlo programs relying on the VBS approximation should not be used to extrapolate measurements from the inclusive region to more VBS-enriched regions. Nonetheless, they still constitute reliable theoretical predictions in the typical VBS regions (with M j 1 j 2 500 GeV and/or ∆y j 1 j 2 3), where the NLO QCD corrections are of the usual size. The present work also shows that the proper inclusion of tri-boson contributions is of great relevance. We strongly advocate for simpler and more physical measurements where tri-boson contributions are not subtracted from the signal. This has the advantage to be clearly gauge invariant and not to rely on conventions and theory predictions. Finally, as argued in previous work, comparing predictions including all QCD, interference, and EW contributions with the measurements is the cleanest way of testing the SM in such processes.
We have also analysed the composition of the LO process by comparing contributions at orders: O α 6 (EW contribution), O α s α 5 (interference), and O α 2 s α 4 (QCD contribution). The findings are in line with known observations that VBS ZZ is a challenging channel due to the very large irreducible QCD background. In particular, in the fiducial region chosen, the LO EW contributions amount to less than 10%. In the LO analysis, we have also further included the loop-induced contribution with gluon-gluon initial state at order O α 4 s α 4 . It has been found to be of the order of 10% in the fiducial region, but becomes relatively negligible in the high-energy limit of differential distributions. Imposing an extra VBS cut M j 1 j 2 > 500 GeV enhances the EW contribution to more than 30%.
We would like to point out that this calculation constitutes a further leap in complexity with respect to previous calculations for VBS processes. The fact that eight charged external particles occur increases significantly the complexity of the virtual and, in particular, of the real corrections. The number of partonic channels is amplified by a factor of 5 with respect to like-sign W scattering and 1.5 relative to WZ scattering. As a consequence, the JHEP11(2020)110 required CPU time is increased and efficient book-keeping, automation and parallelisation are crucial.
Finally, we hope that these results will be useful in the current and upcoming measurements of ZZ VBS at the LHC. In particular, we believe that this article contains valuable information for the experimental collaborations when conducting their analysis.