Evidence for electroweak production of four charged leptons and two jets in proton-proton collisions at $\sqrt{s} =$ 13 TeV

Evidence is presented for the electroweak (EW) production of two jets (jj) in association with two Z bosons and constraints on anomalous quartic gauge couplings are set. The analysis is based on a data sample of proton-proton collisions at $\sqrt{s}= $ 13 TeV collected with the CMS detector in 2016-2018, and corresponding to an integrated luminosity of 137 fb$^{-1}$. The search is performed in the fully leptonic final state ZZ$\to\ell\ell\ell'\ell'$, where $\ell,\ell' = $ e, $\mu$. The EW production of two jets in association with two Z bosons is measured with an observed (expected) significance of 4.0 (3.5) standard deviations. The cross sections for the EW production are measured in three fiducial volumes and the result is $\sigma_{\mathrm{EW}}$(pp $\to$ ZZjj $\to\ell\ell\ell'\ell'$jj) = 0.33$^{+0.11}_{-0.10}$ (stat) $^{+0.04}_{-0.03}$ (syst) fb in the most inclusive volume, in agreement with the standard model prediction of 0.275 $\pm$ 0.021 fb. Limits on anomalous quartic gauge couplings are derived in terms of the effective field theory operators T0, T1, T2, T8, and T9.


