Measurement of the electroweak production of Z$\gamma$ and two jets in proton-proton collisions at $\sqrt{s} =$ 13 TeV and constraints on anomalous quartic gauge couplings

The first observation of the electroweak (EW) production of a Z boson, a photon, and two forward jets (Z$\gamma$jj) in proton-proton collisions at a center-of-mass energy of 13 TeV is presented. A data set corresponding to an integrated luminosity of 137 fb$^{-1}$, collected by the CMS experiment at the LHC in 2016-2018 is used. The measured fiducial cross section for EW Z$\gamma$jj is $\sigma_{\mathrm{EW}}$ = 5.21 $\pm$ 0.52 (stat) $\pm$ 0.56 (syst) fb = 5.21 $\pm$ 0.76 fb. Single-differential cross sections in photon, leading lepton, and leading jet transverse momenta, and double-differential cross sections in $m_{\mathrm{jj}}$ and $\lvert\Delta\eta_{\mathrm{jj}}\rvert$ are also measured. Exclusion limits on anomalous quartic gauge couplings are derived at 95% confidence level in terms of the effective field theory operators $\mathrm{M}_{0}$ to $\mathrm{M}_{5}$, $\mathrm{M}_{7}$, $\mathrm{T}_{0}$ to $\mathrm{T}_{2}$, and $\mathrm{T}_{5}$ to $\mathrm{T}_{9}$.


Introduction
Vector boson scattering (VBS) processes are purely electroweak (EW) interactions at leading order (LO). In a proton-proton (pp) collision where two vector bosons radiated from the incoming quarks scatter, the two outgoing quarks appear as jets widely separated in pseudorapidity (η) and with a large dijet mass (m jj ), providing a unique signature. The VBS is of great interest because of the role of the Higgs boson in restoring unitarity to the VBS cross section. Studies of VBS complement direct Higgs boson measurements [1][2][3][4][5] and open a window to beyond the standard model (BSM) scenarios at energy scales outside the reach of direct searches. Standard model (SM) EW production of Zγjj, which can proceed through VBS, can be extracted by exploiting the unique features of the VBS signature. A precise measurement of EW Zγjj production is sensitive not only to quartic gauge couplings (QGCs) in the SM as well as possible anomalous QGCs (aQGCs) [6], but also to triple gauge couplings (TGCs) and anomalous TGCs. Since the latter is well constrained in diboson production [7], they are not explored in this paper.
In this paper, we present the first observation of EW Zγjj production. Events corresponding to the + − γjj final states where = e or µ, are selected, and the dijet system is required to satisfy the typical VBS signature. Figure 1 shows the representative Feynman diagrams (upper left and center), in which a Z boson and a photon are produced in a scattering interaction between two W bosons, an example of the production of Zγjj induced by quantum chromodynamics (upper right), one of the important backgrounds in this study, and examples of non-VBS EW production of Zγjj (lower). Previous experimental results for EW Zγjj production have been reported by the ATLAS [8] and CMS [9] Collaborations based on data collected in 2016 at √ s = 13 TeV, corresponding to integrated luminosities of 35.9 fb −1 and 36 fb −1 , respectively. The observed (expected) significance reported by ATLAS was 4.1 (4.1) standard deviations (SD), and by CMS was 4.7 (5.5) SD.
The measurements reported here are based on pp collision data at √ s = 13 TeV, collected from 2016 to 2018 by the CMS experiment at the CERN LHC, corresponding to an integrated luminosity of 137 fb −1 . A simultaneous maximum likelihood (ML) fit is used to extract the signal significance, the signal strength, fiducial cross sections, and unfolded differential cross sections for both EW and EW+QCD production of Zγjj. Unfolded differential cross sections are measured as functions of three 1-dimensional observables (the transverse momenta (p T ) of the leading lepton, photon, and jet) and one 2-dimensional observable (m jj and |∆η jj |). Using the selected + − γjj events with high photon p T , constraints on BSM contributions to the VVZγ vertices, where V = W, Z, or γ, are determined in an effective field theory (EFT) framework using dimension-8 operators [6]. This paper is arranged as follows. The CMS detector and the event simulation are summarized in Secs 2 and 3. The object reconstruction and the event selection are described in Sec 4. The estimation of the main backgrounds is given in Sec 5. The systematic uncertainties are discussed in Sec 6 and the results are presented in Sec 7. The paper is summarized in Sec 8.

