Measurements of $Z\gamma$ and $Z\gamma\gamma$ production in $pp$ collisions at $\sqrt{s}=$ 8 TeV with the ATLAS detector

The production of $Z$ bosons with one or two isolated high-energy photons is studied using $pp$ collisions at $\sqrt{s}$ = 8 TeV. The analyses use a data sample with an integrated luminosity of 20.3 fb$^{-1}$ collected by the ATLAS detector during the 2012 LHC data taking. The $Z\gamma$ and $Z\gamma\gamma$ production cross sections are measured with leptonic ($e^{+}e^{-}$, $\mu^{+}\mu^{-}$, $\nu\bar{\nu}$) decays of the $Z$ boson, in extended fiducial regions defined in terms of the lepton and photon acceptance. They are then compared to cross-section predictions from the Standard Model, where the sources of the photons are radiation off initial-state quarks and radiative $Z$-boson decay to charged leptons, and from fragmentation of final-state quarks and gluons into photons. The yields of events with photon transverse energy $E_T>$ 250 GeV from $\ell^{+}\ell^{-}\gamma$ events and with $E_T>$ 400 GeV from $\nu\bar{\nu}\gamma$ events are used to search for anomalous triple gauge-boson couplings $ZZ\gamma$ and $Z\gamma\gamma$. The yields of events with diphoton invariant mass $m_{\gamma\gamma}>$ 200 GeV from $\ell^{+}\ell^{-}\gamma\gamma$ events and with $m_{\gamma\gamma}>$ 300 GeV from $\nu\bar{\nu}\gamma\gamma$ events are used to search for anomalous quartic gauge-boson couplings $ZZ\gamma\gamma$ and $Z\gamma\gamma\gamma$. No deviations from Standard Model predictions are observed and limits are placed on parameters used to describe anomalous triple and quartic gauge-boson couplings.


Introduction
This paper is organized as follows. The ATLAS detector is briefly described in Section 2. The signal and background simulation is presented in Section 3. The object and event selections and the background estimation are described in Section 4 and Section 5, respectively. The results of cross-section measurements and their comparison with the Standard Model predictions are presented in Section 6 and Section 7, respectively. The limits on the anomalous triple and quartic gauge-boson couplings are presented in Section 8. Section 9 provides the conclusions.

The ATLAS detector and LHC data sample
The ATLAS detector has been described in detail elsewhere [13]. A short overview is presented here with an emphasis on the electromagnetic calorimeter needed for precision measurement of the high-energy photons. The major components of the ATLAS detector are an inner tracking detector (ID) surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer (MS). The ID is composed of three subsystems. The pixel and silicon microstrip detectors cover the pseudorapidity range |η| < 2.5, while the transition radiation tracker (TRT) has an acceptance range of |η| < 2.0. The TRT provides identification information for electrons by the detection of transition radiation. The MS is composed of three large superconducting air-core toroid magnets, a system of three stations of chambers for tracking measurements with high precision in the range |η| < 2.7, and a muon trigger system effective over the range |η| < 2.4.
The electromagnetic calorimeter is a lead/liquid-argon detector composed of a barrel (|η| < 1.475) and two endcaps (1.375 < |η| < 3.2). For |η| < 2.5 the calorimeter has three layers, longitudinal in shower depth, with the first layer having the highest granularity in the η direction, and the second layer collecting most of the electromagnetic shower energy for high-p T objects. A thin presampler layer covering the range |η| < 1.8 is used to correct for the energy lost by EM particles upstream of the calorimeter. The hadronic calorimeter system, which surrounds the electromagnetic calorimeter, is based on two different detector technologies, with scintillator tiles or liquid-argon as the active medium, and with either steel, copper, or tungsten as the absorber material. Photons are identified as narrow, isolated showers in the EM calorimeter with no penetration into the hadronic calorimeter. The fine segmentation of the ATLAS calorimeter allows efficient rejection of jets fragmenting to high-energy π 0 or η mesons that could be misidentified as isolated direct photons.
Collision events are selected using a three-level trigger system. The first-level trigger is based on custombuilt electronics that use a subset of the total detector information to reduce the data rate to below the design value of 75 kHz. The subsequent two trigger levels run on a processor farm and analyze detector information with greater precision. The resulting recorded event rate from LHC pp collisions at √ s = 8 TeV during the data-taking period in 2012 was approximately 400 Hz. After applying criteria to ensure nominal ATLAS detector operation, the total integrated luminosity useful for data analysis is 20.3 fb −1 . The uncertainty in the integrated luminosity is determined to be 1.9%. It is derived, following the same methodology as that detailed in Ref. [14], from a calibration of the luminosity scale obtained from beam-separation scans.
The + − γ and + − γγ events are selected using single-lepton or dilepton triggers. The p T thresholds are 24 GeV for single-lepton triggers, and 12 GeV (13 GeV) for dielectron (dimuon) triggers. A dimuon trigger with asymmetric muon p T thresholds of 8 GeV and 18 GeV is also used. The ννγ and ννγγ events are selected using a single-photon trigger with a threshold of E T > 120 GeV and a diphoton trigger with a threshold of E T > 20 GeV, respectively. For the events falling within the acceptance of the measurement, the trigger efficiency is close to 100% for e + e − γ(γ) and ννγ final states, about 99% for ννγγ final states, and about 95% for µ + µ − γ(γ) final states.

Simulation of signals and backgrounds
Simulated signal and background events are produced with various Monte Carlo event generators, processed through a full ATLAS detector simulation [15] using Geant4 [16], and then reconstructed with the same procedure as for data. Additional pp interactions (pileup), in the same and neighboring bunch crossings, are overlaid on the hard scattering process in MC simulation. The MC events are then reweighted to reproduce the distribution of the number of interactions per bunch crossing in data. The mean number of interactions per bunch crossing in the dataset considered is 20.7.

Monte Carlo generation of SM pp → Zγ(γ) + X and anomalous gauge-boson couplings processes
The efficiency of the event selection is studied using a MC simulation of the Zγ and Zγγ signals using the Sherpa 1.4 generator [17] with the CT10 parton distribution function (PDF) set [18], and leadingorder (LO) matrix elements with up to three additional final-state partons for Zγ and up to one additional final-state parton for Zγγ. Sherpa uses the CKKW scheme [19,20] to merge matrix elements and parton showers. This "multileg" approach ensures that the first few hardest emissions are modeled by the realemission matrix elements. Sherpa was found to adequately characterize the distributions of selected Zγ candidates in a previous publication [8]. Theoretical uncertainties in the Sherpa predictions in  are taken to be the same as those estimated with MCFM in Section 7.1.
Signal samples with anomalous triple and quartic gauge-boson couplings are generated using Sherpa for aTGC and Vbfnlo 2.7.0 [21][22][23] interfaced to Pythia 8.175 [24] for parton showering, hadronization, and the underlying event for aQGC. More details are given in Section 8.

Monte Carlo generation of background processes
In the measurements of the e + e − γ(γ), µ + µ − γ(γ), and ννγ(γ) production cross sections, backgrounds are estimated either from simulation or from data. The main backgrounds arise from object misidentification and are obtained using data-driven techniques, as described in Section 5. MC simulated backgrounds are used for validation in this case. Smaller backgrounds are estimated directly from simulation.
The WZ and ZZ backgrounds are generated with Powheg-box [25,26] and the CT10 PDF set, with parton showering, hadronization, and the underlying event modeled by Pythia 8.165 with the AU2 set of tuned parameters [27]. The background arising from ttγ is generated with MadGraph5_aMC@NLO 5.2.1.0 [28] and the CTEQ6L1 [29] PDF set, with parton showering, hadronization, and the underlying event modeled by Pythia 8.183. Sherpa 1.4 with the CT10 PDF set is used to simulate τ + τ − γ(γ), γ+jets, and Wγ(γ) events. An alternative MC sample of simulated γ + jet events is generated using Pythia 8.165 with the CTEQ6L1 PDF set. An alternative MC sample of simulated Wγ events is generated using Alpgen 2.14 [30] with the CTEQ6L1 PDF set, interfaced to Herwig 6.520 [31] with Jimmy 4.30 [32] and the AUET2 set of tuned parameters [33] for parton showering, hadronization, and the underlying event.

