Measurement of 𝒁𝜸𝜸 production in 𝒑 𝒑 collisions at √ 𝒔 = 13 TeV with the ATLAS detector

Cross-sections for the production of a 𝑍 boson in association with two photons are measured in proton–proton collisions at a centre-of-mass energy of 13 TeV. The data used correspond to an integrated luminosity of 139 fb − 1 recorded by the ATLAS experiment during Run 2 of the LHC. The measurements use the electron and muon decay channels of the 𝑍 boson, and a ﬁducial phase-space region where the photons are not radiated from the leptons. The integrated 𝑍 (→ ℓℓ ) 𝛾𝛾 cross-section is measured with a precision of 12% and diﬀerential cross-sections are measured as a function of six kinematic variables of the 𝑍𝛾𝛾 system. The data are compared with predictions from MC event generators which are accurate to up to next-to-leading order in QCD. The cross-section measurements are used to set limits on the coupling strengths of dimension-8 operators in the framework of an eﬀective ﬁeld theory.


Introduction
Processes involving the production of three electroweak (EW) gauge bosons from proton-proton collisions are typically rare. Some of these triboson processes are only just becoming accessible due to the unprecedented integrated luminosity provided during Run 2 of the LHC. Measurements of such processes provide a direct probe of non-Abelian quartic gauge couplings, both those that are predicted by the Standard Model (SM) and those that could only be due to new physics. The production of a boson in association with two prompt photons provides an opportunity to test the electroweak sector of the SM and to constrain any potential new physics effects. Neutral quartic gauge couplings are not allowed in the SM, and hence the production of has no tree-level contribution involving quartic couplings. In this paper, the leptonic decay channels of the boson, i.e. → ℓ + ℓ − where ℓ = , , are considered. Despite having lower branching fractions than the quark and neutrino decay channels, the leptonic channels benefit from having a cleaner final state and smaller backgrounds. The measurement of is also crucial for our understanding of the irreducible background to (→ ℓℓ) (→ ) production, and for searches for resonances in the ℓℓ final state.
The production of ℓℓ from proton-proton collisions proceeds at leading order by diagrams of the first three types given in Figure 1. The production of ℓℓ via three on-shell bosons is shown in Figure 1 to this signal arise from processes involving jets which are misidentified as photons. Previous studies of the ℓℓ final state with the ATLAS detector at √ = 8 TeV [1], and the CMS detector at 8 TeV [2] and 13 TeV [3], were performed in phase spaces which included both the ISR and FSR production of photons. The measurement presented in this paper suppresses the FSR contribution, which allows a simpler interpretation of the measurements since the triboson process produces the dominant signal contribution. The ℓℓ final state is also sensitive to new physics via anomalous quartic couplings, an example of which is shown in Figure 1(d).
The measurements use → + − + and → + − + events recorded by the ATLAS detector at √ = 13 TeV. The ATLAS detector is described in Section 2. The full Run 2 dataset corresponding to an integrated luminosity of 139 fb −1 is used. It is described in Section 3 along with simulated event samples. The event selection is given in Section 4. Background processes are estimated using a combination of data-driven techniques and simulation, which are described in Section 5. The systematic uncertainties are discussed in Section 6. The yields of signal events are corrected (unfolded) to a fiducial volume where the integrated (differential) cross-section is measured; the unfolding procedures and results are described in Section 7. Differential cross-sections are measured as a function of the transverse energy 1 T of the leading (highest T ) photon, the transverse energy 2 T of the subleading photon, the transverse momentum ℓℓ T of the dilepton system, the transverse momentum ℓℓ T of the four-body system, the invariant mass of the diphoton system, and the invariant mass ℓℓ of the four-body system. The distribution is a measure of the hadrons recoiling against the ℓℓ system and is hence sensitive to QCD modelling. The distribution is useful for constraining backgrounds to resonances in the ℓℓ final state, and the ℓℓ spectrum describes the scale of the full four-body system. The data are compared with predictions from Monte Carlo (MC) event generators with matrix elements calculated to up to next-to-leading order (NLO) in perturbative QCD. The differential measurements are also used to constrain new physics effects arising through anomalous neutral quartic couplings. This is done via an effective field theory (EFT) approach, where limits are set on the coupling strengths of dimension-8 operators. This procedure and the measured limits are presented in Section 8. The conclusions are stated in Section 9.

