Measurements of the pp $\to$ W$^\pm\gamma\gamma$ and pp $\to$ Z$\gamma\gamma$ cross sections at $\sqrt{s} =$ 13 TeV and limits on anomalous quartic gauge couplings

The cross section for W or Z boson production in association with two photons is measured in proton-proton collisions at a centre-of-mass energy of 13 TeV. The data set corresponds to an integrated luminosity of 137 fb$^{-1}$ collected by the CMS experiment at the LHC. The W $\to$ $\ell\nu$ and Z $\to$ $\ell\ell$ decay modes (where $\ell =$ e, $\mu$) are used to extract the W$\gamma\gamma$ and Z$\gamma\gamma$ cross sections in a phase space defined by electron (muon) with transverse momentum larger than 30 GeV and photon transverse momentum larger than 20 GeV. All leptons and photons are required to have absolute pseudorapidity smaller than 2.5. The measured cross sections in this phase space are $\sigma$(W$\gamma\gamma$) = 13.6 $^{+1.9}_{-1.9}$ (stat) ${}^{+4.0}_{-4.0}$ (syst) $\pm$ 0.08 (PDF + scale) fb and $\sigma$(Z$\gamma\gamma$) = 5.41 $^{+0.58}_{-0.55}$ (stat) ${}^{+0.64}_{-0.70}$ (syst) $\pm$ 0.06 (PDF + scale) fb. Limits on anomalous quartic gauge couplings are set in the framework of an effective field theory with dimension-8 operators.


Introduction
The measurement of the associated production of a vector boson V (= W, Z) and two photons in proton-proton (pp) collisions is a powerful test of the standard model (SM). The nonabelian nature of the electroweak interaction predicts the presence of self-interacting vector boson vertices. The strength of the interaction is set by the values of triple and quartic gauge couplings predicted by the SM. The measurement of possible deviations from the theoretical predictions could provide indirect evidence of new particles or new interactions. Discrepancies at high photon momentum, where new physics might give a measurable deviation from the SM cross section, would produce evidence for the possible existence of anomalous quartic gauge couplings (aQGCs). A parametrisation of predictions involving anomalous couplings, independent of any specific new physics model, can be calculated in an effective field theory (EFT) framework [1]. Triboson production is also an important background for several SM and beyond the SM processes, such as the Higgs boson production in association with vector bosons (with H → γγ). Thus, studies of the self-couplings of the electroweak gauge bosons provide an excellent opportunity for a deeper understanding of electroweak interactions.
Some of the elementary processes resulting in the production of a massive vector boson in association with two photons at the CERN LHC are presented in the leading-order (LO) Feynman diagrams of Fig. 1. Topologies with final states that originate from events where the V boson is produced in the hard interaction between the two partons, and the photons come from either initial or final state radiation processes or from quartic gauge vertices together with the V boson, are studied, and are referred to as Vγγ in the text.
Previous measurements of the Vγγ production cross sections have been performed by the ATLAS and CMS Collaborations at the CERN LHC in pp collisions at a centre-of-mass energy of √ s = 8 TeV [2][3][4]. Limits on the presence of aQGCs were also reported in these papers.
In this paper, the first measurements of the pp → Wγγ and pp → Zγγ cross sections at √ s = 13 TeV are presented using data collected between 2016 and 2018 by the CMS experiment, corresponding to an integrated luminosity of 137 fb −1 . For these measurements, only direct decays into electrons or muons are considered. Measurements are compared with the latest available calculations at next-to-LO (NLO) in perturbative quantum chromodynamics (QCD) [5][6][7]. Limits on the aQGCs are presented in the framework of an electroweak EFT with dimension-8 operators.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity η coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionisation chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system [8]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz with a latency of about 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimised for fast processing that reduces the event rate to around  1 kHz before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, is reported in Ref. [9].

