Evidence for electroweak production of four charged leptons and two jets in proton-proton collisions at

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 √ s = 13 TeV collected with the CMS detector in 2016–2018, and corresponding to an integrated luminosity of 137fb − 1 . The search is performed in the fully leptonic ﬁnal state ZZ → (cid:3)(cid:3)(cid:3) (cid:4) (cid:3) (cid:4) , where (cid:3), (cid:3) (cid:4) = e , μ . The EW production of two jets in association with two Z bosons is measured with an observed (expected) signiﬁcance of 4.0 (3.5) standard deviations. The cross sections for the EW production are measured in three ﬁducial volumes and the result is σ EW ( pp → ZZjj → (cid:3)(cid:3)(cid:3) (cid:4) (cid:3) (cid:4) 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 ± 0 . 021fb. Measurements of total cross sections for jj production in association with two Z bosons are also reported. Limits on anomalous quartic gauge couplings are derived in terms of


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 E-mail address: cms -publication -committee -chair @cern .ch. 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 QCDinduced 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 highmass 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 pp 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].

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. The first level of the CMS trigger system, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events of interest with a latency of 3.2 μs. The high-level trigger processor farm further decreases the event rate from around 100 kHz to less than 1 kHz, before data storage [19]. A more detailed description of the CMS detector, together with a def-inition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [20].

Signal and background simulation
Several Monte Carlo (MC) event generators are used to simulate signal and background contributions. The simulated samples are employed to optimize the event selection, evaluate the 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 contribution of electrons and muons from τ decays to the signal is very small and is therefore neglected.
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 ). The resulting corrections range from 9%, at values of m ZZ close to 180 GeV, to 5%, for high m ZZ values. Additional NLO EW corrections are applied for m ZZ > 2m Z , following the calculations from Ref. [28].
These corrections become larger with increasing values of m ZZ and are below 5% for m ZZ < 600 GeV.
The interference between the EW and QCD diagrams is evaluated using dedicated samples produced with MG5 at LO, via the direct generation of the interference term between the two processes.
The loop-induced production of two Z bosons from a gluongluon (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 [29]. For the 1-and 2-jet contributions, a pp initial state instead of gg is specified in MG5 to also include initialstate 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 [30,31]. An NLO/LO K factor, which is extracted from Refs. [32,33], 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 [34] package versions are used for parton showering, hadronization and the underlying event simulation, with parameters set by the CUETP8M1 tune [35] (CP5 tune [36]) for the 2016 (2017 and 2018) data-taking period. The NNPDF3.0 (NNPDF3.1) set of parton distribution functions, PDFs [37], is used for the 2016 (2017 and 2018) data-taking period. Unless specified otherwise, the simulated samples are normalized to the cross sections obtained from the respective event generator.
The detector response is simulated using a detailed description of the CMS detector implemented in the Geant4 package [38,39]. 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 [40].
The primary triggers require the presence of a pair of loosely isolated leptons, whose exact requirements depend on the datataking year. Triggers requiring three leptons with low transverse momentum (p T ), as well as isolated single-electron and singlemuon 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 [41] 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 [42,43] 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 [44].
Muons are reconstructed by combining information from the silicon tracker and the muon system [45]. The matching between the muon-system and tracker tracks proceeds either outside-in, 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 threedimensional 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 scalar 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. [46]. 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 [47] 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 [42], as implemented in the Fast-Jet package [43], 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 [48].
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 [49,50] 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 [51].
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 signal-enriched region that requires m jj > 400 GeV and | η jj | > 2.4 and corresponds to a signal purity of ≈20%, and a tight VBS signalenriched 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.

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 heavyflavor 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 sameflavor 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 con-tribution from asymmetric photon conversions, and p miss T < 25 GeV to suppress the WZ contribution.
We validate the procedure using a second control region from opposite-sign same-flavor leptons that fail the selection criteria. The procedure is identical to that used in Ref. [40].

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 from each source of uncertainty is summarized below. All quoted ranges correspond to variations for the different leptonic final states and fiducial analysis regions.
Renormalization and factorization scale uncertainties are evaluated by varying both scales independently. The following vari- , taking the largest variation as the systematic uncertainty, which is about 6% for the EW signal, 11% for the interference term, and ranges from 10 to 12% for the qq → ZZjj QCD background, which is described at a higher QCD order.
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. [32,33]. The PDF and related α S variations are evaluated from the variations of the respective eigenvalues set following the NNPDF prescription [37], 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.
Although in all simulated samples additional partons are described at the LO in matrix-elements or better, we investigate residual uncertainties from parton-shower modeling. Following the prescription from Ref. [52], the renormalization scales are varied independently for the initial-and final-state radiations by factors of 0.5 and 2, and alternative samples are simulated using herwig 7 [53] with the CH3 tune [54]. The largest deviation from the nominal value is used as the uncertainty. On average it ranges from 4%, for the gg → ZZjj background and EW signal, to 5% for qq → ZZjj, and is up to 16% at the lowest values of K D .
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 [50] is 2.2-6.3% (0.2-0.4%). The uncertainty in the trigger as well as the lepton reconstruction 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 [55][56][57]. 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%, while theory uncertainties range 9-12%.

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]: it utilizes matrix element calculations for the EW ZZjj and qq → ZZjj processes from mcfm [58] and employs both the kinematical distributions of leptons and jets to separate signal from background.
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. 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 and signal strength 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-dominated region of the discriminant distribution. The signal strength of the EW signal in the ZZjj inclusive selection is also determined from the same fit. Separate fits are used to determine the EW signal strengths in the other two analysis regions. Fits that only use the event counts in the three regions are performed to determine the signal strengths of the EW+QCD ZZjj production.