Introduction
In the standard model (SM), the electroweak (EW) vector bosons, like the other fundamental particles, acquire their masses through the coupling to the Brout-Englert-Higgs field. The photon remains massless, with only two degrees of polarization (i.e., transverse), whereas the W and Z bosons acquire an additional degree of freedom (i.e., longitudinal), as a consequence of the electroweak symmetry breaking (EWSB) [1,2]. Thus, the scattering of massive vector bosons is at the heart of the EWSB mechanism and its study can lead to significant insight into the origin of particle masses. Moreover, if the couplings between the Higgs boson and vector bosons (HVV) differ from their SM values, the subtle interplay between HVV, triple, and quartic gauge couplings as predicted in the SM is incomplete, and the cross section for the longitudinal scattering diverges at large scattering energies, eventually violating the unitarity.
At the CERN LHC, vector boson scattering (VBS) is the interaction of two EW vector bosons emitted by quarks (q) from the two colliding protons. The VBS process is generally labeled by the type of outgoing vector bosons. The two jets (jj) originating from the scattered quarks are typically emitted in the forward-backward region of the detector, giving rise to events whose signature in the detector is characterized by a region in rapidity (so-called "rapidity gap") [3,4], where no additional hadronic activity is expected from the hard scattering. The decay of the vector bosons into fermions defines the final signature of the VBS-like event. The pure VBS contributions, however, are embedded into a wider set of possible two-to-six processes, with which they interfere (Fig. 1). All processes at the order of α 6 EW (tree level) are considered as EW production ( Fig. 1 upper panels and bottom left panel), whereas the processes at the order α 4 EW α 2 S where at tree level the jets are induced by quantum chromodynamics (QCD) (lower right panel in Fig. 1), constitute a background referred to as QCD-induced background. Kinematic requirements on the dijet system are used to define fiducial regions enriched in VBS-like events and where QCD-induced backgrounds are suppressed.
Both the ATLAS and CMS Collaborations have performed searches for the scattering of massive vector bosons, using data from proton-proton (pp) collisions at the center-of-mass energy of 13 TeV. The ATLAS Collaboration reported the observation of EW production of two jets in association with a same-sign W boson pair [5], with a WZ boson pair [6], and, recently, with a Z boson pair [7]. Results were also reported on the measurement of the EW diboson production (WW, WZ, ZZ) in association with a high-mass dijet system in semileptonic final states [8], with an observed significance of 2.7 standard deviations. The CMS Collaboration observed the production of two EW-induced jets with two same-sign W bosons [9, 10] and with WZ pairs [10], and measured the EW production of jets in association with ZZ [11] with an observed significance of 2.7 standard deviations. This paper presents evidence for the EW production of two jets in association with two Z bosons, where both Z bosons decay into electrons or muons, ZZ → ( , = e, µ). Despite a low cross section, a small Z → branching fraction, and a large QCD-induced background, this channel provides a clean leptonic final state with a small experimental background, where one or more reconstructed lepton candidates originate from the misidentification of jet fragments or from nonprompt leptons.
The search for the EW-induced production of the jj final state is carried out using p p collisions at √ s = 13 TeV recorded with the CMS detector at the LHC. The data set corresponds to an integrated luminosity of 137 fb −1 collected in 2016, 2017, and 2018. A discriminant based on a matrix element likelihood approach (MELA) [12-16] is used to extract the signal significance and to measure the cross sections for the EW and the EW+QCD production of the jj final state in a fiducial volume. Finally, the selected jj events are used to constrain anomalous quartic gauge couplings (aQGC) described in the effective field theory approach [17] by the operators T0, T1, and T2, as well as the neutral-current operators T8 and T9 [18]. signal efficiency and acceptance, and to model the signal and irreducible background contributions in the signal extraction fit.
The EW production of two Z bosons and two final-state quarks, where the Z bosons decay leptonically, is simulated at leading order (LO) using MADGRAPH5 aMC@NLO v2.4.6 (abbreviated as MG5 in the following) [21]. The leptonic Z boson decays are simulated using MADSPIN [22]. The sample includes triboson processes, where the Z boson pair is accompanied by a third vector boson that decays hadronically, as well as diagrams involving the quartic gauge coupling vertex. The predictions from this sample are cross-checked with those obtained from the LO generator PHANTOM v1.2.8 [23] with agreement in the yields and the distributions exploited for the signal extraction.
The leading QCD-induced production of two Z bosons in association with jets, whose contribution with two jets in the final state is referred to as qq → ZZjj, is simulated at next-to-leading order (NLO) with MG5 with up to two extra parton emissions, and merged with the parton shower simulation using the FxFx scheme [24]. Next-to-next-to-leading order corrections calculated with MATRIX v1.0.0 [25][26][27] are applied as K factors, differentially as a function of the invariant mass of the ZZ system (m ZZ ). Additional NLO EW corrections are applied for m ZZ > 2m Z , following the calculations from Ref. [28]. The interference between the EW and QCD diagrams is evaluated using dedicated samples produced with MG5 at LO. Its contribution is added in the cross section fits via a linearized scaling law. The loop-induced production of two Z bosons from a gluon-gluon (gg) initial state, whose contribution with two jets in the final state is referred to as gg → ZZjj, is simulated at LO with up to two extra parton emissions using MG5 by explicitly requiring a loop-induced process. For the 1-and 2-jet contributions, a pp initial state instead of gg is specified in MG5 to also include initial-state radiation contributions where a gluon involved in the hard process is emitted from an initial quark. Finally, the samples with 0 to 2 extra partons are merged with parton shower simulation using the MLM matching scheme [29,30]. An NLO/LO K factor, which is extracted from Refs. [31,32], is used to normalize this process.
Background processes that contain four prompt, isolated leptons and additional jets in the final state, namely ttZ and VVZ (V = W, Z), are simulated with MG5 at NLO.
The simulation of the aQGC processes is performed at LO using MG5 and employs matrix element reweighting to obtain a finely spaced grid for each of the five anomalous couplings probed by the analysis.
The PYTHIA 8.226 and 8.230 [33] package versions are used for parton showering, hadronization and the underlying event simulation, with parameters set by the CUETP8M1 tune The detector response is simulated using a detailed description of the CMS detector implemented in the GEANT4 package [37,38]. The simulated events are reconstructed using the same algorithms used for the data, and include additional interactions in the same and neighboring bunch crossings, referred to as pileup. Simulated events are weighted so that the pileup distribution reproduces that observed in the data, which has an average of about 23 (32) interactions per bunch crossing in 2016 (2017 and 2018).