Event simulation
The associated production of a W (Z) boson and at least two photons is searched for in events with one lepton (two opposite-sign, same-flavour leptons) and two photons. Only electron and muon decay channels are used, while the τ decays are treated as a background. The Vγγ signal samples are generated at NLO with MADGRAPH5 aMC@NLO [10] (2.2.2 for 2016 samples, 2.2.6 for 2017 and 2018 samples) with no additional jets in the matrix element calculation. The NLO NNPDF 3.0 set [11] (for 2016 samples) and the next-to-NLO NNPDF 3.1 set [12] (for 2017 and 2018 samples) are used as parton distribution function (PDF) sets.
The main background contribution comes from the misidentification of jets as photons, which is estimated in single-photon control regions. Thus, various background samples involving a single photon are needed, including V(V)γ samples. Other single-photon and diphoton processes (such as tt produced in association with one or two photons or a photon and a jet) contribute as backgrounds and are estimated using Monte Carlo (MC) simulations. The background contribution from the associated production of a W(Z) boson and a Higgs boson decaying to two photons is negligible, and is not considered. is used to separate genuine electrons from misidentified ones. A tight identification is used to select prompt electrons (produced at the primary vertex) and isolated electrons in the final state [20]. Background contributions from misidentified jets or electrons inside a jet are rejected applying electron isolation criteria, which exploit the PF-based event reconstruction. The electron isolation variables are obtained by summing the p T of charged hadrons compatible with the primary vertex I chg , of neutral hadrons I neu , and of photons I pho inside a cone of radius ∆R = √ (∆η) 2 + (∆φ) 2 = 0.3 around the electron direction, where φ is the azimuthal angle in radians. Additional photons and neutral hadronic contributions to the isolation variable, coming from pileup, are subtracted using the jet area approach [21].
Muons candidates are required to have p T > 15 GeV in the pseudorapidity range |η| < 2.4. Muon identification criteria are based on the fit quality for tracks measured in the tracker and muon detectors. A tight muon identification is used to reconstruct muons in the final state [22]. To distinguish between prompt muons and those from hadron decays within jets, muons are required to be isolated with respect to all nearby PF reconstructed particles. For the computation of the PF isolation, the I chg , I neu and I pho components are summed in a cone of ∆R = 0.4 around the muon direction. The corrected energy sum is obtained by subtracting the pileup contribution to I neu and I pho , which is estimated as half of the corresponding charged hadronic component.
Photons are selected with p T > 20 GeV in the pseudorapidity range |η| < 1.44 and 1.57 < |η| < 2.5. Photon identification is based on the sequential application of several selections. A medium photon identification is used to reconstruct prompt photons (i.e. not from hadron decays) in the final state [20]. The average efficiency for this selection is 80%. Photons selected in this analysis are required to have a narrow transverse shape of the electromagnetic shower, a minimal energy deposit in the HCAL, and to be isolated with respect to other particles. The same isolation variable previously described for the electron selection is used.
The reconstruction, identification, and isolation efficiencies of leptons and photons and the trigger efficiencies of leptons are measured with the "tag-and-probe" technique [23], as a function of particle η and p T in both data and simulation. A sample of events containing a Z boson decaying into e + e − or µ + µ − is used for these measurements. The photon efficiencies are derived using a sample of electrons from Z decays with no requirement on the track and charge of the candidate. These efficiencies are used to correct for the differences between data and simulation.
An event is categorised as a W decaying to leptons if exactly one electron (muon) with p T > 35 (30) GeV is selected. The selected lepton must match the one that triggered the event and must be associated with the primary vertex of the collision. If the event contains any additional different-flavour leptons or opposite-sign same-flavour leptons, it is excluded from the W boson candidate sample, but is further checked for the presence of a Z boson candidate.
An event is categorised as a Z boson candidate decaying to leptons if two opposite-sign leptons of the same flavour are selected. The leading p T electron (muon) is required to have p T > 35 (30) GeV, and the subleading electron or muon is required to have p T > 15 GeV. Only the leading lepton is required to match the one that triggered the event, although both are required to be associated with the primary vertex of the collision. The invariant mass of the dilepton system is required to be m > 55 GeV.
Events are selected if they have a single W or Z boson candidate and at least two photons. All reconstructed photons must be separated from each other and from each reconstructed lepton by ∆R > 0.4. Photons are discarded if |m e,γ − m Z | < 5 GeV (where m e,γ is the invariant mass of an electron and the leading or subleading photon and m Z is the Z boson invariant mass) or if |m eγγ − m Z | < 5 GeV (where m eγγ is the invariant mass of an electron and the two photons). In this way, photons likely coming from final-state bremsstrahlung radiation are removed and, therefore, the contribution from electrons misidentified as photons is reduced as well.