Selection of Zγ and Zγγ signal events
The event selection criteria are chosen to provide precise cross section measurements of Zγ and Zγγ production, and to provide good sensitivities to anomalous gauge-boson couplings between photons and the Z bosons. The selections are optimized for each of these measurements to obtain high signal efficiency together with good background rejection.

Physics object reconstruction and identification
Collision events are selected by requiring at least one reconstructed primary vertex candidate with at least three charged-particle tracks with p T > 0.4 GeV. The vertex candidate with the highest sum of the p 2 T of the associated tracks is chosen as the event primary vertex. This criterion may choose the wrong primary vertex in ννγ(γ) events. The effect of such a wrong choice was studied in simulation and found to have negligible impact on the photon transverse energy resolution for this analysis.
Electron candidates are reconstructed within the fiducial acceptance region |η| < 2.47 from an energy cluster in the EM calorimeter associated with a reconstructed track in the ID [35]. Photon candidates are reconstructed from energy clusters with |η| < 2.37 [36]. The EM cluster of the electron/photon candidate must lie outside the transition region between the barrel and endcap EM calorimeters, thus electrons and photons with 1.37 < |η| < 1.52 are rejected. The cluster energies are corrected using an in situ calibration based on the known Z boson mass [37]. Clusters without matching tracks are classified as unconverted photon candidates, whereas clusters that are matched to one or two tracks that originate from a conversion vertex are considered as converted photon candidates. Both the unconverted and converted candidates are used in the analysis. Electron tracks are required to be matched to the event primary vertex. The electron d 0 significance, defined as the ratio of the absolute value of the transverse impact parameter, d 0 , with respect to the primary vertex, to its measured uncertainty, must be less than 6.0, and the weighted electron longitudinal impact parameter with respect to the primary vertex |z 0 × sin θ| must be less than 0.5 mm. Reconstructed electrons are required to have p T > 25 GeV. The photon E T threshold depends on the analysis channel.
Muon candidates are identified, within pseudorapidity |η| < 2.5, by matching complete tracks or track segments in the MS to tracks in the ID [38]. Similarly to electrons, the muon candidates are required to be matched to the primary vertex with a transverse impact parameter significance of less than 3.0, and a weighted longitudinal impact parameter |z 0 ×sin θ| of less than 0.5 mm. Reconstructed muons are required to have p T > 25 GeV.
Photons and electrons are required to meet identification criteria based on shower shapes in the EM calorimeter, leakage into the hadronic calorimeter, and ID tracking information. The resulting selected photons are classified as "loose" or "tight" and the electrons as "medium" as defined in Refs. [35,36,39]. The "tight" identification criterion for photons is used to suppress the background from multiple showers produced in meson (e.g., π 0 , η) decays [36]. The electron identification criteria are used to suppress background electrons (primarily from photon conversions and Dalitz decays) and jets faking electrons [37].
Photons, electrons, and muons are required to be isolated from nearby hadronic activity. Photons are considered isolated if the sum of transverse energy calculated from clusters of calorimeter energy deposits [40] in an "isolation" cone of size ∆R = 0.4 around the candidate, E iso T , is smaller than 4 GeV after subtracting the contribution from the photon itself, and corrected for the leakage of the photon energy and the effects of underlying event and pileup [41,42]. For electrons to be isolated, the calorimeter transverse energy deposits and the sum of the transverse momenta of tracks in a cone of size ∆R = 0.2 around the candidate after subtracting the contribution from the electron itself must be below 0.14× p e T and 0.13× p e T , respectively, where p e T is the electron transverse momentum. Muons are considered isolated if the sum of the transverse momenta of ID tracks excluding the track associated with the muon in a cone of size ∆R = 0.2 is below 0.1 × p µ T , where p µ T is the muon transverse momentum. All lepton and photon efficiencies of the trigger, reconstruction, and identification are corrected in the simulation with data-derived correction factors.
Jets are reconstructed from clustered energy deposits in the calorimeter using the anti-k t algorithm [43] with radius parameter R = 0.4 and are required to have p T > 30 GeV and |η| < 4.5. Reconstructed calorimeter jets are corrected for effects of noncompensating response, energy losses in the dead material, shower leakage, as well as inefficiencies in energy clustering and jet reconstruction by applying a simulation-based correction derived in bins of η and E. An in situ calibration corrects for differences between data and simulation in the jet response. This jet energy scale calibration is thoroughly discussed in Ref. [44]. In order to reduce pileup effects, for jets with p T < 50 GeV and |η| < 2.4 the jet vertex fraction (JVF), defined as the ratio of the summed scalar p T of tracks associated with both the R = 0.4 jet and the primary vertex to that of all tracks associated with the jet, must be greater than 0.5.
To reject electrons reconstructed from a bremsstrahlung photon emitted by a muon traversing the calorimeter, any electron candidate within a ∆R = 0.1 cone around a selected muon is removed. Jets are removed if they are found within a ∆R = 0.3 cone around a selected lepton or photon.
The missing transverse momentum vector p miss T is the vector of momentum imbalance in the transverse plane. The reconstruction of the direction and magnitude of the missing transverse momentum vector is described in Ref. [45]. The p miss T is calculated from the vector sum of the calibrated transverse momenta of all jets with p T > 20 GeV and |η| < 4.5, the transverse momenta of electron and muon candidates, and all calorimeter energy clusters not belonging to a reconstructed object (soft-term). Selection criteria based on p miss T or its magnitude E miss T are used only in the neutrino channels, as described in Section 4.3.
4.2 Selection of + − γ and + − γγ event candidates Selected + − γ or + − γγ event candidates must contain exactly one pair of same-flavor, opposite-charge isolated leptons (electrons or muons) and at least one or two isolated photons with E γ T > 15 GeV, respectively. In the case of additional photon candidates, those with the highest E γ T are selected. The dilepton invariant mass m + − is required to be greater than 40 GeV. The reconstructed photons are removed if they are found within a ∆R = 0.7 (0.4) cone around a selected lepton for + − γ ( + − γγ) events. A further requirement on the photon-photon separation of ∆R(γ, γ) > 0.4 is applied in + − γγ events. The selected events are categorized as inclusive events, referring to those with no requirement on the jets, and exclusive events, which are defined to be those with no selected jet with p T > 30 GeV and |η| < 4.5.

Selection of ννγ and ννγγ event candidates
The ννγ event candidates are selected by considering events with E miss T > 100 GeV and at least one isolated photon with E γ T > 130 GeV. The separation between the reconstructed photon direction and p miss T in the transverse plane is required to be ∆φ( p miss T , γ) > π/2, since in signal events the Z boson should recoil against the photon. The ννγγ event candidates are selected by considering events with E miss T > 110 GeV and at least two isolated photons with E T > 22 GeV and ∆R(γ, γ) > 0.4. The directions of the diphoton system and the p miss T are required to be separated in the transverse plane by ∆φ( p miss T , γγ) > 5π/6. In the case of additional photon candidates in ννγ/ννγγ events, one/two photons with the highest E γ T are selected. To suppress W(γ)+jets and Wγ(γ) backgrounds, events containing an identified muon or electron (as defined in Section 4.1 without isolation requirement) are rejected. The selected events are categorized as inclusive events and exclusive events, as described in Section 4.2.

Estimation of backgrounds
This section describes the background estimation in each of the final states. The backgrounds in the + − γ and + − γγ final states are discussed in Section 5.1. The dominant backgrounds in these final states are Z+jets and Zγ+jets with jets misidentified as photons. The backgrounds in the ννγ and ννγγ final states are discussed in Section 5.2. The dominant backgrounds in these final states are those with jets misidentified as photons, those with electrons misidentified as photons, as well as W( ν)γ and W( ν)γγ where the lepton from the W decay is not detected. 5.1 Backgrounds to + − γ and + − γγ Backgrounds in the selected + − γ and + − γγ samples are dominated by events in which hadronic jets, which contain photons from π 0 or η decays, are misidentified as prompt photons. In the + − γ measurement, the background from jets misidentified as photons originates from the production of Z bosons in association with jets (Z+jets), while in the + − γγ measurement this background originates from both Zγ in association with jets (Zγ+jets) and Z+jets events with one or two jets misidentified as photons, respectively. The backgrounds from jets misidentified as photons are estimated using data-driven methods as described in Sections 5.1.1 and 5.1.2. Smaller backgrounds originate from ttγ, WZ, and τ + τ − γ for + − γ, and from WZ, ZZ, and τ + τ − γγ for + − γγ. The backgrounds from ttγ and τ + τ − γ(γ) yield the same final states as the signals, while the backgrounds from WZ and ZZ meet the selection criteria when the electrons from the W or Z decay are misidentified as photons or when final-state photons are radiated. These are expected to contribute in total less than 1.5% of the selected event yield in both the + − γ and + − γγ final states, and are derived from simulation as described in Section 5.1.3.