Event reconstruction and selection
The final state consists of at least two pairs of oppositely charged isolated leptons and at least two hadronic jets. The ZZ selection is similar to that used in the CMS H → ZZ → measurement [39].
The primary triggers require the presence of a pair of loosely isolated leptons, whose exact requirements depend on the data-taking year. Triggers requiring three leptons with low transverse momentum (p T ), as well as isolated single-electron and single-muon triggers, help to recover efficiency. The overall trigger efficiency for events that satisfy the ZZ selection described below is > 98%.
Events are reconstructed using a particle-flow algorithm [40] that identifies each individual particle with an optimized combination of all subdetector information. The candidate vertex with the largest value of summed physics-object p 2 T is the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [41, 42] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum (p miss T ), taken as the negative vector sum of the p T of those jets (which include the leptons). Electrons are identified using a multivariate classifier, which includes observables sensitive to bremsstrahlung along the electron trajectory, the geometrical and energy-momentum compatibility between the electron track and the associated energy cluster in the electromagnetic calorimeter, the shape of the electromagnetic shower, isolation variables, and variables that discriminate against electrons originating from photon conversions [43].
Muons are reconstructed by combining information from the silicon tracker and the muon system [44]. The matching between the muon-system and tracker tracks proceeds either outsidein, starting from a track in the muon system, or inside-out, starting from a track in the silicon tracker. The muons are selected from the reconstructed muon track candidates by applying minimal requirements on the track in both the muon system and silicon tracker.
To further suppress electrons from photon conversions and muons originating from in-flight decays of hadrons, the three-dimensional impact parameter of each lepton track, computed with respect to the primary vertex position, is required to be less than four times the uncertainty in the impact parameter.
Leptons are required to be isolated from other particles in the event. The relative isolation is defined as where the sums run over the charged and neutral hadrons, as well as the photons, in a cone defined by ∆R ≡ (∆η) 2 + (∆φ) 2 = 0.3 around the lepton trajectory, where η and φ denote the azimuthal angle and pseudorapidity of the particle, respectively. To minimize the contribution of charged particles from pileup to the isolation calculation, charged hadrons are included only if they originate from the primary vertex. The contribution of neutral particles from pileup p PU T is evaluated for electrons with the jet area method described in Ref. [45]. For muons, p PU T is taken as half the p T sum of all charged particles in the cone originating from pileup vertices. The factor of one-half accounts for the expected ratio of charged to neutral particle production in hadronic interactions. Muons with R iso < 0.35 are considered isolated, whereas for electrons, the R iso variable is included in the multivariate classifier.
The lepton reconstruction and selection efficiency is measured in bins of p T and η using the tag-and-probe technique [46] on events with single Z bosons. The measured efficiencies are used to correct the simulation. The muon (electron) momentum scales are calibrated in bins of p T and η using the J/ψ meson and Z boson (Z boson only) leptonic decays. Jets are reconstructed from particle-flow candidates using the anti-k T clustering algorithm [41], as implemented in the FASTJET package [42], with a distance parameter of 0.4. To ensure a good reconstruction efficiency and to reduce the instrumental background, as well as the contamination from pileup, loose identification criteria based on the multiplicities and energy fractions carried by charged and neutral hadrons are imposed on jets [47]. Only jets with |η| < 4.7 are considered.
Jet energy corrections are extracted from data and simulated events to account for the effects of pileup, uniformity of the detector response, and residual differences between the jet energy scale in data and simulation. The jet energy scale calibration [48, 49] relies on corrections parameterized in terms of the uncorrected p T and η of the jet, and is applied as a multiplicative factor, scaling the four-momentum vector of each jet. To ensure that jets are well measured and to reduce the pileup contamination, all jets must have a corrected p T > 30 GeV. Jets from pileup are further rejected using pileup jet identification criteria based on the compatibility of the associated tracks with the primary vertex inside the tracker acceptance and on the topology of the jet shape in the forward region [50].
A signal event must contain at least two Z candidates, each formed from pairs of isolated electrons or muons of opposite charges. Only reconstructed electrons (muons) with p T > 7 (5) GeV are considered. At least two leptons are required to have p T > 10 GeV and at least one is required to have p T > 20 GeV. All leptons are required to be separated by ∆R ( 1 , 2 ) > 0.02, and electrons are required to be separated from muons by ∆R (e, µ) > 0.05.
Within each event, all permutations of leptons giving a valid pair of Z candidates are considered. For each ZZ candidate, the lepton pair with the invariant mass closest to the nominal Z boson mass is denoted Z 1 . The other dilepton candidate is denoted Z 2 . Both m Z 1 and m Z 2 are required to be in the range 60-120 GeV. All pairs of oppositely charged leptons that can be built from the ZZ candidate, regardless of flavor, are required to satisfy m > 4 GeV to suppress backgrounds from hadron decays. If multiple ZZ candidates in an event pass this selection, the one with the largest scalar p T sum of the Z 2 leptons is retained. Finally, the invariant mass of the four leptons is required to satisfy m 4 > 180 GeV. This selection is referred to as the ZZ selection.
The search for the EW production of two Z bosons is performed on a subset of events that pass the ZZ selection, namely those with at least two jets. The jets are required to be separated from the leptons of the ZZ candidate by ∆R > 0.4. The two highest p T jets are referred to as the tagging jets and their invariant mass (m jj ) is required to be > 100 GeV. This selection is referred to as the ZZjj inclusive selection and is used to measure the signal significance, the total fiducial cross sections, and to perform the aQGC search. Additionally, two VBS signal subregions are defined for fiducial cross section measurements in signal-enriched regions: a loose VBS signalenriched region that requires m jj > 400 GeV and |∆η jj | > 2.4 and corresponds to a signal purity of ≈20%, and a tight VBS signal-enriched region that requires m jj > 1 TeV and |∆η jj | > 2.4 and corresponds to a signal purity of ≈50%. Finally, a background control region is defined from events that satisfy the ZZjj inclusive selection but fail at least one of the criteria that define the loose VBS signal-enriched region.