Background estimation
The backgrounds in both the Wγγ and Zγγ signal regions are categorised as events with a genuine photon or with another object misidentified as a photon. The largest contribution in both channels comes from the misidentification of jets as photons. Another important source of background originates from electrons that are reconstructed as photons because the deposit in the calorimeter is not associated with a track in the tracker. This contribution is particularly relevant in the Wγγ electron channel, and is dominated by the Zγ electron channel. Both of these background processes are estimated by exploiting a control sample in data. The remaining minor contributions from processes that have at least one genuine photon (tγ, tt γ, tt γγ, Vγ and VVγ) are estimated using MC simulations and referred to as "others".
The background from events containing nonprompt photons is estimated following a method similar to the ones described in Refs. [3,4]. A W or a Z boson is selected together with one photon that passes the standard selection except for the isolation requirement both in data and simulation. Events are categorised as "tight" or "loose" if the photon candidate passes or fails the isolation requirement. The probabilities for photon and for a jet f to be isolated are computed as where N T γ, MC (N L γ, MC ) is the number of simulated events with a tight (loose) photon while N T γ, data (N L γ, data ) is the number of events with a tight (loose) photon candidate in data after the subtraction of the prompt photon contribution from simulation. These probabilities are calculated separately for photons in the ECAL barrel and endcap regions as a function of the photon p T . The jet-photon misidentification background in the diphoton phase space is then estimated by solving the system where the indices of the and f coefficients refer to the leading and the subleading photon. The N XY vector contains the number of events where two (TT), one (TL and LT) or zero (LL) photon candidates pass the isolation requirement. The α AB vector contains the number of signal (γγ) and background (γj, jγ and jj) events. This method is validated in a control region enriched in the jet-photon misidentification background where both photons fail the isolation selection. By solving the matrix for each combination of the single photons p T and η and by applying it to the number of events in the double-photon phase space in each diphoton p T bin one determines the number of background events in the signal region as Because of the large contamination from Zγ → eeγ events where an electron from the Z boson decay is misclassified as a photon, the jet misidentified as a photon contribution for the Wγγ electron channel is estimated from the and f coefficients that are evaluated from the singlephoton Wγ muon sample. For the same reason, the Zγ → eeγ background contribution has been subtracted from data in the Wγγ electron channel before computing the number of background events in the signal region.
To estimate the contribution where an electron is misclassified as a photon, a correction factor is computed from Zγ events and is then applied to the simulation. The invariant mass of an electron and a photon is reconstructed in data and MC simulation while removing the requirement |m e, lead γ − m Z | < 5 GeV. This mass distribution is fitted with the sum of a signal template, derived from MC simulation, and a background function, which has an exponential decay distribution at high mass (above the Z peak) and a turn-on (linked to the electron and photon p T thresholds) described by an error function at low mass. A correction factor, obtained in intervals of the photon p T and η, is then computed as where N data inv is the number of events in the electron-photon invariant mass peak obtained by integration of the fitted signal shape and N data Z is the number of events for the Z → ee invariant mass distribution obtained by integration of a fitted double-sided Crystal-Ball function [24] in the data. The same procedure is used to calculate the number of events in MC simulation. By fitting all the distributions in the different (p T , η) bins, a set of correction factors is computed. All the MC simulations are then corrected for these factors on an event-by-event basis whenever a reconstructed photon matches a generator-level electron. The event-by-event correction factors are on average about 20%.
The pre-fit (i.e. before the fitting procedure described in Section 7) diphoton p T distributions for the Wγγ and Zγγ analyses, separated in the electron and muon channels, are shown in Fig. 2. The same distributions obtained in the control region enriched in jets misidentified as photons are shown in Fig. 3 for the Wγγ and Zγγ electron channels. The data and prediction agree thus validating the jet-photon misidentification background estimation procedure. A similar level of agreement is observed in the muon channel.