Object reconstruction and event selection 4.1 Object reconstruction
The final state of interest consists of a pair of oppositely charged isolated leptons, a photon, and two jets. We employ both a global object reconstruction algorithm and dedicated particletype (e.g., muon) reconstruction algorithms, for different purposes, as described in the following paragraphs. The global object reconstruction algorithm, called the particle-flow (PF) algorithm [26], reconstructs and identifies all particle candidates (photons, electrons, muons, charged and neutral hadrons, and missing transverse momentum) in an event, based on a combination of information from the various elements of the CMS detector; the result is a set of physics objects called PF candidates.
The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [27,28] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets, which includes the leptons.
Electron candidates are reconstructed by combining information from the ECAL and the tracker within |η| < 2.5 and p T > 25 GeV. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. Reconstructed electrons are re-quired to satisfy additional identification requirements [29] as follows: the relative amount of energy deposited in the HCAL; a matching procedure between the trajectory in the inner tracker and that in the supercluster [30] of the ECAL; the number of missing hits in the inner tracker; the consistency between the track and primary vertex; a shower shape variable σ iηiη , which quantifies the transverse spread in η of the electromagnetic shower in the ECAL (discussed in Sec 5); and a photon conversion rejection algorithm. An appropriate working point, referred to as the stringent electron selection, which has an average per-electron efficiency of 80%, is used to identify the electron candidates from the Z boson decays in the signal process. A far less restrictive working point, referred to as the minimal electron selection, which has an average per-electron efficiency of 95%, is used to remove events that contain additional electrons.
Muon candidates are reconstructed by combining information from the silicon tracker and the muon system within the region |η| < 2.4 and p T > 20 GeV [31]. The combined information is used to produce a global track fit, and the muon momenta are obtained from the track curvatures. Muons are selected from the reconstructed muon track candidates by applying additional identification requirements as follows: the number of hits in the muon system and the inner tracker; the quality of the combined fit to a track; the number of matched muon-detector planes; and the consistency between the track and primary vertex. Different muon identification working points are defined according to their efficiency. An appropriate working point, referred to as the stringent muon selection, is used to identify the muon candidates from the Z boson decays in the signal process. The efficiency to reconstruct and identify muons is greater than 96%. A far less restrictive working point, referred to as the minimal muon selection, is used to remove events with additional muons.
Leptons are required to be isolated from other particles in the event. The relative isolation is used and defined as: where the sums run over the charged and neutral hadrons, as well as the photons, in a cone of size ∆R = √ (∆η) 2 + (∆φ) 2 = 0.3 (0.4) for the electron (muon) trajectory, where η and φ denote the pseudorapidity and azimuthal angle. The quantity ∑ p charged T is the scalar p T sum of charged hadrons originating from the primary vertex; ∑ p neutral T and ∑ p γ T are the respective scalar p T sums of neutral hadrons and photons. The contribution from PU in the isolation cone, p PU T , is subtracted using the FASTJET v3.0.2 technique [28]. For electrons, p PU T is evaluated using the "jet area" method described in Ref. [32]. For muons, p PU T is assumed to be half of the scalar p T sum deposited in the isolation cone by charged particles not associated with the primary vertex. The factor of one half corresponds approximately to the ratio of neutral to charged hadrons produced in the hadronization of PU interactions. Electrons are considered isolated for the stringent (minimal) working points if R Iso < 0.0695 (0.1750) in the barrel and R Iso < 0.0821 (0.1590) in the endcap detector regions. Muons are considered isolated for the stringent (minimal) working points if R Iso < 0.15 (0.25).
The efficiencies of lepton reconstruction and selection are measured as a function of p T and η for both data events and MC events. The "tag-and-probe" technique [33] is used on events containing a single Z boson. The ratio of efficiencies from events in data and MC are used to correct the simulation. The momentum scales of both muons and electrons are calibrated in bins of p T and η [29,34].
Photon reconstruction and selection are similar to electron reconstruction and selection. Photons located in the barrel region, |η| < 1.442, and the ECAL endcap region, 1.566 < |η| < 2.500, with p T > 20 GeV, are referred to as γ barrel and γ endcap , respectively. The region 1.442 < |η| < 1.566 is a transition region between the barrel and endcaps and is not used for photon reconstruction, because the reconstruction of a photon object in this region is less precise. Reconstructed photons are required to satisfy further quality criteria [29] based on the following quantities: the relative amount of deposited energy in the ECAL and HCAL; isolation variables constructed separately for the charged and neutral hadrons, photons other than the signal photon; and a procedure that quantifies the likelihood for a photon to originate from electron bremsstrahlung. An appropriate working point, referred to as the stringent photon selection, with an average per-photon efficiency of 80%, is used to reconstruct prompt photons (not from hadron decays) in the final state. A second working point, which will be referred to as the nonprompt-enriched photon selection, is used to reconstruct nonprompt photons that are mainly products of neutral pion and η meson decays [29,35], which constitute an important background in this study. Jets are reconstructed from particle-flow candidates using the anti-k T jet clustering algorithm [27], as implemented in the FASTJET package v3.0.2, with a distance parameter of 0.4. The energies of charged hadrons are determined from a combination of their momenta measured in the tracker and the matching of ECAL and HCAL energy deposits, corrected for the response of the calorimeters to hadronic showers. The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies. To reduce the instrumental background, as well as the contamination from PU, jets are selected by a stringent jet selection based on the multiplicities and energy fractions carried by charged and neutral hadrons.
Jet energy corrections are extracted from data and simulated events to account for the effects of PU, non-uniformity of the detector response, and residual differences between the jet energy scale (JES) in data and simulation. The JES calibration [36] relies on corrections parametrized 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 PU contamination, all jets must satisfy |η| < 4.7 and have a corrected p T > 30 GeV. Jets from PU are further rejected using PU jet identification criteria based on a multivariate technique [37].