5 Background estimation
The dominant background arises from the production of two Z bosons in association with QCD-induced jets. The yield and shape of the matrix element discriminant for this irreducible background are taken from simulation, but ultimately constrained by the data in the fit that extracts the EW signal, as described in Section 7. Other irreducible backgrounds arise from processes that produce four genuine high-p T isolated leptons, pp → ttZ+jets and pp → VVZ+jets. These small contributions feature kinematic distributions similar to that of the dominant background and are estimated from simulation.
Reducible backgrounds arise from processes in which heavy-flavor jets produce secondary leptons or from processes in which jets are misidentified as leptons. They are referred to as Z+X and are predominately composed of Z+jets events, with minor contributions from tt +jets and WZ+jets processes. The lepton identification and isolation requirements significantly suppress this background, which is only 2-3% after the ZZjj inclusive selection and is even smaller in the signal region. This reducible contribution is estimated from data by weighting events from a control region by a lepton misidentification rate, which is also determined from data. Events in the control region satisfy the ZZjj inclusive selection, with the exception that the Z 2 is composed of same-sign same-flavor leptons (SS-SF). The SS-SF leptons are required to originate from the primary vertex without any identification or isolation requirement.
The lepton misidentification rate is measured by selecting events that feature one Z boson candidate and a third reconstructed lepton. The fraction of events for which the third lepton satisfies the identification and isolation criteria is the lepton misidentification rate. The misidentification rates are evaluated using the tight requirement |m Z 1 − m Z | < 7 GeV to reduce the contribution from asymmetric photon conversions, and p miss T < 25 GeV to suppress the WZ contribution. The procedure is identical to that used in Ref. [39].

