Evidence for electroweak production of two jets in association with a Zγ pair in pp collisions at s=13 TeV with the ATLAS detector

Evidence for electroweak production of two jets in association with a Z γ pair in √ s = 13 TeV proton–proton collisions at the Large Hadron Collider is presented. The analysis uses data collected by the ATLAS detector in 2015 and 2016 that corresponds to an integrated luminosity of 36 . 1 fb − 1 . Events that contain a Z boson candidate decaying leptonically into either e + e − or µ + µ − , a photon, and two jets are selected. The electroweak component is measured with observed and expected signiﬁcances of 4 . 1 standard deviations. The ﬁducial cross-section for electroweak production is measured to be σ Z γ jj − EW = 7 . 8 ± 2 . 0 fb, in good agreement with the Standard Model prediction.


Introduction
At the Large Hadron Collider (LHC), the electroweak (EW) sector of the Standard Model (SM) can be probed by studying the self-coupling of vector bosons, which is precisely predicted through the SU(2) L × U(1) Y gauge symmetry. Vector-boson scattering (VBS), VV → VV with V = W/Z/γ, is one of the most interesting processes to look at to study such effects, as it is sensitive to both the triple and the quartic gauge boson couplings, where physics beyond the SM can significantly alter the SM predictions [1-3].
In hadron collisions, VBS events are characterised by the presence of two bosons and two jets, VV j j, that are created in a purely electroweak process [4], through the interaction of bosons that have been radiated from the initial-state quarks. For such events, the scattered quarks are not colour-connected and little hadronic activity is expected in the gap between the two jets. These events are also characterised by a large invariant mass of the dijet system and a large separation of the two jets in rapidity, while the decay products of the bosons are typically produced in the central region.
The VV j j final states, however, are produced mainly through a combination of strong and electroweak interactions in pp collisions. The former are of order α 2 S α 2 EW at the Born level, where α S is the strong coupling constant and α EW is the electroweak coupling constant. In the following, such processes are referred to as QCD-induced backgrounds, Zγ j j−QCD. The production of VV j j events through pure α 4 EW interactions at the Born level are more rare; they are referred to as electroweak processes, Zγ j j−EW, in the following. It is not possible to study VBS diagrams independently from other electroweak processes as only the ensemble is gauge invariant [5]. There is also interference between the SM electroweak and QCD-induced processes. Some example Feynman diagrams for these processes are shown in Figure 1.
Among all the possible electroweak production modes, only the W ± W ± j j and W Z j j production processes have been observed [6][7][8][9]. Evidence for the Zγ j j channel has been reported [10], while limits on electroweak cross-sections have been reported for the Z Z j j channel [11], the W γ j j channel [12], and the VV channel [13,14] where at most one of the two bosons decay to two jets.
The Zγ j j channel [15] is particularly interesting to study since it allows the measurement of the neutral quartic gauge couplings, as for the Z Z final state but with a larger expected cross-section. The Zγ j j cross-sections have been computed at next-to-leading order (NLO) in α S for both the QCD-induced backgrounds [16] and the electroweak signal [17]. These corrections are not included in the following for the electroweak signal, as they were not evaluated for the phase space used in this paper.
This Letter reports the evidence for, and the most precise measurements to date of, electroweak Zγ j j production, where the Z boson decays leptonically into either e + e − or µ + µ − , exploiting the data collected with the ATLAS detector in 2015 and 2016 at a centre-of-mass energy of √ s = 13 TeV, and corresponding to an integrated luminosity of 36.1 fb −1 .

ATLAS detector
The ATLAS detector [18][19][20] is a multipurpose detector with a cylindrical geometry1 and nearly 4π coverage in solid angle.
The collision point is surrounded by inner tracking detectors, collectively referred to as the inner detector (ID), located within a superconducting solenoid providing a 2 T axial magnetic field, followed by a calorimeter system and a muon spectrometer (MS).
The inner detector provides precise measurements of charged-particle tracks in the pseudorapidity range |η| < 2.5. It consists of three subdetectors arranged in a coaxial geometry around the beam axis: a silicon pixel detector, a silicon microstrip detector and a transition radiation tracker.
The electromagnetic (EM) calorimeter covers the region |η| < 3.2 and is based on high-granularity, lead/liquid-argon (LAr) sampling technology. The hadronic calorimeter uses a steel/scintillator-tile detector in the region |η| < 1.7 and a copper/LAr detector in the region 1.5 < |η| < 3.2. The most forward region of the detector, 3.1 < |η| < 4.9, is equipped with a forward calorimeter, measuring electromagnetic and hadronic energies in copper/LAr and tungsten/LAr modules, respectively.
The muon spectrometer comprises separate trigger and high-precision tracking chambers to measure the deflection of muons in a magnetic field generated by three large superconducting toroidal magnets arranged with an eightfold azimuthal coil symmetry around the calorimeters. The high-precision chambers cover the range |η| < 2.7. The muon trigger system covers the range |η| < 2.4 with resistive-plate chambers in the barrel and thin-gap chambers in the endcap regions.
A two-level trigger system [21] is used to select events in real time. It consists of a hardware-based first-level (L1) trigger and a software-based high-level trigger. The latter employs algorithms similar to those used offline to identify electrons, muons, photons and jets.