Event selection
The events of interest are selected by dilepton triggers. We denote the lepton having the larger p T as 1 and the other one as 2. The p T thresholds in the HLT are p 1 T > 23 GeV, p 2 T > 12 GeV for electrons, and p 1 T > 17 GeV, and p 2 T > 8 GeV for muons. In 2016 and 2017, a timing shift in the ECAL endcap was not properly accounted for in the trigger logic, resulting in the trigger decision sometimes mistakenly being assigned to the previous bunch crossing. This led to a sizable decrease in the L1 trigger efficiency for events with high energy deposits in the ECAL endcaps. The loss of efficiency for EW Zγjj events associated with this effect is ≈8% for the invariant mass of two jets m jj > 500 GeV, and increases to ≈15% for m jj > 2 TeV. A correction is therefore applied as a function of jet p T and η using an unbiased data sample with correct timing and is implemented through a factor that represents the probability of the event to avoid having mistimed signals.
Selected events are required to contain two oppositely charged same-flavor leptons, either a pair of electrons or a pair of muons. Both leptons must pass the stringent working points and the corresponding isolation requirements described in Sec 4.1 and must satisfy p T > 25 (20) GeV, and |η| < 2.5 (2.4) in the electron (muon) case. The invariant mass of the dilepton system (m ) is required to be within the window 70 < m < 110 GeV. To reduce the WZ and ZZ backgrounds, events are rejected if there is any additional lepton passing the less restric-tive identification criteria mentioned in Sec 4.1. At least one photon satisfying the stringent identification criteria is required. The photon with the largest p T in the event is used if there is more than one photon passing the stringent identification criteria. The photon is required to have p T > 20 GeV. The ∆R between the selected photon and each of the selected leptons is required to be larger than 0.7. The invariant mass of the dilepton-photon system must satisfy m Zγ > 100 GeV to reduce the contribution from final-state radiation in Z boson decays. The event must also contain at least two jets that satisfy the jet identification criteria described in Sec 4.1 and that are separated from selected leptons and photons by ∆R > 0.5. The jet with the largest p T is referred to as the leading jet and is denoted as j1, and the jet with the lower p T is denoted as j2. The jets are required to satisfy p T > 30 GeV, |η| < 4.7, and ∆R(j1, j2) > 0.5. The dilepton selection with a photon and two jets is henceforth referred to as the "common" selection.
The signal region is defined by the common selection, and by requiring m jj > 500 GeV and |∆η jj | = |η j1 − η j2 | > 2.5. Two additional criteria are used for the signal significance measurement. First, the Zeppenfeld variable [38] η * = |η Zγ − (η j1 + η j2 )/2|, where η Zγ represents the η of the Zγ system, needs to satisfy η * < 2.4. Second, the magnitude of the difference between the azimuthal angle of the Zγ system (φ Zγ ) and the azimuthal angle of the dijet system (φ j1j2 ), ∆φ Zγ,jj = |φ Zγ − φ j1j2 |, which should be large in signal events because the two systems are recoiling against each other, must satisfy ∆φ Zγ,jj > 1.9. A low-m jj control region, in which the EW signal is negligible compared with the contribution from QCD-induced Zγjj production, is defined by the common selection and the requirement 150 GeV < m jj < 500 GeV.
The total and differential cross sections for EW Zγjj and EW+QCD Zγjj production are measured in a fiducial region (see "Fiducial volume" in Table 1) that closely mirrors the EW signal in the VBS at the particle level, in which the requirements of η * and ∆φ Zγ,jj for the optimization of signal significance used in the EW signal region (see "EW signal region" in Table 1) are not applied. The particle-level leptons and photons are required to be prompt, which means that the photon should be from the VBS process and lepton should be from the Z decay, and the momenta of prompt photons with ∆R γ < 0.1 are added to the lepton momenta to correct for final-state photon radiation.
The selection for the aQGC search is similar to the EW signal selection, but targets the characteristic high-energy behavior of aQGC processes by requiring p γ T > 120 GeV. In this high photon p T region, the background has been suppressed as much as possible, so the requirements of η * and ∆φ Zγ,jj are not applied.
A summary of all the selection criteria is displayed in Table 1.