The ATLAS detector
The ATLAS experiment [4] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4 coverage in solid angle. 1 It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range | | < 2.5. It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (referred to as the barrel), covering | | < 1.7. The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to | | = 4.9. The muon spectrometer surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The muon spectrometer (MS) includes a system of precision tracking chambers and fast detectors for triggering. A two-level trigger system is used to select events. The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz. This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions. An extensive software suite [5] is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Data and simulated event samples
The data used in this analysis were collected in proton-proton collisions at √ = 13 TeV from 2015 to 2018. After applying criteria to ensure normal ATLAS detector operation [6], the total integrated luminosity useful for data analysis is 139 fb −1 . The uncertainty in the total Run 2 integrated luminosity is 1.7% [7], obtained using the LUCID-2 detector [8] for the primary luminosity measurements. The average number of inelastic interactions produced per bunch-crossing for the dataset considered is = 33.7. 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of Simulated event samples are used to correct the background-subtracted data yield for detector effects and to estimate several background contributions. The simulated samples were produced with various MC event generators, processed through a full ATLAS detector simulation [9] based on G 4 [10], and reconstructed using the same algorithms as used for data. All simulated samples were corrected with data-driven correction factors to account for differences in trigger, reconstruction, identification and isolation performance between data and simulation. Additional interactions (pile-up) occurring in the same or neighbouring bunch-crossings were modelled by overlaying each simulated event with minimum-bias events generated using P 8.186 [11] with the A3 set of tuned parameters [12] and the NNPDF2.3 [13] set of parton distribution functions (PDFs). The simulated events were then reweighted to reproduce the distribution of the number of interactions per bunch-crossing observed in the data.
Samples of simulated + − and + − events were generated using S 2.2.10 [14] at NLO accuracy in QCD for zero additional partons and LO accuracy for up to two additional partons. Matrix elements were matched and merged with the S parton shower [15] based on Catani-Seymour dipole factorisation [15,16] using the MEPS@NLO prescription [17][18][19][20]. The NNPDF3.0 [21] set of PDFs were used.
For studies of systematic uncertainties and cross-checks, an alternative signal sample is considered. It was produced with the M G 5_ MC@NLO 2.7.3 [22] generator with up to one additional final-state parton at NLO accuracy, using the NNPDF3.0 PDF set. Events were interfaced to P 8.244 [23], via the FxFx merging procedure [24], for modelling of the parton shower, hadronisation and underlying event. Both the baseline and alternative signal samples utilise smooth-cone photon isolation [25], with the parameters 0 = 0.1, = 0.1 and = 2, to remove contributions from fragmentation photons.
The + jets ( + jets) samples used in the estimation of misidentified-photon backgrounds were generated with S 2.2.4 (P B v1 [26][27][28][29]). The¯sample used in the estimation of the background where the leptons originate from top quark decays was generated with M G 5_ MC@NLO 2.3.3. The backgrounds arising from electrons misidentified as photons were modelled with S 2.2.2 ( → ℓℓℓℓ) and S 2.2.5 ( → ℓ ℓℓ ). The contribution from (→ ℓℓ) (→ ) was generated with P B v2. The single-photon and diphoton samples used in the estimation of the pile-up overlay backgrounds were modelled using S 2.2.2. Further details are given in Table 1, along with a summary of the signal samples.

Event reconstruction and selection
Events are selected in the electron and muon channels using unprescaled single-lepton triggers [31,32] with the lowest T threshold available. From 2016 to 2018, this was 26 GeV for both electrons and muons, and in 2015 it was 24 GeV for electrons and 20 GeV for muons. The low-T triggers are supplemented by higher ones with relaxed identification and isolation requirements which improve the overall trigger efficiency. Reconstructed tracks in the ID and clusters of energy deposits in the EM calorimeter are used as inputs in the reconstruction of electrons and photons [33]. Reconstructed track segments in the ID and MS are used as inputs in the reconstruction of muons [34].
Electron candidates are seeded by EM calorimeter energy clusters and must have T > 20 GeV, | | < 2.47, and a matching ID track. The transition region between the barrel and endcap of the EM calorimeter (1.37 < | | < 1.52) is excluded. Electrons are identified using a likelihood discriminant formed from shower shape variables, track variables and a measure of how well the track matches the cluster. All electron Table 1: Summary of simulated MC event samples for the ℓℓ signal process and those used in the estimation of backgrounds. The third and fourth columns give the order in perturbative QCD and the PDF set used in the hard-scattering matrix element calculations. The rightmost column specifies the generator used to model parton showering, hadronisation, the underlying event and multiple parton interactions.