Signal and background simulation
Signal and background processes are modelled using Monte Carlo (MC) simulations. All MC events were passed through the G 4-based [22] ATLAS detector simulation [23]. They are reconstructed and analysed using the same algorithm chain as the data. These events also include additional proton-proton interactions (pile-up) generated with P 8.186 [24] using the MSTW2008LO [25] parton distribution functions (PDFs) and the A2 set of tuned parameters [26] (A2 tune). The simulated events are reweighted to reproduce the mean number of interactions per bunch crossing observed during 2015 and 2016 data-taking. The average number of pp interactions per bunch crossing is found to be around 13 in 2015 and 25 in 2016. Scale factors are applied to simulated events to correct for the differences observed between data and MC simulation in the trigger, reconstruction, identification, isolation and impact parameter efficiencies of photons, electrons and muons [27,28]. The photon energy, electron energy and muon momentum in simulated events are smeared to account for differences in resolution between data and MC simulation [28,29].
In this analysis, the signal is the electroweak production of Z(→ )γ j j (with = e, µ), with jets originating from the fragmentation of partons arising from electroweak vertices. Signal events ( γ j j) were generated [30] at leading-order (LO) accuracy (at order α 5 EW ) using M G 5_ MC@NLO 2.3.3 [31] with no extra parton in the final state. The nominal renormalisation and factorisation scales were set to the transverse mass of the diboson system. The NNPDF30 LO PDF set [32] was used for the generation of the events, and the hadronisation and parton shower of the events was modelled using P 8.212 with the A14 tune [33].
An alternative sample of simulated events at LO accuracy was generated using S 2.2.4 [34] with Comix [35] and merged with the S parton shower [36]. The factorisation and renormalisation scales were set to the invariant mass of the diboson system. The NNPDF30 next-to-next-to-leading order (NNLO) PDF set [32] was used for the generation of this sample.
The main background comes from the QCD-induced production of the Z(→ )γ j j final state, which is estimated at leading order, i.e. at order α 2 S α 3 EW . The MC sample used to model this background was obtained with the S 2.2.2 event generator for the γ process with up to one parton at NLO generated with OpenLoops [37] and up to three partons at LO using Comix. The different final-state multiplicities were merged with the S parton shower using the ME+PS@NLO prescription [38]. This sample was generated using the NNPDF30 NNLO PDF set.
An alternative QCD sample was used to study the dependence of the kinematic distributions on the matrix-element generator. This sample was obtained with M G 5_ MC@NLO 2.3.3 with up to one parton at NLO for the γ process. The second jet is obtained from the real emission at matrix-element level, i.e. with LO accuracy. The event generator was interfaced to P 8.212 for the inclusion of the parton showers and hadronisation of the events; the FxFx merging prescription was used [39] with a merging scale of 20 GeV. The NNPDF30 NLO PDF set and the A14 tune were used for the generation.
Interference between the electroweak and QCD processes was estimated at LO accuracy in QCD using the M G 5_ MC@NLO 2.3.3 MC event generator with the NNPDF30 LO PDF set including only contributions to the squared matrix-element at order α S α 4 EW . The event generator was interfaced to P 8.212 for the parton showers and hadronisation of the events with the A14 tune. These interference effects are found to be positive and are of the order of 3% of the Zγ j j−EW cross-section in the fiducial phase space studied by this analysis.
The second-largest background in this analysis is due to Z+jets processes. In this case, one of the jets is misidentified as a photon. This contribution is estimated using a data-driven method, as explained in Section 5. Cross-checks are performed using a S 2.2.1 MC sample [40], produced with up to two jets at NLO accuracy and up to four jets at LO accuracy, using Comix and OpenLoops, and merged with the S parton shower according to the ME+PS@NLO prescription. The NNPDF30 NNLO PDF set was used for the generation of this sample. This sample is scaled to the NNLO predictions for inclusive Z production using a normalisation factor [40] derived from FEWZ [41] and the MSTW2008NNLO PDF set [25].
The third-largest background affecting this analysis is from tt + γ events. These events were generated at LO accuracy using the M G 5_ MC@NLO 2.3.3 generator interfaced to P 8.212 for the parton showering and hadronisation. The NNPDF23 LO PDF set and the A14 tune were used for the generation. The normalisation of the predictions is corrected for NLO QCD effects [42].
All other backgrounds are smaller and are obtained from MC simulations. The most important of these are the W Z background, which was simulated with S 2.2.1, and the Wt background, which was simulated with the P -Box [43][44][45] generator interfaced to P 6.428 [46] and the 2012 tune [47]. The CT10 NLO PDF set [48] was used for the generation of this sample.