Estimation of the background from jets misidentified as photons in + − γ measurements
For the + − γ measurement, a two-dimensional sideband method is used to measure the background from jets misidentified as photons, as described in Refs. [8,41]. In this method, a looser photon selection is considered, in which the isolation and some identification requirements on the photon are discarded. After this selection, the + − γ events are separated into one signal and three control regions, defined by varying the photon identification and isolation requirements. Photon candidates failing a subset of requirements on the photon shower-shape variables but satisfying all other requirements in the "tight" photon identification are considered as "nontight". Events in the signal region (A) have the photon satisfying the nominal photon isolation and "tight" identification requirements as described in Section 4.1. The three control regions are defined as: i. Control region B: the photon candidate meets the "tight" identification criteria and is not isolated (E iso T > 4 GeV); ii. Control region C: the photon candidate meets the "nontight" identification criteria and is isolated (E iso T < 4 GeV); iii. Control region D: the photon candidate meets the "nontight" identification criteria and is not isolated (E iso T > 4 GeV).
The shower-shape requirements that the "nontight" photons are required to fail are chosen to enhance the Z+jets background events in the control regions while minimizing the correlation with the photon isolation. The number of Z+jets events in the signal region, N j→γ A , can be derived from the number of observed events in the control regions N i (i=B,C,D) : The coefficients c i (i=B,C,D) are equal to the ratio of the + − γ yields in the control regions to the signal region, and are estimated from simulation. The R factor accounts for a potential correlation between the photon identification and isolation variables for the Z+jets background. The central value of R is taken to be one, as would be the case for no correlation. Its uncertainty of 20% is determined by the deviation of the R value from one as determined from simulation studies of the Z+jets background. The yields The uncertainty in the value of R represents the dominant systematic uncertainty of 24% in the estimate of the Z+jets background. The second largest systematic uncertainty of 10% arises from the inaccuracy in modeling of the coefficients c i , mainly due to the uncertainties in photon identification and isolation efficiencies. An additional Z+jets background uncertainty of 5% arises from uncertainties in the estimates of the N Other BKG i in each of the control regions.

Estimation of the background from jets misidentified as photons in + − γγ measurements
A matrix method as described in Ref. [46] is used to estimate the background from jets misidentified as photons in + − γγ events from Z( + − )γ+jets and Z( + − )+jets events with one or two jets misidentified as photons. The method uses as inputs the jet-to-photon misidentification rate (fake rate), f , which is the probability for a jet satisfying "loose" photon identification criteria [36] to be identified as a "tight" and isolated photon, and the real photon identification efficiency, , which is the probability for "loose" prompt photons to be identified as "tight" and isolated photons. The fake rate and the real photon identification efficiency are estimated from data and from MC simulation, respectively. A 4×4 matrix is constructed from the fake rate and the real photon identification efficiency, relating the observed number of events, N TT , N TL , N LT , N LL , to the unknown number of each type of event, N γγ , N γjet , N jetγ , N jetjet , by a set of linear equations: In the subscripts TT, TL, LT, LL, the first (second) subscript refers to the leading (subleading) reconstructed photon candidate; T means that it is "tight" and isolated while L corresponds to a "loose", not "tight" or not isolated candidate. Similarly, the subscripts γγ, γjet, jetγ, and jetjet correspond to the cases of two photons, leading photon and subleading jet, leading jet and subleading photon, and two jets, respectively. The subscripts 1 and 2 refer to the leading and subleading photon candidates, respectively. The number of each type of event, N γγ , N γjet , N jetγ , N jetjet , is obtained by solving Equation (3), from which the number of background events with jets misidentified as photons in the signal region, N j→γ TT , is then obtained: The fake rate is estimated from data using a sample enriched in Z( + − )+jets with one jet misidentified as a photon. To suppress the contribution from Z → + − γ, the invariant mass of opposite-charge dilepton pairs in the events is required to be within 8 GeV of the Z boson mass. A two-dimensional sideband method similar to that described in Section 5.1.1 is used to estimate the number of + − +jets events in which the "loose" jets satisfy the "tight" identification and isolation requirements. As the fake rate depends on the photon E T , a fake rate as a function of the photon E T is used in the matrix method. The real photon identification efficiency, which is also a function of the photon E T , is estimated from MC simulation.
The systematic uncertainty related to the background from jets misidentified as photons is dominated by the potential bias of the two-dimensional sideband method to estimate the fake rate. It is evaluated from Z+jets MC simulation to be about 23%, by comparing the fake rate calculated by the two-dimensional sideband method to the fake rate calculated using the generator-level information in the MC simulation. Other systematic uncertainties, arising from possible inaccuracy in modeling of the real photon identification efficiency, other electroweak backgrounds, as well as the dependence of and f on photon η, sum to about 10%.

Results of the background estimation for + − γ and + − γγ
The backgrounds other than those from jets misidentified as photons are estimated using MC simulation. The systematic uncertainties in these backgrounds consist of the experimental uncertainties described in Section 6.2 and the cross-section uncertainties, which are 22% (ttγ [34]), 10% (WZ [47,48]) and 15% (ZZ [47,49]). The cross-section uncertainties in the τ + τ − γ and τ + τ − γγ backgrounds are evaluated to be 7% using MCFM, as described in Section 7.1. An additional uncertainty of 30% (60%) is assigned to the WZ (ZZ) background to account for the mismodeling of the electron-to-photon fake rate. This uncertainty is estimated by comparing the fake rate predicted by simulation to that estimated in data, using the method described in Section 5.2.3.
The number of events observed in data, N obs Zγ , as well as the estimated background yields in the + − γ and + − γγ measurements, are summarized in Tables 1 and 2, respectively.
The E T distributions of photons selected in the e + e − γ and µ + µ − γ inclusive measurements are shown in Figure 2. The highest-E T photon is measured as E  ) in the e + e − γ and µ + µ − γ channels with the inclusive (N jets ≥ 0) and exclusive (N jets = 0) selections. The first uncertainty is statistical and the second is the sum of all contributions to the systematic uncertainty. The statistical uncertainties arise from the numbers of events in the control regions and the simulation. The systematic uncertainties in the signal include both the experimental uncertainties described in Section 6.2 and the theoretical uncertainties in the cross sections evaluated using MCFM, as described in Section 7.1.