Process
Generator candidates must satisfy the Medium identification working point [33]. To suppress the contribution from jets misidentified as electrons, the electron candidates must be isolated from other activity in the tracking and calorimeter systems. The calorimeter and tracking isolation variables are constructed, respectively, from the sums of cluster energies and track momenta falling within a cone of size Δ = 0.2 around the electron, which are then required to satisfy the Loose working point described in Ref. [33].
Muon candidates are formed from tracks and must have T > 20 GeV and | | < 2.5. Identification requirements comprise selections on track quality and a measure of how well the ID track matches the MS track. Muon candidates must pass the Medium identification working point [34]. The muon candidates are also required to be isolated in the tracking and calorimeter systems using variables similar to those for electrons. The Loose isolation working point is used and is similar to that described in Ref. [34].
Reconstructed tracks matched to common points of origin along the beam axis serve as candidates for the location of proton-proton collisions (vertices). The vertex with the largest sum of the track 2 T is chosen as the primary vertex (PV). Electrons and muons must be consistent with originating from the PV, requiring that their transverse impact parameter significance satisfies | 0 |/ 0 < 3 (5) for muons (electrons) and the longitudinal distance 0 from the PV to the point where 0 is measured satisfies | 0 sin( )| < 0.5 mm.
Photon candidates are seeded by EM calorimeter energy clusters which have T > 20 GeV and | | < 2.37. The transition region between the barrel and endcap of the EM calorimeter (1.37 < | | < 1.52) is excluded. Converted photon candidates are formed from clusters which are matched to a conversion vertex. The conversion vertices are formed from one or two tracks which are consistent with a massless particle converting within the ID volume. Unconverted photon candidates are formed from clusters which are not matched to an electron track or conversion vertex. Photon candidates are required to pass a number of selections on shower shape variables which correspond to the Loose identification working point [33].
Overlap removal requirements are applied to the preliminarily selected objects to prevent the same particle from being reconstructed as two separate physics objects. Photons are removed if they are within Δ = 0.4 of a selected electron or muon. Electrons are removed if they are within Δ = 0.2 of a selected muon.
Correction factors are applied to the selected objects to account for object trigger, reconstruction, identification and isolation efficiency differences between data and simulation.
Candidate events are considered further if they contain at least one opposite-sign same-flavour lepton pair and at least two photons. One of the leptons must be matched to the trigger object which fired the event, and the highest-T (leading) lepton must have T > 30 GeV to be well above the trigger threshold. If an electron (muon) pair is selected, the leading electron (muon) must pass the Tight identification [33] (Tight isolation) requirement. The invariant mass of the dilepton pair must be above 40 GeV in order to remove contributions from low-mass resonances.
The two highest-T photons in the event that pass the Tight identification and Loose isolation [33] requirements are selected. The two selected photons must be separated from each other by at least Δ = 0.4. Finally, the contribution from FSR photons is suppressed by requiring that the sum of the dilepton invariant mass and the smaller of the two three-body invariant masses, formed from the dilepton system and each of the two photons, is greater than twice the boson rest mass. This selection ( ℓℓ + min( ℓℓ 1 , ℓℓ 2 ) > 2 ) is illustrated in Figure 2.  2.10 simulation of the dilepton invariant mass versus the smaller of the two three-body masses formed from the dilepton system and each of the two photons. The events are subject to the full set of signal region requirements, with the exception of the FSR removal selection, which is indicated by the red line ( ℓℓ + min( ℓℓ 1 , ℓℓ 2 ) > 2 ).