Background estimation
The dominant background arises from the QCD-induced production of Zγjj. The yield and shape of this irreducible background are taken from simulation, but are ultimately constrained by the data in the ML fit mentioned in Sec 7.1 that extracts the EW signal. The second most important background arises from events in which the selected photon is not prompt and is mainly from Z+jets events. This background cannot be simulated accurately and is estimated from data, as described in the following paragraph. Other small contributions feature kinematic distributions similar to that of the dominant background and are estimated from simulation including single top quark events in the sand t-channels that are normalized to their respective NLO cross sections; associated single top quark and W boson production normalized to its next-to-next-to-leading order (NNLO) cross section [39]; WW production normalized to Table 1: Summary of the five sets of event selection criteria used to define events in the fiducial cross section measurement region, control region, EW signal extraction region, and the region used to search for aQGC contributions.
The background from events containing a nonprompt photon is estimated using data to calculate the event weight of the corresponding nonprompt photon event shown in Eq. (2), in a region similar to our common selection with the jet requirements removed. The numerator n data represents the number of events passing the stringent photon selection. The denominator N unweighted fake represents the number of events passing the nonprompt-enriched photon selection mentioned in Sec 4.1. The contribution from the signal region is removed for both the numerator and denominator, and the prompt contribution in the denominator is also removed by subtracting the small prompt photon contribution from the data based on simulated samples. The factor fake-fraction is the fraction of nonprompt photons in the region where the stringent photon selection is applied, which is obtained from a template fit of the photon σ iηiη distribution, since the variable σ iηiη quantifies the width of the photon electromagnetic shower in η, which is narrow for prompt and broad for nonprompt photons. The prompt-photon template is obtained from simulated Zγ events and the nonprompt-photon template is obtained from a sideband method of inverting the charged hadron isolation variable of the photon in data. The event weight of an event containing a nonprompt photon can be then calculated as a function of p γ T for photons in the barrel and endcap regions. The nonprompt-photon background estimate in the signal region is determined by these event weights and the rate of events passing the signal region selection with the stringent photon selection replaced by the nonprompt-enriched photon selection.