Systematic uncertainties
The uncertainties in the QCD renormalization and factorization scales for the signal and in the jet energy scale are the two dominant systematic uncertainties in the measurement. The impact of the variation for each source of uncertainty is summarized below.
Renormalization and factorization scale uncertainties are evaluated by varying both scales independently. The following variations from the default scale choice , taking the largest variation as the systematic uncertainty, which is about 6% for the EW signal and ranges from 10 to 12% for the qq → ZZjj QCD background. All quoted ranges correspond to variations for the different final states. Since the uncertainty in gg → ZZjj that relates to missing higher order corrections are accounted for using a K factor, an uncertainty in the normalization of 11% is used, as derived from Refs. [31,32]. The PDF and related α S variations are evaluated from the variations of the respective eigenvalues set following the NNPDF prescription [36], and are 3.2% (6.6%) for the qq → ZZjj QCD background (EW signal). Although the PDFs used are different in the various years (see Section 3), the associated uncertainties are very similar. Given the small dependence on the discriminant value, a constant value of 3-6% is used for these uncertainties, depending on the sample considered.
The impact of the jet energy scale uncertainty ranges from 4.9 to 11.4% (0.7 to 1.2%) for the qq → ZZjj QCD background (EW signal) and the impact of the jet energy resolution uncertainty [49] is 2.2-6.3% (0.2-0.4%). The uncertainty in the trigger as well as the lepton recon-struction and selection efficiencies ranges from 2.5-9%. The uncertainty in the integrated luminosity is 2.3-2.5% depending on the data-taking period [51][52][53]. The uncertainty in the estimate of the reducible background from control samples ranges from 33% to 45%, depending on the final state. This uncertainty includes the limited number of events in the control regions as well as differences in background composition between the control regions used to determine the lepton misidentification rates and those used to estimate the yield in the signal region. The uncertainty from the limited size of the MC samples amounts to 2.5-4.2% for the qq → ZZjj QCD background, 3.2% for the gg → ZZjj QCD background, and is < 1% for the EW signal. For ttZ and VVZ, the limited MC sample size is the dominant source of uncertainty, ranging from 19 to 24%.

Search for the EW production of ZZ with two jets
After the ZZjj inclusive selection, the expected EW signal purity is about 6% with 85% of events coming from the QCD-induced production. Additional kinematic selections are therefore necessary to enhance the contribution from EW production. Table 1 presents the expected and observed event yields for the ZZjj inclusive selection, as well as for the loose and tight VBS signal-enriched selections. The determination of the signal strength for the EW production, i.e., the ratio of the measured cross section to the SM expectation µ = σ/σ SM , utilizes a matrix element discriminant (K D ) to separate the signal and the QCD background. The discriminant is constructed following the approach described in Refs. [13][14][15]. The performance of the K D discriminant was checked against a multivariate discriminant based on a boosted decision tree (BDT) employing seven input variables (m jj , ∆η jj , m 4 , η * , R(p jets T )) as defined and used in Ref.
[11]. Furthermore, a BDT using up to 28 input variables, including the above as well as those used in Ref.
[7], was studied and no significant gain was obtained. This confirms that the K D discriminant captures the differences between the kinematical distributions of signal and background events. Figure 2 presents the m jj and |∆η jj | distributions in the ZZjj inclusive region. The distribution of the K D discriminant for all events in the ZZjj inclusive selection is shown in Fig. 3. The high signal purity contribution is visible at large discriminant values.
The distribution of the K D discriminant for the backgrounds is validated in the background control region defined by selecting events with m jj < 400 GeV or |∆η jj | < 2.4. A good agreement is observed between the data and the SM expectation.  The K D discriminant distribution for events in the ZZjj inclusive selection is used to extract the significance of the EW signal via a maximum-likelihood fit. The expected distributions for the signal and the irreducible backgrounds are taken from the simulation while the reducible background is estimated from the data. The shape and normalization of each distribution are allowed to vary in the fit within the respective uncertainties. This approach constrains the yield of the QCD-induced production from the background-enriched region of the discriminant distribution.
The systematic uncertainties are treated as nuisance parameters in the fits and profiled [54]. The measured signal strengths from the fits in the three analysis regions are used to determine the fiducial cross sections for the EW and the EW+QCD production. The fiducial volumes are almost identical to the selections imposed at the reconstruction level, and is detailed in Table 2. The generator-level lepton momenta are corrected by adding the momenta of generator-level photons within ∆R( , γ) < 0.1. The kinematic requirements to select Z boson candidates and the final ZZjj candidate are the same as those used for the reconstruction-level analysis. Table 3 reports the measured cross sections and their SM predictions in the three ZZjj fiducial regions. For the SM predictions we report those extracted from generated events in MC samples adopted for the analysis, as well as higher-order calculations at NLO in QCD [55,56]. In addition, we compare with a theoretical prediction at LO in QCD, but including NLO EW corrections; in the ZZjj baseline region it amounts to 0.242 +0.015 −0.013 fb, where the uncertainty comes from variations of the factorization and renormalization scales.
The measured (expected) EW signal strength in the ZZjj inclusive region is µ EW = 1.21 +0.47