Backgrounds
The dominant background contributions to ℓℓ arise from processes involving jets misidentified as photons (referred to as → backgrounds), and are estimated using a data-driven method. These backgrounds account for approximately 20% of the data yield in the signal region. A small contribution is expected from electrons misidentified as photons (referred to as → backgrounds), and is estimated from simulation. The remaining backgrounds, which are small, come from processes involving prompt photons, and are also estimated from simulation. The largest of the prompt-photon backgrounds arises from¯events, which contribute approximately 5% of the data yield in the signal region.
The number of data events selected in each channel is given in Table 2 along with the estimated background yields. The data yield in the muon channel is higher than in the electron channel because the muons have a higher reconstruction efficiency than electrons, and also a larger detector acceptance. Table 2 also shows the extracted number of signal events in data (i.e. data minus background), which is compared with predictions from both signal event generators for each channel. The detector-level 1 T , 2 T and ℓℓ data distributions, in both the electron and muon channels, are compared with the signal-plus-background predictions in Figure 3. The estimation of the different backgrounds is described in the following subsections.

→ backgrounds
Processes involving jets misidentified as prompt photons populate the ℓℓ signal region. They typically involve light-hadron decays into a pair of photons within jets. In these processes, one or both of the photon candidates are misidentified jets; these are divided into , and categories, the first two according to whether the lower-or higher-T photon candidate is a misidentified jet. The probability of a jet being misidentified as a photon is poorly modelled in simulation, so a data-driven method is used. Such methods utilise jet-enriched control regions, defined by using photon candidates which fail either the photon identification or isolation selections, or both. A loosened data sample is used, where two photons are selected with a loose identification requirement and with no isolation requirement. Each of the two photons in an event can be assigned to one of four categories defined by the signal region's photon identification and isolation requirements: pass identification and pass isolation (A), pass identification and fail isolation (B), fail identification and pass isolation (C), or fail identification and fail isolation (D). The signal region is hence denoted by AA, and the other 15 combinations define the control regions. In the following, the number of data (signal) events falling into each control region is denoted by data ( signal ) where , = A, B, C, D.
The number of events involving jets misidentified as photons in the signal region can be computed from the relevant yields in the control regions using a matrix method that has been employed in previous diphoton analyses [1,35]. The method uses as inputs the prompt-photon isolation efficiencies ( 1 , 2 ), which are the probabilities for Tight identified photons to be isolated, and the jet-to-photon fake rates ( 1 , 2 ), which are the probabilities for photon candidates which fail the Tight identification selection to be isolated. The real photon efficiencies are measured in simulation and the jet-to-photon fake rates are calculated in data as where the indices 1 and 2 refer respectively to the leading and subleading photons and = A D / B C is a correlation parameter which accounts for the bias due to requiring the photon candidate to fail the identification requirement. The parameter is determined from simulation for each of the photon candidates, and the combined average of = 1.18 ± 0.18 is used for the calculated values of both fake rates. The signal leakage into the jet-enriched control regions is corrected for by subtracting from data the number of signal events predicted by the simulation to fall in the control region. The number of events from each process in the loosened sample ( where , = , ) can be mapped onto the signal region (AA) and the control regions AB, BA and BB by using the matrix AA AB BA BB = 1 2 1 2 1 2 1 2 . (1) The matrix can then be inverted to determine the unknown yields, . The contributions of the four processes to the signal region are determined from the first row of the matrix in Eq. (1): The → background fractions in the signal region are determined after combining the data events from the + − and + − channels. The statistical uncertainty of the fractions is derived using 1000 sets of 'toy' data. For each set, the data yield in each region is randomly drawn from a Poisson distribution with a mean value equal to the observed data yield in that region. Each set of toy data is propagated through the matrix inversion, and the standard deviation of the 1000 extracted background fractions is taken as the statistical uncertainty. The systematic uncertainties related to the choice of control regions, the correlation parameter , and the variation of the photon isolation efficiencies with the photon T are considered. The largest contribution comes from the statistical uncertainty, but is similar in size to the total systematic uncertainty. The method is validated on a 'pseudo-dataset' formed from the signal and → background simulation samples, where the fractions of the four contributing processes are known and are reproduced accurately by the matrix method.
For the differential distributions, the shapes of the → backgrounds are taken from simulation and normalised to the overall fractions found in data. The shapes are derived in a slightly loosened signal region where one of the two photons is allowed to fail either the Tight identification or Loose isolation requirements, in order to increase the number of events selected from simulation. The ability of the simulation to describe the shapes in data is checked in a jet-enriched control region. Two sources of uncertainty affecting these shape templates are considered. The first is related to differences between data and simulation, and is estimated by testing the ability of the simulation to model the data distributions in a jet-enriched control region. The second is related to the choice of control region in which the shapes are derived, and is evaluated by reweighting the shapes to a harder T spectrum because the photons which fail the identification or isolation requirements are typically softer ones.