Event selection
Events are selected if they were recorded during stable beam conditions and if they satisfy detector and data quality requirements, which include all relevant subdetectors functioning normally. Events were collected using single-lepton or dilepton triggers [21,49] that require, respectively, at least one or two electrons or muons.
For single-lepton triggers, the transverse momentum (p T ) threshold in 2015 was 24 GeV for electrons and 20 GeV for muons satisfying a loose isolation requirement based only on ID track information. In 2016, due to the higher instantaneous luminosity, these thresholds were increased to 26 GeV for both the electrons and muons, and tighter isolation requirements were applied. Inefficiencies for leptons with large transverse momenta were reduced by including additional electron and muon triggers that do not include any isolation requirements; these had transverse momentum thresholds of 60 GeV and 50 GeV, respectively. Finally, a single-electron trigger requiring p T > 120 GeV with less restrictive electron identification criteria was used to increase the selection efficiency for high-p T electrons. The threshold was increased to 140 GeV in 2016. In 2015, the dilepton trigger p T threshold was 12 GeV for electrons, and either 10 GeV for muons in conjunction with two muons producing L1 triggers, or 18 GeV and 8 GeV in the case that only one of the two muons produced an L1 trigger. In 2016, these thresholds were raised to 17 GeV for electrons, 14 GeV for the symmetric dimuon trigger, and 22 GeV and 8 GeV for the asymmetric muon trigger. The combined efficiency of these triggers is close to 100% for the events that pass the kinematic requirements described below. To ensure that the trigger efficiency is well determined, the lepton candidates must be matched to the leptons that are selected by the trigger and have a transverse momentum at least 1 GeV above the online threshold.
Events must have a primary vertex with at least two charged-particle tracks which must be compatible with the pp interaction region. In cases where multiple vertices are reconstructed in a single event, the vertex with the highest sum of the p 2 T of the associated tracks is selected as the production vertex of the Zγ system.
Muons are reconstructed by matching tracks in the muon spectrometer to tracks in the inner detector. Candidate muons are required to satisfy the 'medium' set of identification criteria [28] based on the number of detector hits and p T measurements in the ID and the MS. The efficiency of selecting such muons averaged over p T and η is larger than 98%. Both the ID and MS measurements are used to compute the muon momentum. The measurement also takes into account the energy loss in the calorimeters. Muons are used in the analysis if they satisfy p T > 20 GeV and |η| < 2.5.
Electrons are reconstructed by matching energy clusters in the electromagnetic calorimeter to an ID track. Candidate electrons are required to pass a 'medium' likelihood requirement that is built from information about the electromagnetic shower shape in the calorimeter, track properties and the track-to-cluster matching [27]. The efficiency of selecting such electrons is of about 90%. Calorimeter energy and the track's direction are combined to compute the electron momentum. Electrons are required to satisfy p T > 20 GeV and |η| < 2.47.
To ensure that electron and muon candidates originate from the primary vertex, the significance of the track's transverse impact parameter relative to the beam line must satisfy |d 0 /σ d 0 | < 3 (5) for muons (electrons), and the longitudinal impact parameter z 0 , the difference between the value of z at the point on the track at which d 0 is defined and the longitudinal position of the primary vertex, is required to satisfy |z 0 · sin(θ)| < 0.5 mm.
Electrons and muons are required to be isolated from hadronic activity by imposing maximum energy requirements in cones in the η-φ plane, defined using the ∆R = (∆η) 2 + (∆φ) 2 distance, around their direction of flight, using calorimeter-cluster and ID-track information, excluding the electron or muon momentum. The muon isolation criterion imposes a fixed upper limit in a cone of size ∆R = 0.2 for calorimeter-based isolation and in a cone of size ∆R = 0.3 for ID-based isolation, with a selection efficiency above 90% for p T > 20 GeV and at least 99% for p T > 60 GeV [28]. For electrons, the requirements vary with p T , allowing a selection efficiency of at least 90% for p T > 25 GeV and at least 99% for p T > 60 GeV [27]. For electrons, both isolation variables are evaluated with a cone of size ∆R = 0.2.
Photon candidates are reconstructed from clusters of energy deposited in the electromagnetic calorimeter and are classified as unconverted (clusters without a matching track or matching reconstructed conversion vertex in the ID) or converted (clusters with a matching reconstructed conversion vertex or a matching track consistent with originating from a photon conversion). Photons are identified using a cut-based selection relying on shower shapes measured with the electromagnetic calorimeter. Nine discriminating variables are used, built from the distribution of energy in different layers of the EM calorimeter and information about energy leaking into the hadronic calorimeter. The shower shape values in the simulation are corrected to improve their agreement with the shower shapes in data. Photons must satisfy a 'tight' requirement [50]. Photons selected in this analysis must satisfy E T > 15 GeV, where E T is the photon transverse energy, and the pseudorapidity of the cluster must be in the range |η| < 1.37 or 1.52 < |η| < 2.37, to avoid the transition region between barrel and endcap calorimeters. Photon candidates are required to be isolated according to calorimeter-cluster and ID-track information in a cone of size ∆R = 0.2, excluding the photon energy. The requirement varies with E T following the 'fixed loose' requirement [50], providing an average efficiency of about 95%. The calorimetric photon isolation is corrected to account for leakage and ambient transverse energy [27].
Jets are reconstructed from clusters of energy depositions in the calorimeter [51] using the anti-k t algorithm [52,53] with a radius parameter R = 0.4. Calibration of the jet energy is based on simulation and in situ methods from data [54]. Jets originating from non-collision backgrounds or detector noise are removed [55]. Pile-up jets are suppressed using a multivariate combination of track-based variables in the ID acceptance (|η| < 2.5) [56]. All jets considered must have p T > 25 GeV and must be reconstructed in the pseudorapidity range |η| < 4.5.
Jets containing a b-hadron are identified in the ID volume using a multivariate algorithm [57,58] that uses the impact parameters and secondary vertices of the tracks contained in the jet. The selection is done using a working point chosen to provide a 70% selection efficiency for b-jets in an inclusive tt MC sample. The 70% working point has rejection factors of 12 and 381 for charm and light-flavour jets, respectively. Correction factors are applied to the simulated event samples to compensate for differences between data and simulation in the flavour-tagging efficiencies for b-jets, c-jets and light-flavour jets. The correction for b-jets is derived from tt events with final states containing two leptons, and the corrections are consistent with unity with uncertainties at the level of a few percent over most of the jet p T range.
To avoid cases in which a lepton, photon, or jet is reconstructed as two separate final-state objects, several steps are followed to remove such overlaps. Bremsstrahlung radiation by a muon can result in ID tracks and a calorimeter energy deposit that are reconstructed as an electron candidate. Therefore, in cases for which an electron candidate and a muon candidate share an ID track, the object is considered to be a muon, and the electron candidate is rejected. The overlap of objects is measured using the ∆R distance. Due to the isolation requirements placed on electron candidates, any jets that closely overlap an electron candidate within a conical region ∆R < 0.2 are likely to be reconstructions of the electron and so are rejected. When jets and electrons are found within a larger hollow conical region 0.2 < ∆R < 0.4, it is more likely that a real hadronic jet is present and that the electron is a non-prompt constituent of the jet, arising from the decay of a heavy-flavour hadron. Hence an electron candidate found within ∆R < 0.4 of any remaining jet is rejected. Muons can be accompanied by a hard photon due to bremsstrahlung or collinear final-state radiation, and the muon-photon system can then be reconstructed as both a jet and a muon candidate. Non-prompt muons can arise from decays of hadrons in jets; however, these muons are associated with a higher ID track multiplicity than those accompanied by hard photons. In order to resolve these ambiguities between nearby jet and muon candidates, first any jets having fewer than three ID tracks and within ∆R < 0.4 of any muon candidate are rejected, then any muon candidates within ∆R < 0.4 of any remaining jet are rejected. The ambiguity between photons and leptons is resolved by rejecting any photons found to be within ∆R < 0.4 of an electron or a muon. Jets found within ∆R < 0.4 of any photons are rejected.
Events are required to contain exactly two lepton candidates (electrons or muons) with the same flavour and opposite charge and at least one photon candidate satisfying the selection criteria described above. When there are two or more photons selected, the one with the highest E T is used for further processing.
The invariant mass of the two leptons, m , must be at least 40 GeV. The sum of the dilepton mass and the three-body γ invariant mass (m + m γ ) is required to be larger than 182 GeV, which is twice the Z boson mass. This requirement ensures that the three-body invariant mass is larger than the Z boson mass, thus suppressing cases where the Z boson decays into γ [15].
To select electroweak processes, events are further required to include at least two 'VBS-tagging' jets. A symmetric selection of p T > 50 GeV is used to select both tagging jets, which are classified as leading and subleading in transverse momentum. These two jets are required to have an invariant mass m j j > 150 GeV, in order to minimise the contamination from triboson background in which a hadronically decaying W or Z boson is produced in association with a Z boson and a photon, i.e. Zγ + Z/W(→ j j).
The final VBS signal region (SR) is defined by requiring that the pseudorapidity difference between the two leading jets |∆η( j j)| be larger than 1; that no b-tagged jet be present in the event; and that the centrality of the γ system relative to the tagging jets, defined as be lower than 5, where y corresponds to the rapidity, and j 1 and j 2 label the selected leading and subleading jets.