Backgrounds to ννγ and ννγγ
Backgrounds to the ννγ and ννγγ signals originate from several sources (listed in decreasing order of significance): events with prompt photons and mismeasured jet momenta causing missing transverse momentum (dominant for the inclusive measurement); nonsignal electroweak processes, such as W( ν)γ, with partial event detection; events with real E miss T from neutrinos (such as Z(νν) or W(eν)) and misidentified photons from electrons or jets. The largest contributions are determined using data-driven techniques.
25.7 ± 0.5 ± 1.6 29.5 ± 0.6 ± 1.7 18.9 ± 0.5 ± 1.5 21.8 ± 0.5 ± 1.7 ) in the e + e − γγ and µ + µ − γγ channels with the inclusive (N jets ≥ 0) and exclusive (N jets = 0) selections. The first uncertainty is statistical and the second is the sum of all contributions to the systematic uncertainty. The statistical uncertainties arise from the numbers of events in the control regions and the simulation. The systematic uncertainties in the signal include both the experimental uncertainties described in Section 6.2 and the theoretical uncertainties in the cross sections evaluated using MCFM, as described in Section 7.1.   The photon transverse energy (E γ T ) distributions from inclusive (N jet ≥ 0) + − γ events for the electron (left) and muon (right) channels. The numbers of candidates observed in data (points with error bars) are compared to the sum of the SM signal predicted from Sherpa and the various backgrounds discussed in Section 5.1. The uncertainty band on the sum of expected signal and backgrounds includes both the statistical and systematic uncertainties in the MC simulations and the data-driven background estimate added in quadrature. The signal is normalized using the cross sections predicted by Sherpa. The theoretical uncertainties in the signal cross sections are evaluated bin-by-bin using MCFM, as described in Section 7.1. The ratio of the numbers of candidates observed in data to the sum of expected signal and backgrounds is also shown. 18 -1 = 8 TeV, 20.3 fb s Figure 3: The four-body invariant mass (m + − γγ ) distributions from inclusive (N jet ≥ 0) + − γγ events for the electron (left) and muon (right) channels. The numbers of candidates observed in data (points with error bars) are compared to the sum of the SM signal predicted from Sherpa and the various backgrounds discussed in Section 5.1. The uncertainty band on the sum of expected signal and backgrounds includes both the statistical and systematic uncertainties in the MC simulations and the data-driven background estimate added in quadrature. The signal is normalized using the cross sections predicted by Sherpa. The theoretical uncertainties in the signal cross sections are evaluated using MCFM, as described in Section 7.1.  Figure 4: The diphoton invariant mass (m γγ ) distributions from inclusive (N jet ≥ 0) + − γγ events for the electron (left) and muon (right) channels. The numbers of candidates observed in data (points with error bars) are compared to the sum of the SM signal predicted from Sherpa and the various backgrounds discussed in Section 5.1. The uncertainty band on the sum of expected signal and backgrounds includes both the statistical and systematic uncertainties in the MC simulations and the data-driven background estimate added in quadrature. The signal is normalized using the cross sections predicted by Sherpa. The theoretical uncertainties in the signal cross sections are evaluated using MCFM, as described in Section 7.1.

ATLAS
The procedures used to estimate these backgrounds follow closely those in a previous ATLAS measurement [8]. Smaller backgrounds originate from τ + τ − γ for ννγ and τ + τ − γγ for ννγγ. These are expected to contribute less than 1.5% of the selected event yield and are derived from MC simulation. The backgrounds from multijet and + − γ processes are negligible. Each source of background is discussed in detail together with the method used for its estimation in the following subsections.

γ+jets background to ννγ
An imprecise measurement of jet activity in the calorimeter can cause the appearance of fake E miss T in the event. Photon+jets events are one of the dominant background contributions to the ννγ channel. Although the high-E miss T requirement reduces the γ+jets background, a residual contamination from this background remains for the inclusive measurement and is estimated with the following data-driven method.
In order to measure this background from data, a control sample enriched in γ+jets events is selected by applying all the signal region (SR) selection criteria, but inverting the angular separation requirement such that ∆φ( p miss T , γ) < π/2. The data yield in this control region (CR), after subtraction of signal and other backgrounds obtained using the MC simulation, is then extrapolated to the signal region with a transfer factor determined from a γ+jets simulation. The transfer factor equals the ratio of the numbers of γ+jets events in the SR to the CR. The nominal transfer factor is determined to be 1.1 from Sherpa and a 30% uncertainty is estimated using an alternative prediction from Pythia.

W( ν)γ background to ννγ
Misidentified events from W( ν)γ production are one of the dominant background contributions to the ννγ signal. A large fraction (about 60%) of this contamination originates from W(τν)γ events. A scale factor is defined to correct the yield of Wγ events estimated by MC simulation to match the Wγ event yield measured in a control data region constructed by requiring exactly one identified electron or muon instead of the charged-lepton veto. Since the control region contains some amount of signal leakage and other background contaminations, these contributions are estimated using the methods described in Sections 5.2.1, 5.2.3 as well as with MC simulation and then subtracted. With equal branching fraction of the W boson leptonic decays, the MC scale factor for the dominant W(τν)γ events in the signal region and its uncertainty are taken from the measurement of W( ν)γ events in the control region. The main uncertainty of 34% in this background prediction is due to the extrapolation transfer factor from the control region to the signal region. This is estimated by comparing transfer factors between two MC samples generated with Sherpa and Alpgen, respectively. The transfer factor between the control and the signal regions is taken from Sherpa as the baseline and equals 2.2 ± 0.7 for the inclusive selection and 1.8 ± 0.7 for the exclusive selection.

W(eν) background to ννγ
Misidentification of electrons as photons also contributes to the background yield in the signal region. The estimation of this background is made in two steps. The first is the determination of the probability for an electron to be misidentified as a photon using Z(e + e − ) decays reconstructed as e + γ, as described in Ref. [50]. The probability of observing an e + γ pair with invariant mass near the Z boson mass is used to determine an electron-to-photon fake factor f e→γ . This increases from 2% to 6% as |η| increases from 0 to 2.37. The second step is the construction of a control region with nominal ννγ selection criteria, except that an electron is required instead of the photon in the final state. This control region contains W(eν)+jets as the dominant process and some fractions of other processes containing genuine electrons and jets. The estimated W(eν) background is then the product of the electron-to-photon fake factor by the number of events in the chosen control sample. The total uncertainty in this background varies from 10% to 30% as a function of photon E T and η and is dominated by the number of events in the e + γ control sample used to measure the electron misidentification probability.

Z(νν)+jets backgrounds to ννγ
Misidentification of jets as photons gives a non-negligible background contribution to the ννγ signal. A data-driven method similar to the one described for Z( + − )+jets in Section 5.1.1 is used to determine the background contribution from Z(νν)+jets events. A systematic uncertainty of 25% in this background is assigned, dominated by the uncertainty in the correlation factor between identification and isolation of jets reconstructed as photons.

γ+jets and γγ+jets backgrounds to ννγγ
The estimation of γ+jets and γγ+jets backgrounds to the Z(νν)γγ signal uses a two-dimensional sideband method. Four regions are constructed using two orthogonal selections: different E miss T requirements (E miss T < 20 GeV or E miss T > 110 GeV) and different identification requirements for photons (two "tight" photons or one "tight" photon and one photon meeting the looser criteria but not the "tight" ones). Since the correlations between these regions are small, the number of background events in the signal region can be estimated by scaling the number of events in the high-E miss T control region by the ratio of the events from control samples in the low E miss T region. Corrections are applied for the Z(νν)γγ signal and other backgrounds leaking into the control samples. The largest uncertainty in this procedure is due to the number of events in the control regions. Systematic uncertainties for this background are evaluated with alternative low E miss T control regions (5 < E miss T < 25 GeV) and from the uncertainty in the correlation between control regions (15%).

W( ν)γγ background to ννγγ
The background from W( ν)γγ events is dominated by the τνγγ contribution and is estimated using techniques similar to those described above in Section 5.2.2. A control region is defined by requiring exactly one identified electron or muon instead of the charged-lepton veto. After accounting for signal leakage and other background contributions, the control region yield is compared to the Wγγ simulation. Good agreement is found, as in the recent measurement of the Wγγ cross section [51], although in the high-E miss T region considered here the size of the control sample leads to a 100% uncertainty in the transfer factor.

W(eν)γ background to ννγγ
One of the dominant backgrounds in the ννγγ channel originates from the misidentification of electrons as photons. This background is estimated by selecting a control sample in which an electron is required instead of one of the photons in the ννγγ final state. The electron fake rate is estimated as described in Section 5.2.3. The estimated background in the signal region is then obtained by rescaling the yield in the control sample by the electron-to-photon fake rate. The largest uncertainty in this background is 20% and is derived from MC events in a closure test of the method.

Z(νν)γ+jets background to ννγγ
The Z(νν)γ+jets background falls into the signal region when one jet is misidentified as a photon. This background contributes less than 5% of the total event yield and is estimated from the MC simulation. The systematic uncertainty arises from the mismodeling of the jet-to-photon misidentification rate in the MC simulation. It is evaluated to be 127% (106%) in the inclusive (exclusive) channel, based on Z( + − )γ+jets events with one jet misidentified as a photon, by comparing its estimate from data (as described in Section 5.1.2) with the prediction from MC simulation.