Other backgrounds
The second largest background contribution comes from¯events where the top quarks decay leptonically. The normalisation factor for this background is determined in a control region with the same selection requirements as the signal region, except that an opposite-sign lepton pair is selected. The¯process dominates in this region, but the contribution from → events is also considered using the matrix method described above. The ratio of data, with the → background subtracted, to the¯simulation in the control region is used to define a normalisation factor which is applied to¯simulation events entering the signal region. The considered sources of systematic uncertainty are the same as for the signal. The normalisation factor is 0.81 ± 0.17, where the dominant uncertainty is due to the limited number of data events in the control region.
As there are no vertex requirements placed on photons, a source of background arises when two protonproton interactions in the same bunch-crossing overlap to produce a combined ℓℓ system. A data-driven method, such as the one described in Ref.
[36], is not possible due to the limited number of signal region events, so these backgrounds are estimated using simulation only. Two processes contribute at first order: + and + . Random events from each sample are combined and subjected to the fiducial selection (described later in Section 7.1). The resulting particle-level distributions of the six kinematic variables listed in Section 1 are corrected to the detector level using bin-by-bin factors determined from the simulated signal events. A systematic uncertainty is assigned to account for the different T distributions of signal events and pile-up background events. It is estimated by recalculating the bin-by-bin factors after reweighting the signal simulation to the pile-up background photon T spectra. The uncertainties in the predicted cross-section of the single-photon [37] and diphoton [38] samples are significant and hence are also included as systematic uncertainties. The total systematic uncertainty is 35% (27%) for the + ( + ) pile-up background processes.
Another source of background is misidentification of electrons as photons. This → background is modelled by and simulations. The modelling of electron-to-photon misidentification rates has been tested [39] and is found to disagree with data at a level of up to 50% in some regions. Therefore, a conservative systematic uncertainty of 50% is applied to the and yields in the signal region.
The contribution from (→ ℓℓ) (→ ) is estimated directly from simulation. The contribution from (→ + − ) is estimated from simulation and is found to be negligible.

Systematic uncertainties
Systematic uncertainties in the measured cross-sections are related to the background estimation, the detector-to-fiducial acceptance correction factors (both inclusively and through the unfolding, as described in Section 7.1) and the integrated luminosity. The uncertainties in the backgrounds are discussed in Section 5. The correction factor and response matrix used for the unfolding are affected by the selection efficiency, and therefore variations of the different object reconstruction efficiencies are considered.
The performance of the electron and photon reconstruction, and their associated systematic uncertainties, are studied in Ref. [33]. For electrons, the reconstruction, identification and isolation efficiencies and their uncertainties are measured by applying tag-and-probe methods to events containing → + − or / → + − decays [40]. For photons, the corresponding efficiencies are measured using samples of → ℓ + ℓ − decays, and an inclusive photon sample collected using single-photon triggers. The energy scale and resolution for electrons and photons, and their uncertainties, are obtained from a sample of → + − events. For muons, the efficiencies, the momentum scale and resolution, and their uncertainties, are obtained using samples of → + − or / → + − decays [34].
The uncertainty due to the pile-up reweighting procedure discussed in Section 3 is estimated by varying the amount of pile-up in the signal simulation to cover the uncertainty in the ratio of the predicted and measured inelastic cross-sections [41].
The statistical uncertainty due to the limited number of generated signal events is considered.
Various sources of theoretical uncertainty are considered. The uncertainty due to the choice of PDF is estimated from the standard deviation of the mean of 100 variations of the nominal PDF set (NNPDF3.0 ). The renormalisation and factorisation scales are each varied by factors of 0.5 and 2.0, except for shifts in opposite directions, and the envelope of the effects of these scale variations is taken as an estimate of the uncertainty due to missing higher order corrections. The assumed value of the strong coupling constant, s ( ) = 0.118, is varied by ±0.001 and the average effect is taken as the s contribution to the uncertainty. The effect of these theoretical uncertainties is accounted for in the integrated fiducial cross-section measurements and in the predicted cross-sections from the MC event generators. For the differential cross-section measurements, the theoretical uncertainties are covered by the unfolding uncertainty (described in Section 7.2).
The systematic uncertainties in the integrated cross-section in the fiducial region are summarised in Table 3. The measurement in each channel is dominated by the data statistical uncertainty, and the largest systematic uncertainty comes from the → background estimation.