Background estimation
Irreducible backgrounds are modelled kinematically by MC simulations, and their normalisation is evaluated from the data. Reducible backgrounds, in which a photon originates from a hadron decay, are determined using data-driven techniques. Events in which one or more leptons originate from a hadronic jet are suppressed due to the stringent kinematic requirements applied. They are found to be negligible and are not considered here.
The main reducible background originates from Z+jets events, in which one hadronic jet is misidentified as a photon. This background is not well modelled by the MC simulation and, therefore, is estimated with data using a two-dimensional sideband method [59] similar to that applied in a previous analysis with Run 1 data [15]. The shapes of the distributions of Z+jets events are obtained using data, after applying all the signal requirements, but reversing the photon identification and isolation criteria. The normalisation of this background is estimated in a control region designed to maximise the number of events available in the data while being close to the signal region; all signal region requirements are applied, except the transverse momentum requirement for the jets, which is lowered to 30 GeV, and the dijet invariant mass requirement, which is less than 150 GeV. The ratio of Z+jets to Zγ j j−QCD events is computed in this region, since it is found to be the same in this control region and in the signal region, in both data and MC simulation, and used to extrapolate the normalisation of Z+jets events to the signal region. To estimate both the shape and the normalisation of this background, the contamination from events containing real photons is estimated and corrected using MC predictions.
The most important background is irreducible and comes from the QCD-induced production of Zγ events in association with two jets. The normalisation of this background is constrained by the data, as explained in Section 6.
The other main irreducible background arises from tt + γ events. This process is also constrained by the data in a dedicated control region which is referred to as b-Control Region (b-CR). This control region is obtained by imposing all SR selections but requiring the presence of at least one reconstructed b-tagged jet in the event.
The contribution of backgrounds due to the production of Z+jets and γ+jets events in two different overlaid interactions is evaluated using a data-driven estimate and found to be negligible.
Backgrounds from other sources are reducible and are evaluated using MC simulations. The most important contribution comes from W Z events, in which a charged lepton is misidentified as a photon, and from Wt events.