Systematic uncertainties
Systematic effects can affect the rates and distributions of both data and simulation. To estimate these uncertainties, the full analysis procedure is repeated by varying each quantity by plus or minus its standard deviation uncertainty. In this procedure, correlations between the systematic uncertainties are included where appropriate.
The dominant systematic uncertainties come from the estimation of the backgrounds. To determine the systematic uncertainty coming from the jet-photon misidentification background, the same strategy is applied to a QCD control sample that is obtained using the Wγ selection but inverting the isolation requirement on the leptons while keeping the photon selection identical to the one for the signal region. This sample is used to obtain an alternative estimate of the jet-photon misidentification background contribution in the W channel. For the Z channel, the QCD control sample resulting from the Zγ selection with the inversion of the lepton isolation has insufficient events. Hence, a transfer factor from the Wγ selection is computed and applied for the determination of the alternative estimate of the jet-photon misidentification background contribution in the Z channel. The systematic uncertainty is computed as half the difference in the distributions between the standard method or the one just described.     The black points represent the data with error bars showing the statistical uncertainties. The hatched histogram shows the expected signal contribution. The background estimate for electron (jet) misidentified as photons, obtained from control samples in data, is shown in yellow (purple). The remaining background, derived from MC simulation, is shown in green. In the ratio plots, the grey hashed area is the statistical uncertainty on the sum of signal and backgrounds, while the uncertainty in the black dots is the statistical uncertainty of the data. In blue, the expected distribution for an example value of the anomalous coupling parameters f M3 /Λ 4 and f T0 /Λ 4 is also shown (see Section 8 for the details).
Another source of uncertainty in the jet-photon background is related to the modelling of the initial-and final-state radiation and of the energy spectra of the final state particles. An alternative Vγ MC simulation, obtained with SHERPA, is used to evaluate this uncertainty.
The uncertainty in the correction factor related to the background of electrons misidentified as photons is determined by propagating the estimated uncertainty in the correction factor F of Eq. 1. This has two components: a statistical one, that comes from the uncertainty in the fitting procedure; and a systematic one that is computed by taking half the difference between the F factors obtained by performing the fit of the signal component with a double-sided Crystal-Ball function and with the nominal method using an MC template.
The uncertainties in the lepton and photon reconstruction and selection efficiencies are included by computing the cross section with these efficiencies varied up and down by one standard deviation. The uncertainty related to these data-to-simulation corrections is estimated by   The uncertainty in the value of the theoretically computed cross section is accounted for during the subtraction from data of the background processes estimated from MC. Furthermore, the value of the expected cross section has an impact on the estimation of the jet-photon background because the contribution from prompt photons is subtracted from the distribution in data using the Wγ and Zγ simulations. To estimate these contributions, the cross sections of the Wγ, Zγ, and of the other minor backgrounds are varied independently. The uncertainty in the Zγ cross section is estimated as half the difference between the next-to-NLO and the NLO values computed with MATRIX [25], and amounts to 2.5%. The same uncertainty is assumed for the Wγ cross section, and a value of 7.5% is used for the other simulated backgrounds.
The total inelastic cross section is varied by 4.6% [26] to estimate the impact on the final result of the pileup reweighting procedure. The uncertainty because of the integrated luminosity measurement is equal to 2.5, 2.3, and 2.5% for the 2016, 2017 and 2018 data taking periods, respectively [27-30]. Because of the uncorrelated time evolution of some systematic uncertainties, the total integrated luminosity has an uncertainty of 1.8% and is applied to all the processes estimated with an MC simulation. The effect of the uncertainty in the integrated luminosity affects the estimation of the jet-photon misidentification background as well as the MC contributions in the diphoton distributions.
For the extraction of the results, each systematic uncertainty is represented by a nuisance parameter, which affects the shape and the normalisation of the distribution of the various background contributions. The variation of the nuisance parameter results in a continuous perturbation of the spectrum, following a Gaussian probability density function. The impact of each systematic uncertainty is obtained by freezing the set of associated nuisance parameters to their best-fit values and comparing the total uncertainty in the measured cross section with   Table) and Zγγ (lower Table) selections in the electron and muon channels. The systematic uncertainties of the individual backgrounds and the total background are obtained by summing the contributions of different systematic uncertainties in quadrature. The statistical uncertainties are those related to the MC event samples and control region statistical uncertainties.