Fiducial volume definition
The measured yields for the signal process in data are corrected to a fiducial volume which accounts for detector inefficiency, geometry and resolution. The fiducial volume is defined using particle-level objects in simulation which have a proper decay length longer than 10 mm. To correct for bremsstrahlung, each particle-level lepton is 'dressed' by vectorially adding to its four-momentum the four-momenta of any nearby photons, except those from hadron decays, within a cone of size Δ = 0.1 around the lepton. To minimise the model-dependence of the procedure to correct the data from detector level to particle level, also known as unfolding, the selection requirements placed on the particle-level objects are as close as possible to the detector-level selection outlined in Section 4. An exception is the calorimeter transition region, which is included in the selection of particle-level electrons and photons. The detector-to-fiduciallevel correction procedure includes an extrapolation, over a few percent of the total phase space, which accounts for the loss of detector-level acceptance in this region. The particle-level photons must pass an isolation selection which requires iso T , defined as the summed transverse energy of all particles except muons, neutrinos and the photon itself within a cone of size Δ = 0.2 around the photon, to be less than 7% of the photon T . This value is chosen as it best replicates the performance of the Loose isolation working point used in the detector-level selection. The complete set of requirements is listed in Table 4.

Cross-section extraction
The integrated fiducial cross-section, fid , is calculated from the observed yield in data ( data ), the expected background yield ( bkg ) and the total integrated luminosity ( ) fid = data − bkg × .

Photons Leptons
The correction factor ( ) is defined as the ratio of the number of signal simulation events passing the detector-level selection to the number which pass the fiducial-level selection. The value of is 0.286 ± 0.014 in the electron channel and 0.379 ± 0.017 in the muon channel where the uncertainties include the systematic sources discussed in Section 6.
For the differential cross-section measurements, the detector-to-fiducial-level correction is instead done via an iterative Bayesian unfolding procedure [42], with two iterations, which accounts for migrations between bins. These migrations are typically below 5% but can be as large as 18% in the regions with the finest binning. The reliability of the unfolding procedure is tested by unfolding the detector-level signal distribution from simulation, reweighted such that it better describes the data. The difference between the resulting unfolded distribution and the reweighted fiducial distribution is assigned as an uncertainty of the differential measurements. The uncertainty is negligible in most bins, up to 8% in one bin, but overall has a very small effect on the measurements. The integrated and differential cross-section measurements are performed separately in each channel. The results are combined into (→ ℓℓ) measurements via an averaging procedure which accounts for statistical and systematic uncertainties, and their correlations between the two channels. The technique uses a 2 minimisation procedure which is documented in Ref. [43].