Signal extraction procedure
A boosted decision tree (BDT), as implemented in the TMVA package [60], is used to separate the electroweak signal from all the backgrounds. The BDT is trained and the set of variables retained as input to the BDT is optimised on simulated events to separate Zγ j j−EW events from all background processes, excluding Z+jets, and to maximise the signal-to-background ratio. In total, 13 kinematic variables are combined into one discriminant that can take any value in the range [−1, 1]. The final binning of the BDT discriminant, or BDT score, is optimised to improve the sensitivity of the analysis.
Some of the variables used in the BDT training are related to the kinematic properties of the two tagging jets: the invariant mass of these two jets, m j j , the difference in pseudorapidity between these two jets, ∆η j j , and the pseudorapidity and p T of the leading jet. Other variables are related to the kinematic properties of the photon and the Z boson: the transverse momenta of the leading lepton, the system and the γ system, and the and γ invariant masses. Another set of variables is related to the correlation of the two tagging jets and the Zγ system: the distance in the η-φ plane between the Zγ system and the two-jet system, ∆R(Zγ, j j); the smallest distance in the η-φ plane between a photon and a jet, among all possible combinations, Min(∆R(γ, j)); the difference in azimuth between the Zγ system and the two-jet system, ∆φ(Zγ, j j) and the centrality ζ( γ) as defined in Eq. (1).
The modelling by MC simulations of the distribution shapes and the correlations between all input variables for the BDT is verified in the signal region. Compatibility within statistical uncertainty is observed for all distributions, as exemplified by Figures 2(a) and 2(b), except for the high-mass tail of m j j as shown in Figure 2(c), as was also observed in previous analyses of electroweak processes [61-63]. This mismodelling does not impact the signal significance nor the uncertainty in the measurement reported in Section 9. The BDT score distribution in the signal region is used to measure the Zγ j j−EW signal fiducial crosssection with a maximum-likelihood fit, and to determine the significance of the signal. In this fit the eeγ j j and µµγ j j final states are combined. An extended likelihood function is built from the product of two likelihoods corresponding to the BDT score distribution in the Zγ j j SR, and the multiplicity of reconstructed b-jets in the b-CR. The yields of the tt + γ simulated events are constrained by data using this last distribution. The normalisation of the Zγ j j−QCD sample is also constrained by the data, directly in the signal region, since this background dominates the distribution at low BDT score values. The normalisations of these backgrounds are introduced in the likelihood fit as unconstrained nuisance parameters, labelled as µ Zγ j j−QCD and µ tt+γ for the Zγ j j−QCD and tt + γ backgrounds, respectively. For these backgrounds, the kinematic distributions are obtained from MC simulations. Both the shape and the normalisation distributions of Z+jets background are estimated from data, as described in Section 5.
The other irreducible backgrounds are determined from MC simulations. Background normalisations and shapes can vary within the uncertainties, constrained by Gaussian distributions as described in Section 7.
The fiducial cross-section is derived using the signal strength parameter µ Zγ j j−EW : where N data Zγ j j−EW is the number of signal events measured in the data, while N MC Zγ j j−EW is the number of signal events predicted by the M G 5_ MC@NLO 2.3.3 MC simulation.
The observed cross-section (σ fid. ,data Zγ j j−EW ) is then obtained by multiplying the signal strength µ Zγ j j−EW by the M G 5_ MC@NLO 2.3.3 MC cross-section prediction in the fiducial region (σ fid., MC Zγ j j−EW ). Because the effect of interference between the Zγ j j−QCD and the Zγ j j−EW processes is not accounted for in the Zγ j j−QCD contribution, the observed cross-section σ fid.
Zγ j j−EW formally corresponds to electroweak production plus the interference effects.