Systematic uncertainties
The sources of systematic uncertainty can be divided into experimental and theoretical categories. The experimental sources include uncertainties in corrections to the simulation, the method to estimate the nonprompt-photon background contribution, and corrections for detector effects during data taking not properly accounted by the simulation in the experiment. The sources of theoretical uncertainty include the choice of the renormalization and factorization scales, and the choice of the PDFs.
The uncertainty because of renormalization and factorization scales denoted as µ R and µ F , respectively, is evaluated for the signal and QCD-induced Zγjj background. The different choices for µ F and µ R considered are these six combinations: (µ F , µ R ) = (2µ 0 , µ 0 ), (0.5µ 0 , µ 0 ), (µ 0 , 2µ 0 ), (µ 0 , 0.5µ 0 ), (2µ 0 , 2µ 0 ), and (0.5µ 0 , 0.5µ 0 ), in which µ 0 represents the nominal scale. These six combinations are further divided into three groups according to either µ R or µ F not equalling µ 0 or both of them not equalling µ 0 . The difference in one group is calculated and regarded as one component of the uncertainties in µ R and µ F . The uncertainties of these three components range from 1 to 12% for the EW signal and from 6 to 25% for QCD-induced Zγjj background. All three components are included as independent systematic uncertainties in the ML fit introduced in Sec 7.1. The PDF and related strong coupling α S are evaluated using the eigenvalues of the PDF set following the NNPDF prescription [40]. The size of this uncertainty is 1%-3% for both EW signal and QCD-induced Zγjj background.
The uncertainties in the jet modeling, which include the JES, the jet energy resolution (JER), and the PU jet identification, are calculated in simulated events. The JES and JER are obtained by applying the corrections shifted by ±1 SD. The effect of the updated corrections is propagated to all dependent variables and all selection criteria, and these effects on the yield are determined bin-by-bin for each bin of the m jj -|∆η jj | distribution. The variation, +1 SD or −1 SD, that has the larger absolute effect on the yield is assumed as the uncertainty. The size of the uncertainty varies between 1 and 92%, depending on the m jj -|∆η jj | bin, but the larger values typically correspond to bins where the uncertainty is less important because they are applied to a smaller nominal yield. The uncertainties in PU jet identification are calculated by changing the corresponding scale factors ±1 SD, depending on p T and η. The uncertainties in jet energy corrections and PU jet working points are uncorrelated between different years, but correlated between the electron and muon categories.
The systematic uncertainty in the nonprompt-photon background estimate is the quadratic sum of three components. The first component is the uncertainty in the choice of the isolation sideband. The second component is the uncertainty in the bias in the fitting procedure, calculated by performing the procedure in simulated pseudodata and comparing the fit results with the known fractions. This component, which is larger in the endcap than in the barrel, increases with photon p T . The third component is the uncertainty in the modeling of the prompt-photon events in the true template fit, estimated as the difference between promptphoton events simulated from QCD-induced Zγ and EW Zγjj signal. The overall uncertainty in the nonprompt-photon background estimation ranges from 9 to 37%, uncorrelated between years and channels.
The uncertainties that arise from the finite number of events in simulated samples and data control regions, referred to as statistical uncertainties, are calculated bin-by-bin based on the Poisson distribution. The statistical uncertainties typically increase with increasing m jj and |∆η jj | and are uncorrelated across different processes and bins. The uncertainties in the efficiencies of lepton identification, trigger efficiencies, and photon identification range from 0.5 to 3.0% depending on p T and η and include both statistical and systematic sources. The systematic (statistical) uncertainties are dominant in the low (high) p T range.
All simulated samples are also affected by uncertainties associated with the ECAL timing shift, the reweighting of the PU distribution, and the integrated luminosity. The uncertainties associated with the ECAL timing shift correction factors range from 1 to 4%, depending on the process and the m jj -|∆η jj | bin. The uncertainty from pileup reweighting is evaluated by changing the total inelastic cross section of 69.2 mb [41] by ±5%, which results in an uncertainty in the 1%-10% range. The integrated luminosity of the 2016, 2017, and 2018 data-taking periods have uncertainties in the 1.2-2.5% range [42][43][44]. These uncertainties are partially correlated and correspond to a total uncertainty of 1.6%. The uncertainties of ECAL mistiming, reweighting of the PU, and integrated luminosity uncertainties are correlated across years.
All of the above systematic uncertainties are applied in the calculation of the signal significance, measurements of the cross section, and in the search for aQGCs. They are also applied in the cross section measurement, with the exception of the theoretical uncertainties related to the normalization of the signal.
We organize the systematic uncertainties into the following groups: theoretical uncertainties including the scales µ R and µ F , and the PDF uncertainties; corrections applied to jets including JES and JER; uncertainties in the nonprompt-photon background estimate; statistical uncertainties from simulation or data; corrections applied to electron and photon candidates; corrections applied to muon candidates; uncertainties from PU reweighting; uncertainties in the integrated luminosity; uncertainties from the L1 trigger timing shift; and uncertainties from corrections of PU jet identification (ID) working points. The remaining uncertainties, which are referred to as "other", include the uncertainties in the cross section estimation of diboson and tt γ processes that have an impact <0.1%. The uncertainties in µ R and µ F of the QCD-induced production and in the jet energy correction are the dominant systematic uncertainties in the measurement. The impact of the uncertainty of each group on the signal strength measurement is displayed in Table 2.