Cross section measurements
The cross sections for the Wγγ and Zγγ processes are measured separately in the electron and muon channels using a sample of events corresponding to an integrated luminosity of 137 fb −1 (LHC Run 2 data). The observed and predicted numbers of events are presented in Table 2.
The measured yields in the electron and muon channels are extrapolated to a common fiducial phase space determined from simulated signal events at the generated particle level. Generated particles are considered stable if their mean decay length is larger than 1 cm. Generated leptons are required to have a p T > 15 GeV and |η| < 2.5. The momenta of photons in a cone of ∆R = 0.1, the same cone size as the one applied to reconstructed data, are added to the charged lepton momentum to correct for final-state radiation. Generated photons are required to have p T > 15 GeV and |η| < 2.5. Additionally, the candidate photons are required to have no selected leptons or photons in a cone of radius ∆R = 0.4 and no other stable particles, apart from photons and neutrinos, in a cone of radius ∆R = 0.1. Events are then selected in the Wγγ channel by requiring exactly one electron (muon) with p T > 30 GeV and at least two photons with p T > 20 GeV. Events are selected in the Zγγ channel by requiring two electrons (muons), at least one of them with p T > 30 GeV, and not less than two photons, each of them with p T > 20 GeV. Additionally, the invariant mass of the dilepton system is required to be m > 55 GeV.
The expected theoretical cross sections are predicted at NLO and their uncertainties come from the finite MC sample event count used to compute them, the PDF set, the factorisation and renormalisation scales. Statistical uncertainties are estimated to be of the order of 0.2% in both the Wγγ and Zγγ channels. Uncertainties related to the PDF set are estimated using a set of 100 replicas of the NNPDF 3.1 PDF set, following the Ref. [12] prescription, and are estimated to be of the order of 0.3% in the eνγγ and µνγγ and of 0.8% in the eeγγ and µµγγ channels. Uncertainties related to the renormalisation and factorisation scale choice are estimated by independently varying µ R and µ F by a factor of 0.5 and 2, with the condition that 1/2 < µ R /µ F < 2. The uncertainties are defined as the maximal differences from the nominal values and are estimated to be of the order of 0.6 (0.5)% in the eνγγ (µνγγ) and of 0.6 (0.7)% in the eeγγ (µµγγ) channels. Uncertainties related to the value of the strong coupling are estimated to be of the order of 0.03 (0.02)% in the eνγγ (µνγγ) and of 0.4% in the eeγγ and µµγγ channels.
Binned maximum likelihood fits to the diphoton p T distributions in Fig. 2   The results of the fit for the electron and muon channels separately are compatible within two sigmas. In the combined electron and muon channel, the best fit value for the Wγγ signal strength is 0.73 +0.10 −0.10 (stat) +0.22 −0.22 (syst) and for the Zγγ signal strength is 0.91 +0.10 −0.09 (stat) +0.11 −0.12 (syst). The measured cross sections are: σ(Wγγ) meas = 13.6 +1.9 −1.9 (stat) +4.0 −4.0 (syst) ± 0.08 (PDF + scale) fb, σ(Zγγ) meas = 5.41 +0.58 −0.55 (stat) +0.64 −0.70 (syst) ± 0.06 (PDF + scale) fb. The sensitivity for the Wγγ cross section measurement is dominated by the muon channel. The measured signal strengths are summarised in Fig. 4. The significance of the cross section measurement for both the Wγγ and Zγγ channels is quantified using the background-only hypothesis under the asymptotic approximation [34]. The observed (expected) significance for the Wγγ signal is 0.6 (2.7) σ in the electron channel and 3.0 (4.3) σ in the muon channel; for the Zγγ is 3.4 (5.0) σ in the electron channel and 5.4 (5.1) σ in the muon channel; the combined significance for the Wγγ is 3.1 (4.5) σ and for the Zγγ is 4.8 (5.8) σ.