Results of the background estimation for ννγ and ννγγ
A summary of the number of events observed in data and the background contributions in the ννγ(γ) channels is given in Tables 3 and 4. The photon transverse energy and the missing transverse momentum distributions from the selected events in the ννγ channel are shown in Figure 5. The highest-E T photon is measured as E The number of signal events in each of the four production channels, + − γ, ννγ, + − γγ, and ννγγ, is determined by subtracting the estimated backgrounds from the number of observed events. The signal yields are then corrected for detection efficiencies in the fiducial regions used for the measurements. The cross sections are calculated for slightly extended fiducial regions using SM predictions for the extrapolation. These cross sections allow a combination of data obtained from the Z boson to electron and muon decay channels and are more easily compared to predictions from theory. The extended fiducial regions (see Table 5) are defined at the particle level, as described below. The methods used for the determination of the cross sections and their uncertainties are described in Section 6.2. The integrated and differential cross-section measurement results are presented in Sections 6.3 and 6.4, respectively.
"Particle level" refers to stable particles with a proper decay length cτ > 10 mm which are produced from the hard scattering, including those that are the products of hadronization. The fiducial regions are defined with the same object and event kinematic selection criteria as the reconstruction-level selections described in Section 4. Compared with the fiducial regions, the extended fiducial regions use a unified charged Table 3: Total number of events satisfying the ννγ selection requirements in data (N obs Zγ ), predicted number of signal events from Sherpa (N sig Zγ ), and the expected number of background events for each of the sources and together (N bkg Zγ ) with the inclusive (N jets ≥ 0) and exclusive (N jets = 0) selections. The first uncertainty is statistical and the second is the sum of all contributions to the systematic uncertainty. The statistical uncertainties arise from the numbers of events in the control regions and the simulation. The systematic uncertainties in the signal include both the experimental uncertainties described in Section 6.2 and the theoretical uncertainties in the cross sections evaluated using MCFM, as described in Section 7.1.  Table 4: Total number of events satisfying the ννγγ selection requirements in data (N obs Zγγ ), predicted number of signal events from Sherpa (N sig Zγγ ), and the expected number of background events for each of the sources and together (N bkg Zγγ ) with the inclusive (N jets ≥ 0) and exclusive (N jets = 0) selections. The first uncertainty is statistical and the second is the sum of all contributions to the systematic uncertainty. The statistical uncertainties arise from the numbers of events in the control regions and the simulation. The systematic uncertainties in the signal include both the experimental uncertainties described in Section 6.2 and the theoretical uncertainties in the cross sections evaluated using MCFM, as described in Section 7.1.   lepton pseudorapidity selection criterion |η | < 2.47 for + − γ and + − γγ channels. As for ννγ and ννγγ channels, the extended fiducial regions remove the ∆φ( p miss T , γ) > π/2 and ∆φ( p miss T , γγ) > 5π/6 requirements, respectively. Final-state radiation is incorporated into the particle level definition of the leptons by including the contributions from the photons within a cone of ∆R = 0.1 around the lepton direction. The particle level jets are reconstructed using the anti-k t algorithm with a radius parameter of R = 0.4, including all stable particles except for muons and neutrinos. The photons at particle level are required to satisfy the isolation criterion of p h < 0.5, where p h is the transverse energy carried by the closest particle-level jet in a cone of ∆R = 0.4 around the photon direction, subtracting the photon E T and then divided by the photon E T .

Determination of extended fiducial cross sections
The integrated cross sections for Zγ and Zγγ production in the extended fiducial regions are calculated using : where N is the number of candidate events observed, B is the expected number of background events and Ldt is the integrated luminosity corresponding to the dataset analyzed. The factors C and A correct for detection efficiency and acceptance, respectively: [GeV] [GeV] Cuts Inclusive : N jet ≥ 0, Exclusive : N jet = 0 Table 5: Definition of the extended fiducial regions where the cross sections are measured. The variable p νν T is the transverse momentum of the Z boson decaying to a neutrino pair. The variable p h is the transverse energy carried by the closest particle level jet in a cone of ∆R = 0.4 around the photon direction, excluding the photon and divided by the photon transverse energy.
• C is defined as the number of reconstructed signal events satisfying all selection criteria divided by the number of events that, at particle level, meet the acceptance criteria of the fiducial region.
• A is defined as the number of signal events within the fiducial region divided by the number of signal events within the extended fiducial region, which are both defined at particle level.
The corrections A and C are determined using the Zγ and Zγγ signal events generated with Sherpa. The numerical values are summarized in Table 6.  Systematic uncertainties in the acceptances A are evaluated by varying the PDFs and the renormalization and factorization scales. The uncertainty in the acceptances due to the PDF is taken as the envelope of the internal uncertainties from three different PDF sets, namely, the CT10 PDF set, the MSTW2008NLO PDF set [52], and the NNPDF2.3 PDF set [53]. The internal uncertainty from each PDF set is estimated by comparing the acceptance using the PDF central set with the acceptance estimated using the PDF eigenvector sets. The renormalization and factorization scale uncertainties are assessed by varying these two scales independently by a factor of two from their nominal values, and taking the envelope of the resulting variations. The impact of PDF uncertainties varies from 0.04% to 0.3%, while the renormalization and factorization scale uncertainties cause variations from 0.08% to 1.5%. The total uncertainties in the acceptance factors are summarized in Table 6 .
Systematic uncertainties affecting the correction factors C can be grouped into two categories. The first includes the uncertainties arising from the efficiencies of the trigger, reconstruction, identification, and other selection requirements. The second category stems from the uncertainties of energy and momentum scales and resolutions of the final-state objects and the simulation of pileup events. Table 7 presents all the contributions to the uncertainties in C determined using the methods described below. The total uncertainties in the correction factors are summarized in Table 6 .
The photon identification efficiencies are measured in data using a combination of three methods as described in Ref. [36]. The uncertainties induced by the photon identification efficiency are estimated to be 1.5% and 0.5% for the + − γ and ννγ channels, respectively. For the + − γγ and ννγγ channels, after taking into account the correlations between the two photons, the resulting uncertainties are 2.1% and 1.9%, respectively. The photon isolation efficiencies are determined from data by studying the electron isolation efficiencies using Z → e + e − events. The estimated uncertainty increases from 0.5% for photons with E T around 20 GeV to 8% for photons with E T greater than 350 GeV, dominated by the limited size of the Z → e + e − sample in data.
The reconstruction and identification efficiencies of electrons and muons are derived using a tag-andprobe method with Z and J/ψ events decaying into e + e − or µ + µ − pairs [38,39]. The uncertainties are evaluated to be 1.6% for the electron channels, and 0.9% for the muon channels. The uncertainties arising from the selection efficiencies of lepton isolation and impact parameters requirements are also measured with a tag-and-probe method using Z events. They are found to be 2.2%. The uncertainties due to the modeling of trigger efficiencies are evaluated to be 1.9% for the ννγ channel and no more than 0.5% for the other channels [54,55]. The uncertainty in the jet vertex fraction efficiency is estimated by varying the selection requirement to account for the difference between data and simulation. For exclusive N jets = 0 measurements, they are calculated to be no more than 0.6% for all the channels.
The energy scale and resolution and their uncertainties for electrons and photons are obtained using Z → e + e − events [37]. The systematic uncertainty due to the energy scale varies from 1.2% to 2.7% and that associated with the energy resolution is no more than 0.5% for all the final states. The muon momentum scale and resolution are studied using samples of J/ψ, Υ, and Z decays to muon pairs [38]. The corresponding uncertainties are no more than 0.5% in all the channels.
The exclusive N jets = 0 measurements are affected by the uncertainties in the jet energy scale and resolution, because these uncertainties change the distributions of the number of jets with p T > 30 GeV and |η| < 4.5. They are studied using MC simulation, as well as γ+jet, Z+jet, dijet, and multijet data events [44]. Their systematic effect varies from 0.8% to 2.9% for all channels. The uncertainties in the energy and momentum scales and resolutions of reconstructed physics objects are propagated to the E miss T calculation. The uncertainties arising from the scale and resolution of the energy deposits that are not associated with any reconstructed physics object, named the E miss T soft-term [45], are no more than 0.5% for the ννγ final state, and vary from 0.4% to 1.7% for the ννγγ final state. As mentioned in Section 3.1, the MC events are reweighted so that the pileup conditions in the simulation match the data. The pileup events are modeled by MC simulation. The uncertainties associated with the modeling of the pileup events are estimated to be no more than 1.1% for all the final states.  Table 7: Relative systematic uncertainties, in %, in the signal correction factor C for each channel in the inclusive N jets ≥ 0 (exclusive N jets = 0) measurement.