Results
The pre-fit (before the simultaneous fit described in Sec 7.1) signal and background expected yields, as well as the observed data yields in the signal region, are shown in Figs. 2 and 3 separately for the photon in the ECAL barrel and endcaps, for both the dielectron and dimuon channels.

Measurement of the signal significance
To extract as much information as possible from the data set, we perform a likelihood-based statistical analysis. An optimal binning leads to the best expected signal significance used. The likelihood function is the product of the binned signal and backgrounds probability density functions (pdf), which are expected to follow the Poisson distribution: where n i is the number of observed events in data, s i and b i are the expected event yields for the signal and backgrounds, N represents the number of bins in the signal and control regions, and p(˜ θ| θ) is a Gaussian constrained pdf of the effects of the systematic uncertainties θ called nuisance parameters (NP), in which˜ θ represents the external measurements corresponding to each NP. The parameter of interest (POI) µ is the signal strength, which represents the ratio of observed to expected signal yields. The POI µ is estimated by maximizing the profile likelihood ratio defined as: The numerator of this ratio is the profile likelihood function. The quantityˆ θ denotes the value of θ that maximizes L for the specified µ; it is the conditional ML estimator of θ (and thus is a function of µ). The denominator is the maximized (unconditional) likelihood function,μ .
The p-value is then converted to a significance based on the area in the tail of a normal distribution.
The control and signal regions are each divided separately for photons in the ECAL barrel/endcaps and for the dielectron and dimuon channels. The signal region is divided further into bins in m jj and |∆η jj |, as shown in Figs. 4 and 5. The control region is divided into bins in m jj , as shown in Figs. 6 and 7. The postfit yields for every process and for data are listed in Table 3.
The main contributions to the significance are from the bins in the signal region with an excess of signal relative to background events, such as high m jj bins in each channel.
The observed (expected) significance for EW signal obtained from a simultaneous fit of all bins in the signal region and the control region to the data (Asimov data set [46]) is 9.4 (8.5 SD).  Events /bin