Limits on anomalous quartic gauge couplings
Studies of the anomalous gauge couplings can be performed in the EFT framework [1] by expanding the SM Lagrangian to include terms with dimension higher than four. In particular, both the Wγγ and Zγγ processes are sensitive to the presence of dimension-6 and dimension-8 operators [35]. Because of the available statistics in the Vγγ channel, the sensitivity to dimension-6 operators is expected to be lower than the one in the diboson production. The contribution of each operator is proportional to a coupling constant f and to the inverse of the energy scale Λ at which the new phenomena appear.
In the generation of the anomalous couplings samples, a calculation using 10 (8) different dimension-8 operators was performed for the Wγγ (Zγγ) process. The operators can be divided into two subsets: the L M0 -L M7 ones, that contain both the SU(2) L and U(1) Y field strengths and the covariant derivative of the Higgs doublet, and the L T0 -L T9 ones, that contain only the two field strengths. In particular, the Wγγ channel is especially sensitive to the M2 , M3, T0, T1, T2, T5, T6, and T7 operators, whereas the Zγγ channel is especially sensitive  to the T0, T1, T2, T5, T6, T7, T8, and T9 operators. The distribution of the p T of the diphoton system (shown in Fig. 2) is used to constrain the aQGC parameters under the hypothesis of absence of anomalies in triple gauge couplings. The contribution of aQGCs is enhanced at high values of the p T of the diphoton system. The distribution of the aQGCs as a function of the couplings themselves has a quadratic behaviour, and hence a parabolic fit is implemented to interpolate between the different values obtained via the parameter scan. The fitting procedure is performed bin-by-bin to exploit the shape of the distribution to set the limits and include the different systematic uncertainties. To further increase the sensitivity, electron and muon channels are combined. Each operator coefficient is scanned independently with all other operators set to zero. The extraction of the 95% confidence level upper and lower limits on the aQGCs is performed by exploiting the procedure described in Ref. [33]. The expected and measured limits for both the Wγγ and Zγγ processes are presented in Table 3. In

Summary
The cross sections for both the Wγγ and Zγγ processes are measured in proton-proton collisions by the CMS experiment at a centre-of-mass energy of 13 TeV corresponding to an integrated luminosity of 137 fb −1 .
The cross sections are measured in a fiducial region where simulated signal events are selected at generator level in the Wγγ channel by requiring exactly one electron or muon with transverse momentum p T > 30 GeV and at least two photons, each with p T > 20 GeV. Events are selected in the Zγγ channel by requiring two oppositely charged electrons or muons, at least one of them with p T > 30 GeV, and at least two photons, each with p T > 20 GeV. All leptons and photons are required to have pseudorapidity |η| < 2.5. Additionally, the invariant mass of the dilepton system is required to exceed m > 55 GeV.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid and other centres for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies:  [22] CMS Collaboration, "Performance of the CMS muon detector and muon reconstruction with proton-proton collisions at √ s = 13 TeV", JINST 13 (2018) P06015, doi:10.1088/1748-0221/13/06/P06015, arXiv:1804.04528.