Results
The measured integrated cross-sections in each channel and the combined average are The integrated (→ ℓℓ) cross-section is measured with an overall precision of 12% and is compared with the MC event generator predictions in Figure 4, where good agreement between data and both predictions is seen. The S prediction suffers from larger scale uncertainties due to the matrix element NLO accuracy being at the 0-jet level, whereas it includes 1-jet contributions in the M G 5_ MC@NLO prediction.
The measured (→ ℓℓ) differential cross-sections are compared with the predictions from S and M G 5_ MC@NLO in Figures 5, 6 and 7.
The photon transverse energy ( 1 T , 2 T ) distributions displayed in Figure 5 are well described by the predictions. The ℓℓ T distribution in Figure 6(a) describes the transverse momentum of the boson, which typically recoils against the two photons. This distribution is therefore sculpted by the minimum transverse momentum selections imposed on the two photons, which results in the peak around 40 GeV. The ℓℓ T distribution in Figure 6(b) probes the QCD modelling of the transverse momentum of the system. The description by the MC event generator predictions is good with the exception of a discrepancy in the high-T region that shows the impact of the limited accuracy of the predictions in the presence of jets. The distribution is important in the context of diphoton resonance searches in ℓℓ channels. The measured distribution is shown in Figure 7(a) and the simulations provide a good description, particularly in the fourth bin, which is most relevant for (→ ℓℓ) (→ ) measurements. The ℓℓ distribution (Figure 7(b)) provides a measure of the hard scale of the system, and is described well by the predictions, even for ℓℓ > 500 GeV.    Figure 5: The unfolded differential cross-sections as a function of (a) the leading photon transverse energy and (b) the subleading photon transverse energy are compared with NLO predictions from S and M G 5_ MC@NLO. The black uncertainty bar on the data represents the statistical uncertainty, whereas the slightly taller grey uncertainty band represents the total uncertainty. The uncertainty in the predictions includes both the statistical and theoretical uncertainties. The lower panels show the ratios of the predictions to data, as well as the fractional uncertainty of the data.  Figure 6: The unfolded differential cross-sections as a function of (a) the dilepton transverse momentum and (b) the four-body transverse momentum are compared with NLO predictions from S and M G 5_ MC@NLO. The black uncertainty bar on the data represents the statistical uncertainty, whereas the slightly taller grey uncertainty band represents the total uncertainty. The uncertainty in the predictions includes both the statistical and theoretical uncertainties. The lower panels show the ratios of the predictions to data, as well as the fractional uncertainty of the data.  Figure 7: The unfolded differential cross-sections as a function of (a) the diphoton invariant mass and (b) the four-body invariant mass are compared with NLO predictions from S and M G 5_ MC@NLO. The black uncertainty bar on the data represents the statistical uncertainty, whereas the slightly taller grey uncertainty band represents the total uncertainty. The uncertainty in the predictions includes both the statistical and theoretical uncertainties. The lower panels show the ratios of the predictions to data, as well as the fractional uncertainty of the data.

EFT interpretation
The ℓℓ final state can probe the non-Abelian structure of the SU(2) L × U(1) Y symmetry of the Standard Model, which gives rise to gauge boson self-interactions. Modifications of the self-interactions arising through new physics (NP) are investigated using an effective field theory approach [44]. The SM Lagrangian L SM is expanded with operators of dimension > 4, which are suppressed by the energy scale Λ of NP:  Figure 8. The SM contribution estimated by a M G 5_ MC@NLO simulation at LO is also shown and this can be compared with the NLO prediction to investigate the impact of NLO QCD corrections. A slightly softer ℓℓ T spectrum is observed at LO. The predicted differential cross-section at LO is used in Section 8.1 in the estimation of NLO QCD corrections for the EFT prediction.
Limits are placed on the coupling parameters , /Λ 4 by constructing and scanning a profile likelihood ratio, taking as input the baseline S 2.2.10 production (expected limits) and Run 2 data (observed limits), the contributions of the transverse operators, and all sources of uncertainty. The limits are calculated for the combination of the electron and muon channels. In the fit, the likelihood function is represented by a multivariate Gaussian distribution, where theory uncertainties are modelled by additional Gaussian constraints. All experimental uncertainties are encoded in the covariance matrix accounting for ℓℓ T bin correlations. A description of the experimental uncertainties is given in Section 6. The bin-to-bin correlation of the statistical uncertainty, which can be present after unfolding, is found to be negligible and is not considered further. The shift of each systematic uncertainty is applied in a fully correlated way between all bins; correlations between different sources of systematic uncertainty are not considered. Theory uncertainties, consisting of renormalisation and factorisation scale, PDF, and s uncertainties, are included for SM production and all transverse operators. Gaussian constraints for the limited size of the S 2.2.10 sample and the EFT samples are also added. The largest experimental uncertainties stem from the limited size of the data sample and the estimation of the fake-photon contribution, reaching 18% and 14%, respectively, in certain ℓℓ T bins. The renormalisation and factorisation scales are the sources of the largest theory uncertainties, which reach values of 23% in the last ℓℓ T bin. Limits at 95% confidence level are constructed from the profile likelihood ratio by applying Wilks' theorem [46] and thus assuming that the test statistic is 2 -distributed. The effect of one transverse operator at a time is studied while all other Wilson coefficients are set to zero.