Fiducial cross section
Fiducial cross sections are measured in the fiducial region, which is designed to mirror the signal region as closely as possible, as shown in Table 1. The fiducial cross section is extracted using the same binning of m jj and |∆η jj | as used in the signal significance measurement, and with the same simultaneous fit in all regions and channels, with the following exception: the events that pass the EW signal selection but fail the fiducial selection are regarded as a background. We define the fiducial cross section as: where σ g is the cross section for the generated signal events,μ is the signal strength parameter, and a gf is the acceptance for the events generated in the fiducial region and evaluated through simulation.   A combined EW+QCD Zγjj cross section is also measured in the same fiducial region using the same procedure, except that the control region is excluded. In this measurement, both the EW and QCD contributions are considered signal. The combined Zγjj cross section is defined as  The theoretical fiducial cross section for QCD Zγjj production is 8.93 ± 1.70 (scale) ± 0.08 (PDF) fb. The expected fiducial cross section for the combined QCD and EW Zγjj production is 13.3 ± 1.72 (scale) ± 0.10 (PDF) fb. The best fit value for the combined Zγjj signal strength and the measured cross section are

Unfolded differential cross section distribution
Unfolding is used to correct measured detector-level distributions to the particle-level. It accounts for the limited acceptance and efficiencies of the detector, as well as for bin-to-bin migration between the measured and corrected distribution arising from detector resolution. Simulated EW samples from MC event generators are used to perform the unfolding. Distributions obtained from the generated events correspond to the particle level (which will be referred to as "gen"). The same distributions obtained using simulated events correspond to detector level (which will be referred to as "reco"). In the well-defined fiducial phase space, the events in "reco" and "gen" follow the relation: can also be derived from the principle of ML; the corresponding likelihood is in Eq. (11).
This solution is based on the assumption that R is nonsingular and insensitive to small perturbations. If that is not the case, the problem is ill-posed and R is ill-conditioned, which means that a regularization [47] method is needed. The sensitivity to fluctuations associated with the ML solution can be quantified by the condition number of R, which is a measure of how illconditioned the problem is. In general, if the condition number c(R) is small (less than 10), then the problem can most likely be solved using the unregularized ML estimate. In this paper, an appropriate choice of binning is made to guarantee that the response matrix R is nonsingular and c(R) < 10. The regularization parts thus are not considered in this unfolded differential cross section measurement.
The unfolded differential cross section is measured in the same way as the fiducial cross section, with a simultaneous fit to data in the control region and signal region in bins of the single variables leading photon p T , leading lepton p T , and leading jet p T , and in bins of two variables m jj and |∆η jj |. The off-diagonal elements of the response matrix correspond to 1% for leading photon p T , 1% for leading lepton p T , 3% for leading jet p T , and around 5% for m jj and |∆η jj |. All systematic uncertainties discussed in Sec 6 are also considered for the corresponding process, especially for the evaluation of R ij . The signal strengths for each bin of the EW Zγjj unfolded distribution are listed in Tables 4 and 5 and the corresponding differential cross section distributions are shown in Fig. 8. A combined EW+QCD Zγjj unfolded differential cross section is also measured in the same region using the same procedure, except that the control region is excluded. In this measurement, both the EW and QCD contributions are considered signal. The combined Zγjj unfolded  0.00060 ± 0.00016 0.00112 ± 0.00092 differential cross section is shown in Fig. 9, and Tables 6 and 7. Within the uncertainties, the unfolded distributions agree well with the SM predictions.  Figure 9: Unfolded differential cross section as a function of the leading lepton p T , leading photon p T , leading jet p T , and m jj -|∆η jj | for EW+QCD Zγjj. The black points with error bars represent the data and their statistical uncertainties, whereas the red bands represent the total theoretical uncertainties from the MG5 simulation. The last bin includes overflow events.