Systematic uncertainties
Several sources of systematic uncertainty in the signal and background processes can affect the cross-section measurement. For the background and signal processes, systematic uncertainties that affect the shapes of the BDT score and N b-jets distributions are considered. For the signal, systematic uncertainties that only affect the acceptance are considered, while for the background processes, those affecting the normalisation of the distributions are also included. Systematic uncertainties in the shapes of distributions are not applied if they are consistent with statistical fluctuations.
The uncertainties in the theoretical modelling of the signal and backgrounds that are estimated from MC simulations can affect the cross-section measurement. The uncertainties due to the knowledge of the PDF and the α S value used in the generation of the MC samples are determined using the PDF4LHC prescription [64]. They are found to be of the order of 2% for the normalisation of the EW, QCD and ttγ samples. The uncertainties due to missing higher-order QCD corrections are evaluated using an envelope of the largest deviations obtained by varying the renormalisation and factorisation scales independently by factors of two and one-half. These uncertainties are large for the Zγ j j−QCD and tt + γ backgrounds and their impact on the normalisation is +30 −20 %, while the impact on the Zγ j j−EW signal normalisation is smaller and of about ±5%. The shape uncertainties due to the missing higher-order QCD corrections are of 2-4% for the signal and the Zγ j j−QCD and tt + γ backgrounds.
For the signal and tt + γ background, uncertainties due to the parton shower and underlying-event modelling are evaluated using dedicated tune variations [33] and by re-showering events with Herwig 7.0.1 [65, 66] instead of P 8.212. These uncertainties affect the shape of the distribution by at most 5% at large values of the BDT score and are negligible at low values of the BDT score. For the Zγ j j−QCD background, the modelling uncertainties due to the parton showers, underlying-event and matrix elements are estimated by comparing the predictions of the S 2.2.2 and M G 5_ MC@NLO 2.3.3 MC generators. The difference between these two predictions is taken as a systematic uncertainty applied as both a positive and a negative variation. The effects in the BDT score distribution and the N b-jets distribution due to these two set-ups are taken as uncertainties. The main impact is on the BDT score distribution with an effect ranging from −5% to 20% at low and high values of the BDT score, respectively.
As discussed in Section 6, the interference between the electroweak signal and the Zγ j j−QCD process is not included in the fit. However, this interference can distort the shape of the Zγ j j−EW template. This effect is estimated using the M G 5_ MC@NLO 2.3.3 MC generator at LO in QCD, and is taken as a systematic uncertainty in the shape of the signal. The size of this effect ranges from 5% to 2% at low and high values of the BDT score, respectively.
Another source of uncertainty arises from the sizes of the MC samples and data sample in the regions used in the analysis to model the BDT score and N b-jets distributions. The statistical uncertainty in the shape of the BDT score for Zγ j j−QCD events ranges from 5% to 13% at low and high values of the BDT score, respectively, and is about 2% for Zγ j j−EW events. For the Z+jets background, the statistical uncertainty dominates the shape uncertainty due to the size of the data sample; it is about 10% at low BDT score values, and up to 50% at high values, while uncertainties from other sources are negligible.
Other systematic uncertainties originate from the reconstruction, identification, and energy calibration of electrons, photons, muons, and jets. The largest of these uncertainties is due to the jet energy scale calibration. The uncertainties due to the jet energy resolution and to the suppression of pile-up jets are also considered, as described in Ref. [67]. The total effect of the uncertainties related to jet reconstruction and calibration on the normalisation of the Zγ j j−QCD background is about 8%, and it is about 4% for Zγ j j−EW events in the signal region. The impact on the shape is largest for the Zγ j j−QCD events and ranges from 2% to 18% at low and high BDT values, respectively. The uncertainty due to the heavy-flavour tagging efficiency amounts to 9% (2%) in normalisation in the signal (b-CR) region for the tt + γ background, whereas the impact on other backgrounds and on the shape of the template distributions is found to be negligible. Uncertainties in the lepton identification, reconstruction, isolation requirements, trigger efficiencies, energy scale and energy resolution are determined by using Z → events [27, 28, 68]. The largest experimental uncertainty associated with the photons is the photon identification efficiency [50]. The photon and lepton uncertainties affect only the normalisation of Zγ j j−EW, Zγ j j−QCD, and tt + γ processes. In total, these effects are of 3% to 4% in the signal region.
A 20% yield uncertainty is assigned to the Z+jets reducible background estimate. This uncertainty accounts for the number of events in the control regions used in the two-dimensional sideband measurement, and for the correlation between photon identification and isolation requirements.
A 20% yield uncertainty is assigned to the other backgrounds estimated from simulations.
A variation in the pileup reweighting of MC samples is included to cover the uncertainty on the ratio between the predicted and measured inelastic cross-section [69]. The resulting uncertainty in the measured fiducial cross-section is 5%.
The uncertainty in the combined 2015-2016 integrated luminosity is 2.1% [70], obtained using the LUCID-2 detector [71] for the primary luminosity measurements.
The effect of each of these uncertainties on the fiducial cross-section measurement is shown in Table 1. The individual sources are grouped into either theoretical or experimental categories. The largest uncertainties are due to the jet reconstruction and calibration, followed by the theoretical uncertainties in the modelling of the Zγ j j−EW signal and of the Zγ j j−QCD background. Uncertainties affecting only the normalization of the Zγ j j−QCD and ttγ background have almost no impact on the Zγ j j−EW cross-section measurement since the corresponding normalization parameters are constrained by the fit to the data as explained in Section 6.