Integrated extended fiducial cross sections for Zγ and Zγγ production
The measurements of the cross sections of each final state and the combined charged-lepton final states, along with their uncertainties, are based on the maximization of the profile-likelihood ratio: where L represents the likelihood function, σ is the cross section and θ are the nuisance parameters corresponding to sources of the systematic uncertainties. Theσ andθ terms denote the unconditional maximum-likelihood estimate of the parameters, i.e., where the likelihood is maximized for both σ and θ. Theθ(σ) corresponds to the value of θ that maximizes L for given parameter values of σ. The likelihood function is defined as : It corresponds to the product of the Poisson probability of observing N i events in each final state, given the expectation for the signal S i and background B i , and is multiplied by the Gaussian constraints on the systematic uncertainties θ with central values θ 0 from auxiliary measurements as described in Section 6.2.
The measured cross sections for the Zγ and Zγγ processes in the extended fiducial regions defined in Table 5 are summarized in Table 8. The theoretical predictions in the table are described in Section 7. The significance for the combination of e + e − γγ and µ + µ − γγ processes is 6.3 (6.0) standard deviations for the inclusive (exclusive) selection.
The Zγ inclusive (exclusive) cross sections in the extended fiducial regions are measured with a precision of 6% (6%) in the + − γ final state and 50% (24%) in the ννγ final state. The smaller uncertainty in the exclusive ννγ measurement results from the reduced background fraction as shown in Table 3. The Zγγ inclusive (exclusive) cross sections in the extended fiducial regions are measured with a precision of 16% (19%) in the + − γγ final state and 70% (60%) in the ννγγ final state. The precision of the Zγ cross-section measurements is driven by their systematic uncertainties. For the Zγγ cross sections, the precision of the measurements is dominated by the statistical uncertainty in the + − γγ final state, and is equally affected by statistical and systematic uncertainties in the ννγγ final state.
The systematic uncertainties in the measured cross sections in Table 8 arise from the uncertainties in the acceptances A and correction factors C, as well as from the uncertainties in the estimates of backgrounds.
In the + − γ and + − γγ final states the two sources have effects of comparable size on the measured cross sections, while in the ννγ and ννγγ final states the uncertainties in the estimates of backgrounds dominate.
Compared with the Zγ measurements at √ s = 7 TeV [8], the systematic uncertainty is reduced in the + − γ final state while it becomes larger in the ννγ final state. The reduced systematic uncertainty in the + − γ final state mainly results from the reduced systematic uncertainty from photon identification efficiency, as well as the smaller statistical uncertainty in the data-driven estimate of the Z+jets background. The larger systematic uncertainty in the ννγ final state is largely a result of the increased photon E T threshold requirement due to the increased single-photon trigger E T threshold, which results in generally increased systematic uncertainties in the estimates of backgrounds.
The measurements of the cross sections in the e + e − γ and µ + µ − γ final states agree within one standard deviation. In order to assess the compatibility of the cross-section measurements in the e + e − γγ and µ + µ − γγ final states, a profile-likelihood ratio is constructed, parameterized as a function of the difference in measured cross sections. With this approach, the measurements are found to be compatible within 1 Table 5. The SM predictions from the generator MCFM calculated at NLO, as well as the predictions at NNLO [56] (for Zγ only), are also shown in the table with combined statistical and systematic uncertainties. All MCFM [57] and NNLO predictions are corrected to particle level using parton-to-particle scale factors as described in Section 7.1.

Differential extended fiducial cross section for Zγ production
The measurements of differential cross sections allow the comparison of data results to theory predictions in terms of not only their overall normalizations, but also their shapes. The measurements are performed for Zγ production in several observables that are sensitive to higher-order perturbative QCD corrections.
These include the photon transverse energy E γ T , the invariant mass of the + − γ three-body system, and the jet multiplicity N jets . The differential cross sections are defined in the extended fiducial region, and are extracted with an unfolding procedure to remove measurement inefficiencies and resolution effects from the observed distributions. The procedure described in Ref.
[8] is followed, using an iterative Bayesian method [58]. Events from simulated signal MC samples are used to generate a response matrix for each distribution. Each element of the response matrix is the conditional probability that an event is found in bin i in the measurement given that it is in bin j at the particle level. In the first iteration, the prior distribution of the particle level prediction is given by the signal MC sample. The response matrix and the measured distribution then modify the prior distribution, giving the posterior distribution at the particle level. For each further iteration, the posterior distribution of the previous iteration is used as the new prior distribution. Three iterations are found to be optimal, as too many iterations give rise to large statistical fluctuations, while too few can produce a result that is biased by the dependence on the initial prior distribution.
The statistical uncertainties of the unfolded distribution are estimated using pseudoexperiments, generated by fluctuating each bin of the observed spectrum according to a Poisson distribution with the expected value equal to the observed yield. The shape uncertainties from the number of signal MC events are also obtained by performing pseudoexperiments. The sources of systematic uncertainty are discussed in Section 6.2, with their impact on the unfolded distribution assessed by varying the response matrix for each of the systematic uncertainty sources by one standard deviation and adding up the resulting changes in quadrature. The results from the electron and muon channels are combined with equal weight, taking into account the correlations between the systematic uncertainties in the two channels.
In addition to the systematic uncertainties described in Section 6.2, the differences between the unfolded results with three iterations and the results with two or four iterations are taken as systematic uncertainties associated with the unfolding method.
The differential cross sections are presented as a function of E γ T in Figure 7 for the inclusive and exclusive measurements of the + − γ channel and in Figure 8 for the inclusive and exclusive measurements of the ννγ channel. The differential cross sections are shown in Figure 9 as a function of m + − γ . Figure 10 shows the cross sections in the + − γ channel measured in bins of jet multiplicity. The predictions in the figures are described in Section 7. As with the integrated cross sections shown in Table 8, the differential cross sections of the exclusive measurements in the ννγ channel have smaller uncertainties than the inclusive measurements.

Comparison of measurements to Standard Model predictions 7.1 Estimation of Standard Model expectations
The measurements of Zγ and Zγγ production are compared to SM predictions using the parton shower Monte Carlo Sherpa 1.4 and the NLO parton-level generator MCFM. 2 In addition, parton-level NNLO SM predictions for Zγ are compared to data using the calculations described in Ref. [56]. The theory predictions include off-shell Z bosons and direct photons arising from initial-state radiation (from the quarks) and radiative Z-boson decay in the case of charged-lepton final states, and from fragmentation of final-state quarks and gluons into photons, leading to the production channels pp → + − γ(γ) + X and pp → ννγ(γ) + X. In the Sherpa and MCFM generators, contributions from quark/gluon fragmentation into isolated photons are also included. The CT10 PDF set [18] is used for the Sherpa and MCFM generation, and the MMHT2014 PDF set [59] is used for the NNLO predictions. The renormalization and factorization scales are set equal to m Zγ (m Zγγ ) for the MCFM NLO generation of Zγ (Zγγ) events and to m 2 Z + (E γ T ) 2 for the NNLO Zγ predictions. The other electroweak parameters used are the default values [60] from the authors of the generators.
The events generated with Sherpa as described in Section 3.1 are also compared to the measurements at particle level. For the NLO and NNLO parton-level predictions, parton-to-particle correction factors C * (parton → particle) must be applied in order to obtain the particle level cross sections. These correction factors are computed as the ratios of the pp → Zγ(γ) cross sections predicted by Sherpa with hadronization and the underlying event disabled to the cross sections with them enabled. The systematic uncertainties in the correction factors are evaluated by using an alternative parton-showering method [61] within Sherpa, and are found to be negligible compared to the statistical uncertainties. The particle level cross sections are obtained by dividing the NLO and NNLO parton-level predictions by the C * (parton → particle) correction factors summarized in Table 9. The corrections are a few percent for the inclusive cross sections and reach about 10% for some exclusive channels. The correction factors in Table 9 apply to the predictions made for the Zγ and Zγγ cross sections in the extended fiducial region described in Table 5. The systematic uncertainties in the SM NLO cross sections are estimated by varying the QCD scales by factors of 0.5 to 2.0 (independently for the renormalization and factorization scales) and varying the CT10 PDFs by their uncertainties at 68% confidence level. The uncertainties due to the contribution of photons from fragmentation of quarks or gluons are estimated by varying the fraction of hadronic energy p h in the isolation cone from 0.25 to 0.75. For the NLO exclusive zero-jet cross sections the method suggested in Ref. [62] is used to estimate the additional uncertainty due to the N jet = 0 requirement. The systematic uncertainties in the SM NNLO cross sections are determined as described in Ref. [56]. In all cases the uncertainties in the parton-to-particle correction factors are included.
1.01708 ± 0.00065 0.96809 ± 0.00078 ννγ 0.9987 ± 0.0025 0.9150 ± 0.0030 + − γγ 1.0273 ± 0.0039 0.9755 ± 0.0047 ννγγ 1.0012 ± 0.0076 0.873 ± 0.010 Table 9: Parton-to-particle correction factors C * (parton → particle) obtained from the Sherpa MC samples. For + − γ and + − γγ channels the parton-to-particle level correction factors are the weighted average over both lepton flavors (e, µ). The uncertainties include both the statistical and systematic contributions. The measured (points with error bars) and predicted differential cross sections as a function of m + − γ for the pp → + − γ process in the inclusive N jets ≥ 0 (left) and exclusive N jets = 0 (right) extended fiducial regions. The error bars on the data points show the statistical and systematic uncertainties added in quadrature. The MCFM and NNLO predictions are shown with shaded bands that indicate the theoretical uncertainties described in Section 7.1. The Sherpa predictions are shown with shaded bands indicating the statistical uncertainties from the size of the MC samples. The lower plots show the ratios of the predictions to data (shaded bands). The error bars on the points show the relative uncertainties of the data measurements themselves. The bin size varies from 10 GeV to 1360 GeV.