Limits on anomalous quartic gauge couplings
In an EFT approach to BSM physics, dimension-8 operators are constructed from contractions of the covariant derivative of the Higgs doublet and the charged and neutral field strength tensors associated with gauge bosons. Assuming the SU(2) × U(1) symmetry of the EW gauge field, nine independent charge-conjugate and parity-conserving dimension-8 operators are constructed [6]. The operators affecting the Zγjj channel can be divided into those contain-  ing an SU(2) field strength, the U(1) field strength, and the covariant derivative of the Higgs doublet, L M 0−5,7 , and those containing only the two field strengths, L T 0−2,5−9 . The coefficient of the operator L X Y is denoted by F XY /Λ 4 , where Λ is the unknown scale of BSM physics.
The effects of the aQGCs in addition to the SM EW Zγ process, as well as interference between the EW and QCD-induced processes, are simulated as described in Sec 3. The contribution from aQGCs enhances the production of events at large Zγ mass, so the m Zγ distribution is used to extract limits on the aQGC parameters. To obtain the prediction for the signal as a function of the value of the aQGC parameter, a quadratic function is used to fit the ratio of aQGC and SM yields in each bin of m Zγ in the aQGC region defined in Sec 4.2. From Fig. 10, no statistically significant excess of events relative to the SM prediction is observed.
The likelihood function is the product of Poisson distributions and a normal constraining term with nuisance parameters representing the sources of systematic uncertainties in any given bin. The final likelihood function is the product of the likelihood functions of the electron and muon channels. This test statistic, is essentially the same test statistic as in Sec 7.1 except that the α test represents the aQGC point being tested. The symbol θ represents a vector of nuisance parameters assumed to follow lognormal distributions. The parameterˆ θ corresponds to the maximum of the likelihood function at the point α test . Theα andˆ θ parameters correspond to the global maximum of the likelihood function.
This test statistic is assumed to follow a χ 2 distribution [48]. It is therefore possible to extract the limits immediately from twice the difference in the log-likelihood function 2∆NLL = t α test [49]. The observed and expected 95% confidence level limits for the coefficients, shown in Table 8, are obtained by varying the coefficients of one operator at a time and setting all other anomalous couplings to zero. The statistical uncertainty is dominant in the limits setting. The unitarity bound is defined as the scattering energy at which the aQGC coupling strength is set equal to the observed limit that would result in a scattering amplitude that violates unitarity [50]. These results provide the most stringent limit to date on the aQGC parameter F T9 /Λ 4 .

Summary
This paper presents the first observation of the electroweak (EW) production of a Z boson, a photon, and two jets (Zγjj) in √ s = 13 TeV proton-proton collisions recorded with the CMS detector in 2016-2018 corresponding to an integrated luminosity of 137 fb −1 . Events were selected by requiring two opposite-sign leptons with the same flavor from the decay of a Z boson, one identified photon, and two jets that have a large separation in pseudorapidity and a large dijet mass. The measured cross section in the fiducial volume defined in Table 1 for EW Zγjj production is 5.21 ± 0.52 (stat) ± 0.56 (syst) fb = 5.21 ± 0.76 fb, and the fiducial cross section of EW and QCD-induced production is 14.7 ± 0.80 (stat) ± 1.26 (syst) fb = 14.7 ± 1.53 fb. Both the observed and expected signal significances are well in excess of 5 standard deviations. Differential cross sections for EW and EW+QCD are measured for several observables and compared to standard model predictions computed at leading order. Within the uncertainties, the measurements agree with the predictions. Constraints are set on the effective field theory dimension-8 operators M 0 to M 5 , M 7 , T 0 to T 2 , and T 5 to T 9 , giving rise to anomalous quartic gauge couplings. These constraints are either competitive with or more stringent than those previously obtained.

Acknowledgments
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 centers and personnel of the Worldwide LHC Computing Grid and other centers 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, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: BMBWF and FWF (Austria); FNRS and FWO (Belgium); CNPq,  [4] ATLAS, CMS Collaboration, "Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at √ s = 7 and 8 TeV", JHEP 08 (2016) 045, doi:10.1007/JHEP08(2016)045, arXiv:1606.02266.   [8] ATLAS Collaboration, "Evidence for electroweak production of two jets in association with a Zγ pair in pp collisions at √ s = 13 TeV with the ATLAS detector", Phys. Lett. B 803 (2020) 135341, doi:10.1016/j.physletb.2020.135341, arXiv:1910.09503.