Phase space for cross-section measurements
The Zγ j j electroweak cross-section is measured in a fiducial phase space that is defined to closely follow the selection criteria of the signal region described in Section 4. This region is defined at the particle level, using stable particles (with a proper decay length cτ > 10 mm) which are produced from the hard scattering, including those that are the products of hadronisation, before their interaction with the detector. Leptons produced in the decay of a hadron, a τ-lepton, or their descendants are not considered in the definition of the fiducial phase space. Prompt-lepton four-momenta are obtained from a sum of the leptons four-momenta and the four-momenta of photons not originating from hadron decays within a cone of size ∆R = 0.1 around the leptons ('dressed leptons'). Jets are reconstructed using the anti-k t algorithm with radius parameter R = 0.4 using stable particles, excluding electrons, muons, neutrinos, and photons associated with the decay of a W or Z boson. Photon isolation is defined as the transverse momentum of the system of stable particles within a cone of size ∆R = 0.2 around the photon, excluding muons, neutrinos, and the photon itself.
The following selection requirements are imposed to define the fiducial phase space. The charged leptons from the Z boson decay are required to have transverse momentum p T > 20 GeV and |η| < 2.5, and the invariant mass of the two leptons must be larger than 40 GeV. The photon is required to have transverse momentum p T > 15 GeV and |η| < 2.37. Photons isolated from any hadronic activity are selected by requiring that the photon isolation, as defined above, divided by the photon transverse momentum be less than 5%. The angular distance between each of the charged leptons from the Z decay and the photon is required to be ∆R > 0.4. The sum of the two-lepton invariant mass and two-lepton plus photon invariant mass must satisfy (m + m γ )> 182 GeV in order to exclude final-state radiation events.
In addition to these requirements, at least two jets with p T > 50 GeV and |η| < 4.5 are required. The angular distance between each of the charged leptons from the Z boson decay and each of the jets is required to be ∆R( j, ) > 0.3. The angular distance between the photon and each of the jets is required to be ∆R( j, γ) > 0.4. The invariant mass m j j of the two highest-p T jets is required to be m j j > 150 GeV, and the difference between the pseudorapidities of these two leading jets is required to be |∆η j j | > 1. Finally, the centrality of the γ system relative to the tagging jets must satisfy ζ( γ)<5.