Extended fiducial cross sections compared to SM predictions
The measured extended fiducial cross sections for pp → + − γ+X and pp → ννγ + X production are compared to SM predictions in Table 8. The estimates of the cross section at NLO and NNLO and their systematic uncertainties are obtained as described above. Predictions are made for both inclusive production (no restriction on the system recoil X) and exclusive production of events having no central (|η| < 4.5) jet with p T > 30 GeV. There is generally good agreement between the cross-section measurements for these Zγ channels and the SM predictions; the NNLO calculation of the inclusive cross section for the Z( + − )γ channel gives better agreement with the measurement than the NLO calculation.
Requiring two photons with E T > 15 GeV results in a + − γγ cross section a factor of approximately 400 times smaller than + − γ production. The measurements for both the + − γγ and ννγγ channels are compared to the NLO MCFM predictions in Table 8. The measurements in these channels are statistically limited, but the data are consistent with the predicted SM cross sections. The measured cross sections and the MCFM predictions are compatible within 1.7 (0.9) standard deviations in the inclusive (exclusive) + − γγ channel, and within 1.2 standard deviations in the ννγγ channel.

Differential cross sections compared to SM predictions
The background-subtracted, unfolded differential cross sections for the E γ T spectra from pp → + − γ + X and pp → ννγ + X production are compared to SM expectations in Figures 7 and 8 pp → + − γ + X the NLO calculation underestimates the production of photons at high E T , whereas the NNLO calculation and the Sherpa shower MC both agree with the data. For exclusive pp → + − γ + X production all three SM calculations are in good agreement with the data, as are the SM predictions for the photon E T spectra from pp → ννγ + X production.
The differential spectra of the + − γ invariant mass from pp → + − γ + X are compared to data in Figure 9. For the exclusive channel all three SM predictions agree well with the data. For the inclusive channel the NLO prediction underestimates the cross section at high m + − γ , while the NNLO calculation is in good agreement with the data.
In Figure 10 the measured jet multiplicity spectrum from + − γ events is compared to NLO MCFM predictions for zero and one jet, and to Sherpa for zero to three jets. These SM predictions are in agreement with the data.

Limits on triple and quartic gauge-boson couplings 8.1 Anomalous triple gauge-boson couplings ZZγ and Zγγ
Within the Standard Model, vector-boson self-interactions are completely fixed by the model's SU(2) L × U(1) Y gauge structure [63]. Their observation is thus a crucial test of the model. Any deviation from the SM prediction is called an anomalous coupling. Anomalous triple gauge-boson couplings for Zγ production can be parameterized by four CP-violating (h V 1 , h V 2 ) and four CP-conserving (h V 3 , h V 4 ) complex parameters (where V = Z, γ). All of these parameters are zero at tree level in the SM. Since the CPconserving couplings h V 3,4 and the CP-violating couplings h V 1,2 do not interfere and their sensitivities to aTGCs are nearly identical [63], the limits from this study are expressed in terms of the CP-conserving parameters h V 3,4 . The yields of Zγ events with high E γ T with the exclusive zero-jet selection are used to set the limits. The exclusive selection is used since it significantly reduces the SM contribution at high E γ T and therefore optimizes the sensitivity to anomalous couplings. The contributions from aTGCs increase with the E T of the photon, and the search is optimized to have the highest sensitivity by using the extended fiducial cross sections for Zγ production with E γ T greater than 250 GeV for + − γ and greater than 400 GeV for ννγ. The neutrino channel has the highest sensitivity to aTGCs. The measured cross sections and the SM predictions in these high-E γ T phase-space regions (aTGC regions) are shown in Table 10  for the channels studied. The E γ T threshold is 250 GeV for the electron and muon channels and is 400 GeV for the neutrino channel. The first uncertainty is statistical, the second is systematic.