Limits on anomalous quartic gauge couplings
In an effective field theory approach to physics beyond the Standard Model, dimension-8 operators stem from covariant derivatives of the Higgs doublet and from charged and neutral field strength tensors associated to gauge bosons. The latter generate eight independent operators, corresponding to couplings of the transverse degrees of freedom (Ti) of the gauge fields. The ZZjj channel is particularly sensitive to the charged-current operators T0, T1, and T2, as well as the neutral-current operators T8 and T9 [18]. The m 4 distribution is used to constrain the aQGC parameters f Ti /Λ 4 , corresponding to Wilson coefficients of the aforementioned operators. Figure 4 shows the expected m 4 distributions with postfit normalizations for the SM and for an example aQGC scenario, as well as the observed distribution in the data. The expected yield enhancement at large values of m 4 exhibits a quadratic dependence on the anomalous couplings, and a parabolic function is fitted to the per-mass bin yields, allowing for an interpolation between the discrete coupling parameters of the simulated aQGC signals. The statistical      Table 4 lists the individual lower and upper limits obtained by setting all other anomalous couplings to zero. The unitarity bounds are determined using the results from Ref.
[60] as the scattering energy m 4 at which the aQGC strength set equal to the observed limit would result in a scattering amplitude that violates unitarity.

Summary
A search was performed for the electroweak production of two jets in association with two Z bosons in the four-lepton final state in proton-proton collisions at 13 TeV. The data correspond to an integrated luminosity of 137 fb −1 collected with the CMS detector at the LHC.
The electroweak production of two jets in association with a pair of Z bosons is measured with an observed (expected) significance of 4.0 (3.5) standard deviations. The measured fiducial cross section is σ fid = 0.33 +0.11 −0.10 (stat) +0.04 −0.03 (syst) fb, which is consistent with the standard model prediction of 0.275 ± 0.021 fb.
Limits on anomalous quartic gauge couplings are set at 95% confidence level in terms of effective field theory operators, with units in TeV −4 : These are the most stringent limits to date on the neutral current operators T8 and T9.

Acknowledgments
We are grateful to A. Denner, R. Franken, M. Pellen and T. Schmidt for providing the calculation of the EW cross section with NLO EW corrections in our ZZjj inclusive fiducial region. We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies:   [7] ATLAS Collaboration, "Observation of electroweak production of two jets and a Z-boson pair with the ATLAS detector at the LHC", (2020). arXiv:2004.10612. Submitted to Nature Phys.
[10] CMS Collaboration, "Measurements of production cross sections of WZ and same-sign WW boson pairs in association with two jets in proton-proton collisions at √ s = 13 TeV", (2020). arXiv:2005.01173. Submitted to PLB.