Non-unitarised limits
The expected and observed limits are displayed in Figure 9. Constraints arising from unitarity conservation are not considered. The observed limits are typically 11%-12% less stringent than those expected. This is driven by the larger contribution of fake photons in data and the corresponding uncertainties in the fake-photon normalisation and shape.
Higher-order QCD corrections are not available for the EFT prediction. In order to study the impact of such corrections on the constraints that can be placed on couplings of dimension-8 operators, a test fit was performed assuming that the EFT scales similarly to the SM with respect to higher-order corrections. In this test, the M G 5_ MC@NLO 2.7.3 differential cross-section at NLO (see Table 1) was divided by the M G 5_ MC@NLO LO prediction displayed in Figure 8 to obtain bin-wise correction factors. The parameter settings for the LO simulation were identical to those chosen for the generation of the EFT contributions, except that all Wilson coefficients were set to zero. The differential cross-sections predicted by O 8 , were then multiplied by the correction factors. The results of this study show that the expected and observed constraints on the coupling parameters , /Λ 4 are 13%-15% more stringent. Although such higher-order QCD corrections result in a sizeable impact on the limits, the underlying assumption cannot be validated with the available theoretical calculations in the EFT formalism, therefore the nominal confidence intervals, as shown in Figure 9, are calculated without the bin-wise correction factors.
The confidence intervals presented for the four transverse operators O 8 ,1 , O 8 ,2 , O 8 ,6 , O 8 ,7 are the first published by the ATLAS experiment at a centre-of-mass energy of 13 TeV and are up to two orders of magnitude more stringent than the limits extracted at 8 TeV. The non-unitarised limits derived in this analysis are less stringent than those published by CMS in their ± and analysis at 13 TeV [3], but differ by less than a factor of two. This difference is primarily driven by binning requirements on ℓℓ T . The binning was optimised to have sufficient events in the unfolding procedure. This results in a reduced sensitivity to EFT effects, particularly in the final bin.

Unitarisation treatment
Limits for all transverse operators are also derived as functions of an energy scale cut-off which prevents the violation of unitarity at large energy scales. Various techniques make use of a truncation of the EFT contribution to restore unitarity. This analysis uses a method which is typically referred to as clipping [47]. Any EFT contribution is suppressed above an energy scale c . The ℓℓ invariant mass is used to select various thresholds for the clipping energy. Technically, this is achieved by scanning the simulated events at parton level before the parton showering is performed and suppressing any EFT contribution of events in which ℓℓ > c . The SM contribution is not truncated and is allowed to reach arbitrary energy scales. The evolution of the expected and observed confidence intervals as a function of c for the coupling parameter of transverse operator O 8 ,8 is shown in Figure 10. A similar evolution is observed for the remaining dimension-8 operators, where the confidence intervals become 4-5 more stringent between c = 1.1 TeV and c = ∞.  Figure 8: Comparison of the differential cross-section at particle level as a function of the dilepton transverse momentum between the observation in full Run 2 data, the NLO prediction from S , the LO prediction from M G 5_ MC@NLO, and the EFT prediction of one dimension-8 operator. The black uncertainty line on the data represents the statistical uncertainty, whereas the slightly taller grey uncertainty band represents the total uncertainty. The uncertainty in the S prediction includes both the statistical and theoretical uncertainties.

Conclusions
The production of a boson in association with two photons in a phase-space region dominated by the ISR production of photons is studied in proton-proton collisions. The measurements are performed using 139 fb −1 of 13 TeV proton-proton collision data recorded by the ATLAS detector at the LHC. The electron and muon decay channels of the boson are used, and events where either photon is radiated from one of the final-state leptons are rejected.
The integrated (→ ℓ + ℓ − ) cross-section is measured with a precision of 12%, with approximately equal contributions from statistical and systematic uncertainties. The cross-section is also measured differentially for the first time, and is used to test Standard Model predictions at up to next-to-leading-order accuracy from S and M G 5_ MC@NLO. The description by the MC event generator predictions is good with the exception of a discrepancy in the high-ℓℓ T region that shows the impact of the limited accuracy of the predictions in the presence of jets.
The measurements are also used to set limits on the Wilson coefficients, divided by the new physics scale Λ, of dimension-8 EFT operators. The constraints on four of the eight operators under consideration are tightened by up to two orders of magnitude with respect to previous ATLAS analyses using 8 TeV data.