Form factors (FF) are introduced to avoid unitarity violation at very high parton center-of-mass energy
with the form factor exponent n set to three for h V 3 and four for h V 4 to preserve unitarity [64], where Λ FF is the approximate energy scale at which contributions from physics beyond the SM would become directly observable. The dependencies of the unitarity bounds on the aTGC parameters from the scale Λ FF calculated as in Ref. [65] are shown in Figures 11 and 12, where the observed and expected limits are derived as discussed below. A form factor with Λ FF = 4 TeV is chosen as the lowest scale to preserve unitarity for all the studied parameters. The limits on aTGCs are also given without a form factor (Λ FF = ∞) as a benchmark, although unitarity is not preserved in this case.  The cross-section predictions with aTGCs (σ aTGC Zγ ) are obtained from the MCFM generator. The number of expected Zγ events in the exclusive aTGC region (N aTGC The anomalous couplings influence the kinematic properties of the Zγ events and thus the corrections for event reconstruction (C Zγ ). The maximum variations of C Zγ due to nonzero aTGC parameters within the measured aTGC limits are quoted as additional systematic uncertainties. Since the influence of the anomalous couplings on the acceptance corrections (A Zγ ) and parton-to-particle (C * (parton → particle)) corrections is an order of magnitude smaller than on C Zγ , it is neglected.
The limits on a given aTGC parameter are extracted from a frequentist profile-likelihood test, as explained in Section 6.3. The profile likelihood combines the observed number of exclusive Zγ candidate events for the E γ T threshold mentioned above, the expected signal as a function of aTGC as described in Equation 7, and the estimated number of background events separately for each channel. A point in the aTGC space is accepted (rejected) at the 95% confidence level (C.L.) if fewer (more) than 95% of the randomly generated pseudoexperiments exhibit larger profile-likelihood-ratio values than that observed in data. A pseudoexperiment in this context is a set of randomly generated numbers of events, which follow the Poisson distribution with the mean equal to the sum of the number of expected signal events and the estimated number of background events. The systematic uncertainties are included in the likelihood function as nuisance parameters with correlated Gaussian constraints, and all nuisance parameters are fluctuated in each pseudoexperiment.
The allowed ranges for the anomalous couplings are shown in Table 11 for ZZγ (h Z 3 and h Z 4 ) and Zγγ (h  limits versus Λ FF is shown in Figures 12 and 11. The obtained observed limits are almost a factor of two better than the expected limits, which is due to a downward fluctuation in the region of high E γ T for the ννγ channel. All anomalous couplings considered are found to be compatible with the SM value zero. The observed limits on h  Table 11. These limits are the most stringent to date. The limits on all possible combinations of each pair of aTGC are also evaluated by the same method. The 95% C.L. regions in two-parameter aTGC space are shown as contours on the (h γ 3 , h γ 4 ) and (h Z 3 , h Z 4 ) planes in Figures 14 and 15, since only these pairs are expected to interfere [63].    Since all sensitivity of the measurement of aTGCs is contained in a single measurement of the Zγ cross section in the high-E γ T regions, the likelihood ratio used to obtain the two-parameter limits has one effective degree of freedom. Therefore the results obtained for the aTGC frequentist limits found in the one-parameter fit are identical to the corresponding limits obtained from the two-parameter fits at the points where the other aTGC is zero.

Anomalous quartic gauge-boson couplings ZZγγ and Zγγγ
Triboson Zγγ production in the SM has no contributions from the quartic gauge-boson couplings ZZγγ and Zγγγ. However, physics beyond the SM could induce these anomalous neutral QGCs, enhancing the cross section for Zγγ production and modifying the kinematic distribution of the final-state Z boson and photons. The effect of such new couplings can be modeled using an effective field theory (EFT) [67] that includes higher-dimensional operators [68].
The event generator Vbfnlo is used to produce the Zγγ events with the aQGCs introduced using EFT dimension-8 operators with coefficients f T 0 /Λ 4 , f T 5 /Λ 4 , f T 9 /Λ 4 , f M2 /Λ 4 , and f M3 /Λ 4 in the linear Higgs-doublet representation [68] for the aQGC parameterization [69]. In this formalism, the parity conserving effective Lagrangian, which induces pure quartic couplings of the weak gauge bosons, is introduced by employing the linear representation for the higher order operators and assuming that the recently observed Higgs boson belongs to a SU(2) L doublet [68]. Dimension-8 operators are the lowest-dimension operators that lead to quartic gauge-boson couplings without exhibiting triple gauge-boson vertices. The f T,x operators contain only the field strength tensor while the f M,x operators contain both the Higgs double derivatives and the field strength. A weak boson field is either from the covariant derivative of the Higgs doublet field or from the field strength tensor. In the SM, all these aQGC operator coefficients are equal to zero. The parameters f T 0 /Λ 4 and f T 5 /Λ 4 are most sensitive to production of aQGC effects, f T 9 can only be probed via neutral QGCs such as Zγγ while f M2 /Λ 4 and f M3 /Λ 4 are chosen since they can be related to dimension-6 operators constrained by LEP experiments and CMS [69], which allows further comparisons and future aQGC combinations across different experiments. The corresponding coefficients a 0 and a c in the LEP formalism can be translated in the context of EFT dimension-8 operators (for ZZγγ/Zγγγ vertices) according to the formalism transformation equation as follows [69]: Form factors are introduced to restore unitarity at a very high parton center-of-mass energy The parameter Λ FF is chosen to preserve unitarity up to √ŝ = 8 TeV with the FF exponent n set to two.
In order to have better sensitivities to aQGCs, the measured Zγγ exclusive (zero-jet) fiducial cross section is used with the additional requirement m γγ > 300 (200) GeV for ννγγ ( + − γγ) channel. The SM backgrounds in these aQGC-optimized regions are estimated using the same methods as described in Section 5 for the Zγγ cross-section measurements. Theory predictions for the SM signal and data observations in these aQGC extended fiducial regions are shown in Table 12 Table 12: Theoretical Vbfnlo SM and observed cross sections in chosen aQGC regions (with the exclusive selection) for the channels studied. The m γγ threshold is 200 GeV for the electron and muon channels and is 300 GeV for the neutrino channel. The first uncertainty is statistical, the second is systematic.
The reconstruction efficiency C Zγγ is calculated from simulation samples with nonzero aQGCs using the events generated at LO by Vbfnlo and parton-showered by Pythia8. The deviation of the reconstruction efficiency from that for SM production using Sherpa is taken as an additional uncertainty of 20% (60%) for the νν channel ( + − channels). The differences in A Zγγ and C * (parton → particle) between aQGC and SM simulation samples are at the percent level and were neglected. The expected and observed 95% C.L. limits of each dimension-8 operator coefficient are derived from one-dimensional profile-likelihood fits as described in the aTGC study. The Λ FF -dependent observed/expected limits are obtained using the signal cross-section parameterization produced at LO by Vbfnlo and shown in Figure 16. The unitarity bounds versus Λ FF are also plotted in the figure with the FF exponent n equal to two. Table 13 shows the expected and observed 95% C.L. limits with no unitarization restriction along with those respecting unitarity bounds at the maximum allowed value of Λ FF according to the Vbfnlo estimation. The limits without unitarization are compared to the limits from the most recent CMS results [70][71][72] and ATLAS results [51] in Figure 17. The limits are presented in the formalism as implemented in Vbfnlo [73], except for the ones in Figure 17, which are presented in the formalism as implemented in MadGraph5_aMC@NLO [73] (left plot) and in the LEP formalism [69] (right plot) in order to be compared to other results.     Figure 17: Comparison of the observed limits for f T 0 /Λ 4 , f T 5 /Λ 4 , f T 9 /Λ 4 (on the left) and LEP parameters [74] a 0 /Λ 2 and a c /Λ 2 (on the right) without FF unitarization. The limits are presented in the formalism as implemented in MadGraph5_aMC@NLO [73] (left) and in the LEP formalism [69] (right).

Summary
The production cross section of Z bosons in association with isolated high-energy photons is measured using 20.3 fb −1 of pp collisions at √ s = 8 TeV collected with the ATLAS detector at the LHC. The analyses use the decays Z→νν and Z/γ * → e + e − or µ + µ − with m + − > 40 GeV. The Z/γ * decays to charged leptons are triggered on using electrons or muons with large transverse momentum. The production channels studied are pp → + − γ + X and pp → + − γγ + X where the photons are required to have E T > 15 GeV. The events with Z decays to neutrinos are selected using high-E T photon triggers. The production channels studied are pp → ννγ + X with photon E T > 130 GeV and pp → ννγγ + X where the photons have E T > 22 GeV. In all production channels the photons are required to be isolated and to satisfy tight identification criteria. The dominant backgrounds arise from jets faking photons and these are evaluated using data-driven techniques.
The cross sections and kinematic distributions for channels with Z/γ * decays to electrons and muons are combined assuming lepton universality and presented for a single charged-lepton flavor in fiducial regions defined by the lepton and photon acceptance. For the channels with Z decays to neutrinos, the cross sections and kinematics are quoted for the sum of the three neutrino flavors. This leads to studies of the following four production channels: • pp → + − γ + X, • pp → + − γγ + X, • pp → ννγ + X, • pp → ννγγ + X.
The cross sections are measured in a fiducial region, for both the inclusive case, with no requirements on the recoil system X, and the exclusive case in which there are no jets with p T > 30 GeV within |η| < 4.5.
The data are compared to SM predictions using a parton shower Monte Carlo (Sherpa) and parton-level perturbative calculations carried out at NLO (MCFM) and NNLO, corrected by parton-to-particle scale factors.
There is good agreement between the measurements and the SM predictions. Sherpa reproduces the kinematic spectra, including the jet multiplicity spectrum, in the single-photon production channels + − γ+X and ννγ+ X. The NLO and NNLO matrix element generators are used to predict the photon E T and m + − γ differential spectra in these single-photon channels, and the magnitude of the cross sections. There is good agreement between data and the SM predictions, with the NNLO calculations needed to account for the production of the high-E T photons where the NLO calculation significantly underestimates the data. In the two-photon production channels + − γγ + X and ννγγ + X the cross sections are compared to the NLO predictions. The measurements in these channels are statistically limited, but the data and SM predictions agree within the uncertainties.
Having found no significant deviations from SM predictions, the data are used to set limits on anomalous couplings of photons and Z bosons. These could result from Z/γ * s-channel production coupled to a final-state Z boson and one photon (anomalous triple gauge-boson couplings, or aTGCs), or a finalstate Z boson and two photons (anomalous quartic gauge-boson couplings, or aQGCs). The limits on the aTGCs are determined using a modified SM Lagrangian with operators proportional to parameters conventionally denoted as h V 3 and h V 4 (V = Z or γ). The contributions from aQGCs are introduced using an effective field theory concentrating on those operators most sensitive to the Zγγ final state. Limits are derived for these aTGC and aQGC parameters.