Table 2
Particle-level selections used to define the fiducial regions for EW and EW+QCD cross sections. The systematic uncertainties in shape and normalization are treated as nuisance parameters in the fits and profiled [59]. The size of the interference between the EW and QCD production is very small (9% and 3.5% of the EW signal in the ZZjj inclusive and VBS-enriched tight region, respectively). Its effect is included in the EW signal fits via a square-root scaling of the signal strength, approximated with a linear expansion to simplify the fitting technique, while it is neglected in the EW+QCD fits.

ZZjj inclusive Leptons
The measured signal strengths from the fits 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 are 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 reconstructionlevel 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, including the relative K factors where applicable. For the EW ZZjj prediction, in addition, we compare to higher-order calculations at NLO in QCD [60,61] and with a  theoretical prediction at LO in QCD, but including NLO EW corrections [62]. Uncertainties in all SM predictions come from variations of the factorization and renormalization scales. PDF+α S variation uncertainties are summed in quadrature, except from the prediction from Ref. [62] for which only the uncertainty in the scale variation is available. The

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 the Wilson coefficients of the aforementioned operators, under the hypothesis of absence of anomalies in triple gauge couplings. Fig. 4 shows the expected m 4 distributions in the ZZjj inclusive region, 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 ex- Table 3 Measured cross sections and corresponding SM predictions in the three fiducial regions. The reported SM predictions include those extracted from generated events in MC samples adopted for the analysis (LO), as well as higher-order calculations at NLO in QCD [60]  hibits 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 analysis employs the same methodology used for the signal strength, including the profiling of the systematic uncertainties. Using two different approaches, the distributions of the SM processes, including the EW component, are either normalized to their measured values in the EW signal extraction (as discussed in Section 7) or to their expected values. The Wald Gaussian approximation and Wilks' theorem are used to derive 2σ confidence level (CL) intervals on the aQGC parameters [63][64][65]. The measurement is statistically limited. 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. [66] 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.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
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.  Expected and observed limits of the 2σ CL intervals on the couplings of the quartic operators T0, T1, and T2, as well as the neutral current operators T8 and T9. Observed limits in parentheses are obtained by using the prefit normalization of SM processes. The unitarity bounds are also listed. All coupling parameter limits are in TeV −4 , while the unitarity bounds are in TeV.

Coupling
Exp