Cross-section measurements
A profile-likelihood-ratio test statistic [72] is used to measure the signal strength µ Zγ j j−EW and its associated uncertainties. The systematic uncertainties are treated as nuisance parameters and are constrained with Gaussian distributions in the fit. The statistical uncertainty is determined after having fixed each nuisance parameter to its profiled value. Figure 3 shows the BDT score distribution in the signal region, and the b-tagged jet multiplicity in the b-CR, after having applied the background and signal normalisations and having set the nuisance parameters to the values adjusted by the profile-likelihood fit. The yield of events obtained after the fit is detailed in Table 2.
The signal strength is measured to be:   Table 2: Expected and observed numbers of events in the signal and control regions, after the fit. The expected number of Zγ j j−EW events from M G 5_ MC@NLO 2.3.3 and the estimated number of background events from the other processes are shown. Total post-fit uncertainties, as described in Section 7, are shown. The uncertainty in the total expected yield is smaller than the quadratic sum of the uncertainties in each process due to these being anti-correlated after the fit.

SR
b-CR The fiducial cross-section is measured by computing the product of signal strength and the predicted cross-section used to define it: Zγ j j−EW = 7.8 ± 1.5 (stat.) ± 1.0 (syst.) +1.0 −0.8 (mod.) fb = 7.8 ± 2.0 fb. This cross-section corresponds to the electroweak Zγ j j particle-level cross-section in the fiducial phase space defined in Section 8 using dressed leptons, and including constructive interference between the signal and the QCD-induced processes. To verify the results of the multivariate analysis, a cut-based approach is also used. The centrality of the γ system, ζ( γ), is used as the sensitive variable for the extraction of the electroweak signal. In this alternative method, the signal region is further divided into two regions depending on the dijet invariant mass. The region with m j j < 500 GeV is used to constrain the QCD background, while the region m j j > 500 GeV is kept to extract the signal. The b-CR is also kept in this method as a second control region to constrain the tt + γ background normalisation. The profile-likelihood fit is applied in the same way as described above with the full treatment of systematic uncertainties, allowing exclusion of the background-only hypothesis with a statistical significance of 2.9σ (2.7σ) observed (expected). This is compatible with the result obtained with the multivariate approach and with the SM prediction.
The analysis is also used to measure the Zγ j j cross-section in the same fiducial phase space. This measurement includes the electroweak-induced signal, the QCD-induced process and their interference. The extraction procedure is the same as for the electroweak Zγ j j measurement. A template likelihood fit is performed, using the BDT score distribution in the signal region. For this measurement, the b-CR is not used, to simplify the statistical model. The omission of this control region in the fit does not change the result. The signal template corresponds to the sum of the Zγ j j−EW and Zγ j j−QCD templates. The cross-section is measured to be:

Conclusion
Evidence for electroweak production of two jets in association with a Zγ pair is presented using 36.1 fb −1 of pp collision data at √ s = 13 TeV collected with the ATLAS detector at the LHC. The production cross-section of this process is measured in a fiducial phase space approximating the acceptance of the analysis. This measurement uses the leptonic decay of the Z boson into e + e − or µ + µ − . The measurement is performed using a BDT to enhance the signal to background-ratio. The dominant backgrounds, Zγ j j−QCD, Z+jets and tt + γ are all estimated from the data. The background-only hypothesis is excluded with observed and expected significances of 4.1 standard deviations.
In the fiducial phase space the electroweak cross-section of the Zγ j j process is measured to be: Zγ j j−EW = 7.8 ± 1.5 (stat.) ± 1.0 (syst.) +1.0 −0.8 (mod.) fb, = 7.8 ± 2.0 fb in good agreement with the Standard Model predictions at LO in perturbative QCD.