Search for resonances decaying into photon pairs in 139 fb$^{-1}$ of $pp$ collisions at $\sqrt{s} =$ 13 TeV with the ATLAS detector

Searches for new resonances in the diphoton final state, with spin 0 as predicted by theories with an extended Higgs sector and with spin 2 using a warped extra-dimension benchmark model, are presented using 139 fb$^{-1}$ of $\sqrt{s} = $ 13 TeV $pp$ collision data collected by the ATLAS experiment at the LHC. No significant deviation from the Standard Model is observed and upper limits are placed on the production cross-section times branching ratio to two photons as a function of the resonance mass.


Introduction
Many theories of physics beyond the Standard Model (SM) predict the presence of new high-mass states which can currently only be produced in the high-energy collisions at the Large Hadron Collider (LHC). One powerful experimental signature of these new states is a resonance in the diphoton invariant mass spectrum, as this final state provides excellent invariant mass resolution which can be used to separate the signal from the SM background processes.
This Letter presents a search for new high-mass resonances decaying into two photons using 139 fb −1 of √ = 13 TeV proton-proton ( ) collision data recorded by the ATLAS detector from 2015 to 2018 at the LHC. This analysis searches for a generic resonance using two benchmark signal models: a spin-0 resonant state ( ), which is predicted by many models that include extensions to the Higgs sector [1][2][3][4][5][6][7]; and a spin-2 graviton ( * ), taken here to be the lightest Kaluza-Klein (KK) [8] excitation in a Randall-Sundrum model [9,10] with one warped extra dimension (RS1). Previous searches for high-mass diphoton resonances using the [2015][2016] collision data have been reported by the ATLAS and CMS collaborations [11,12]. Modifications with respect to the analysis in Ref.
[11] include a common event selection for the spin-0 and spin-2 resonance searches, the use of the functional decomposition method [13] to assess the spurious-signal uncertainty, and updates to the photon reconstruction, identification, isolation and energy calibration.
This search uses functional forms to describe the signal and background components in a fit of the diphoton invariant mass spectrum, , to determine the signal yield. In the absence of a significant signal excess, limits are placed on the production cross-section times branching ratio, × ( → ). For the spin-0 case, the limits are computed in terms of the fiducial cross-section, defined as the product of the cross-section times the branching ratio to two photons within a fiducial acceptance which closely follows the selection criteria applied to the reconstructed data. In the spin-2 analysis, limits are instead placed on the total cross-section times branching ratio to the two-photon final state. These limits are provided as a function of resonance mass and width in the spin-0 search and as a function of graviton mass and coupling / Pl for the spin-2 resonance search. As the graviton natural width is related to the coupling via Γ * = 1.44( / Pl ) 2 * , the graviton is expected to be a narrow resonance [14] for the couplings considered in this search, i.e. / Pl < 0.1.

ATLAS detector
The ATLAS experiment [15][16][17] 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 surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadronic 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 hadron sampling calorimeter made of steel and scintillator tiles covers the central pseudorapidity range (| | < 1.7). The endcap and forward regions are instrumented with LAr calorimeters 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 for EM and hadronic energy measurements up to | | = 4.9. The muon spectrometer surrounds the calorimeters and is based on three large air-core toroidal superconducting 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. 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 keep the accepted rate below 100 kHz. A software-based trigger reduces the accepted event rate to an average of 1 kHz depending on the data-taking conditions.

Data and simulated event samples
The search is carried out using the √ = 13 TeV collision dataset with a bunch spacing of 25 ns collected from 2015 to 2018, with stable beam conditions and all ATLAS subsystems operational, which corresponds to an integrated luminosity of 139.0 ± 2.4 fb −1 [18]. The data were recorded using a diphoton trigger that required two electromagnetic clusters with transverse energies T above 35 and 25 GeV, respectively, both fulfilling photon identification criteria based on shower shapes in the electromagnetic calorimeter. The efficiency of the diphoton trigger with respect to the full event selection is above 99% for 2015-2016 and above 98% for 2017-2018 [19].
Simulated Monte Carlo (MC) events are used to optimize the analysis selections and characterize the signal and background. Interference effects between the resonant signal and all background processes are expected to be small for a narrow-width signal and are neglected in this analysis to keep it as model-independent as possible.
To model the spin-0 resonance signals, MC samples were generated for a hypothetical resonance mass in the range 200-3000 GeV assuming a decay width Γ of 4 MeV to describe a hypothetical resonance in the narrow-width approximation (NWA). Four production modes were generated to assess the model dependence of the fiducial selection. The gluon-gluon fusion (ggF) samples were generated at next-toleading order (NLO) in pQCD using M G 5_ MC@NLO [20], with a set-up matching that described in Ref. [11]. Vector-boson fusion (VBF) samples were generated at NLO in pQCD with P B v2 [21-23] interfaced with P 8 [24] using the CT10 parton distribution function (PDF) set [25] and the AZNLO tune [26] of the parameter values. Simulated samples for production of a scalar resonance in association with a vector boson or¯pair were generated at leading order (LO) in pQCD with P 8 using the NNPDF23LO PDF set [27] and the A14 tune [28]. Additional gluon-gluon fusion samples were produced for various fixed decay widths, ranging from 4 MeV up to 10% of .
To model the spin-2 resonance signals, MC samples incorporating the RS1 model were generated at LO in pQCD using P 8 with the NNPDF23LO PDF set and the A14 tune, with masses * in the range 500-5000 GeV and fixed coupling values / Pl ranging from 0.01 to 0.1.
The irreducible background in this search, coming from events with two prompt photons, was simulated using the S [29, 30] event generator, version 2.2.4. Matrix elements were calculated with up to one additional parton at NLO and with two or three partons at LO in pQCD and merged with the S parton-shower simulation using the ME+PS@NLO prescription [31][32][33][34]. The NNPDF3.0 NNLO PDF set [35] was used in conjunction with a dedicated parton-shower tune in the S generator.
The effects of multiple interactions in the same bunch crossing as the hard scatter and in neighbouring ones are included using simulated events generated with P 8. Simulated events were weighted to reproduce the distribution of the average number of interactions per bunch crossing observed in data.
All simulated signal events were processed using a full simulation of the ATLAS detector [36] based on G 4 [37]. The background events were processed using a fast simulation of the ATLAS detector [38], where the full simulation of the calorimeter is replaced with a parameterization of the calorimeter response. All simulated events were reconstructed with the same reconstruction algorithms as those used for data.

Object and event selection
Photon candidates are reconstructed from topological clusters of energy deposits in the EM calorimeter and calibrated as described in Ref. [39]. The event selection requires at least two photon candidates with T > 25 GeV and | | < 2.37, excluding the barrel-to-endcap transition regions of the calorimeter, 1.37 < | | < 1.52. The two highest-T photons and additional information from the tracking systems are used to identify the diphoton production vertex [40].
To reduce the background from jets, photon candidates are required to fulfil tight identification criteria based on shower shapes in the EM calorimeter and energy leakage into the hadronic calorimeter [39]. A looser identification is considered for background estimations that benefit from a larger sample size. The tight identification is optimized in sub-ranges of photon T and | | to have an identification efficiency greater than 90% for the T range considered in this analysis. The calorimeter isolation transverse energy iso T is required to be smaller than 0.022 T + 2.45 GeV, where iso T is defined as the sum of transverse energies of the positive-energy topological clusters [41] within a cone of size Δ = 0.4 around the photon candidate. The core of the photon shower is excluded, and iso T is corrected for the leakage of the photon shower into the isolation cone. The contributions from the underlying event and pile-up are subtracted using the techniques described in Refs. [42,43]. The track isolation is defined as the scalar sum of the transverse momenta of tracks with T > 1 GeV which satisfy some loose track-quality criteria, are not associated with a photon conversion and originate from the diphoton production vertex in a Δ = 0.2 cone around the photon candidate. It is required to be less than 0.05 T .
The diphoton invariant mass is evaluated using the energies of the leading-and subleading-T photons, their separation Δ in azimuthal angle and Δ in pseudorapidity determined from their positions in the calorimeter, and the position of the diphoton production vertex. Additional kinematic requirements are placed on the photon T relative to the diphoton invariant mass: the leading photon must have T / > 0.3 and the subleading photon must satisfy T / > 0.25. These requirements were optimized to retain improved significance with respect to Ref.
[11] while allowing the analysis to have a unified selection for the spin-0 and spin-2 models. For the spin-0 model, the new selections result in a slightly worse expected limit for the NWA scalar models below 700 GeV but up to 30% better expected limits above. For the spin-2 model, these updated selections lead to an improvement in the expected limit by 15% to 50% depending on * for masses up to 2 TeV.
In total, 433 919 data events with > 150 GeV are selected.
The fiducial volume for the spin-0 interpretation is defined by requiring two photons at generator level with | | < 2.37, and T / > 0.3 and T / > 0.25 for the leading and subleading photons, respectively. The particle isolation, defined as the scalar sum of T of all the stable particles (except neutrinos) found within a Δ = 0.4 cone around the photon direction, is required to be less than 0.05 T + 6 GeV.

Signal modelling
Parametric models of the diphoton invariant mass distributions of the signal are used to test for different resonance masses and Γ or / Pl . The detector resolution is modelled by a double-sided Crystal Ball (DSCB) function [44,45], with parameters expressed as a function of . The DSCB function is convolved with the true lineshape, which is parameterized as the product of a relativistic Breit-Wigner (BW) function and mass-dependent factors accounting for the parton luminosity and the matrix elements of the processes. This convolution is performed independently for the spin-0 and spin-2 models. The lineshapes for the spin-0 model are derived from the functional forms defined in the M G 5_ MC@NLO MC simulation of the ggF process, where a BW function of width Γ is combined with a mass-dependent gluon-gluon luminosity functional form. The lineshapes for the spin-2 model are derived from the functional forms defined in the P 8 MC simulation of the RS1 model [46], where a BW function of width 1.44( / Pl ) 2 * is combined with mass-dependent gluon-gluon and quark-antiquark luminosity functional forms.

Background estimates
The unified spin-0 and spin-2 search selections in this analysis allow the use of a common background modelling procedure. The largest background component comes from the non-resonant production of photon pairs ( events); smaller backgrounds come from events containing a photon and a jet ( events) and events with two jets ( events), where the jets are misidentified as photons. The relative contribution of these processes is determined using the two-dimensional sideband method described in Ref.
[47] and shown in Figure 1. The overall purity of the selected events increases with from around 89% at 150-200 GeV to around 97% above 400 GeV with an uncertainty of 0.5-3%. No significant difference in purity is seen between the LHC data-taking periods, and the remaining background is dominated by events.
The analytic form of the continuum background is determined by fitting a template built from reconstructed S MC events and events derived from a dedicated data control region where the photon identification requirements are inverted. The diphoton invariant mass distribution is fitted in the range above 150 GeV. The lower value of the range is chosen to allow for enough events in the side-bands to ensure a good description by the analytic form, while the upper value of the range is chosen such as the procedure remains stable in injection tests . The search region for a resonant signal covers the region 160-3000 GeV for a NWA spin-0 resonance, 400-2800 GeV for a wide spin-0 resonance, and 500-2800 GeV for a graviton resonance. The procedure described in Ref.
[40] is used to ensure that the chosen analytic function is flexible enough to model any potential variations in the background template. The analytic form is chosen from the family of functions previously used to describe the diphoton invariant mass spectrum [48]. Additional background templates are constructed from variations due to the uncertainties in the measured background composition, shape of the component and theoretical uncertainties such as the choice of PDF and the variation of renormalization and factorization scales. The function with the fewest degrees of freedom that maintains enough flexibility to model all variations of the template shape is chosen, and it takes the form where , , 0 , 1 are free parameters and = / √ .  Fits to these background templates are used to estimate the bias due to this choice of analytic function; any fitted signal yield is considered as a 'spurious-signal' systematic uncertainty. Up to masses of 2 TeV this uncertainty is dominated by the statistical fluctuations of the background template due to the limited number of generated MC events, leading to fluctuations that are significant compared to the expected statistical uncertainty on the data. To suppress the impact of these fluctuations, each background template is smoothed using the functional decomposition (FD) method [13]. This method uses a linear combination of orthonormal exponential functions to fit a background template and the fit result is binned and used as a smoothed template for the spurious-signal determination only.
To illustrate the FD process, an example fit to a pseudo-experiment and the resulting spurious-signal uncertainty is shown in Figure 2. This pseudo-experiment dataset is generated from the functional form given in Eq. (1), where the parameters are determined by fitting simulated events. First, the FD method is used to fit the pseudo-experiment, resulting in a smoothed FD template; this template and the pseudo-experiment are shown in Fig. 2(a). The spurious-signal estimated by fitting the smoothed FD template with the signal-plus-background model has fewer fluctuations and a smaller amplitude than the one estimated by fitting the unsmoothed template as shown in Fig. 2(b): it can be larger than 50% of the statistical uncertainty with the unsmoothed template while it decreases from 40% at 160 GeV to 20% at 500 GeV and less than 10% above 1 TeV with the smoothed template. To verify that the smoothing maintains the underlying shape of the distribution, this process is repeated for an ensemble of pseudo-experiments and the mean number of extracted signal events is compared between the smoothed and unsmoothed templates. The resulting bias in the determination of the spurious-signal uncertainty, shown in Fig. 2(c), agrees between the unsmoothed and smoothed template ensembles and is found to be much smaller than the uncertainty from a single pseudo-experiment. It should be noted that the mean spurious-signal is not exactly equal to 0 because of the limited number of pseudo-experiments. Finally, by using this smoothed template approach, an improved description of the spurious-signal uncertainty is obtained, leading to an improvement of 2-28% in the expected sensitivity of this search as shown in Fig. 2(d). (d) Figure 2: An illustration of the smoothing method used to determine the spurious-signal systematic uncertainty. (a) An example FD fit to the spectrum of a pseudo-experiment dataset generated from the functional form fit to simulated events. (b) Number of spurious-signal events relative to the statistical uncertainty of the input dataset calculated using signal-plus-background fits to the input dataset (black) and the smoothed FD template (red). The background function used in these fits is the same function used to generate the pseudo-experiment dataset, and the variance of the unsmoothed template highlights the role of the number of events in the MC sample in the systematic uncertainty determination. (c) Mean of the relative spurious-signal calculated by repeating the signal-plus-background fits on a large ensemble of pseudo-experiments and the FD smoothed templates. The shaded area indicates the size of the spurious-signal uncertainty from the single smoothed template pseudo-experiment of Fig. 2(b). (d) The effect of the reduced systematic uncertainty on the expected fiducial cross-section limits for the NWA spin-0 signal hypothesis. The use of smoothed templates leads to a limit improvement of 2-28%.
The spurious-signal uncertainty for this background function choice is 40-10% of the statistical uncertainty of the fitted signal yield for NWA signals with masses ranging from 160 GeV to 3000 GeV. This uncertainty increases to 70-20% for wider signals with masses ranging from 400 GeV to 2800 GeV.

Statistical procedure
The statistical analysis of the data uses binned maximum-likelihood fits of the distribution. The systematic uncertainties listed in Table 1 are accounted for in the fits by using nuisance parameters constrained by Gaussian penalty terms in the likelihood function. They include the experimental uncertainties in the luminosity determination, pile-up profile, trigger efficiency, and photon identification and reconstruction. Additional systematic uncertainties in the signal shape come from the uncertainties in the photon energy resolution and scale. For the spin-0 model, the correction factors inside the fiducial volume, used to compute the limit on the fiducial cross-section, are determined using gluon-gluon fusion MC events. They are compared with the ones computed for the other production modes, leading to a small mass-dependent systematic uncertainty in the signal yield. In the absence of a signal, the expected and observed 95% confidence level (CL) exclusion limits on the cross-section times branching ratio to two photons are computed using a modified frequentist approach CL s [50,51], assuming the asymptotic approximation for the test-statistic distribution. In regions where very few events are observed and the expected signal yield is also small, the asymptotic approximation is found to deviate from the limits based on pseudo-experiments. To account for this effect, for masses greater than 1000 GeV, pseudo-experiments are used to verify the expected and observed limits, and used in place of the asymptotic limit where differences are observed.
For the spin-2 model, the limits on the total production cross-section can be translated into a lower limit on the graviton mass. For that purpose the spin-2 exclusion limits are extended up to * = 5000 GeV for all coupling hypotheses.

Results
The diphoton invariant mass distribution of the events passing the selection is shown in Figure 3, along with the background-only fit. The highest-mass diphoton candidate selected in data has a mass of 2.36 TeV.
The probability of compatibility with the background-only hypothesis, quantified with the local -value expressed in standard deviations, is shown in Figure 4 as a function of the hypothesized resonance mass and width. The most significant excess is observed for a mass of 684 GeV for the spin-0 NWA model and for the spin-2 / Pl = 0.01 model, corresponding to a 3.29 local significance. The corresponding global significance is 1.30 and 1.36 for the two models, respectively.
The observed and expected limits, shown in Figure 5 for two values of the signal width, are in good agreement, consistent with the absence of a signal.
The observed limits for different masses and different values of Γ / or / Pl are summarized in Table 2. For all signals considered in this analysis, the largest single improvement in the limits with respect Table 1: Summary of the systematic uncertainties. The mass-dependent values are highlighted with a star * and are given for signal masses between 160 and 2800 GeV unless stated otherwise.

Signal yield
Luminosity ±1.7% Trigger ±0.5% Photon identification ±0.5% Photon isolation ±1.5% Photon energy scale/resolution negligible Pile-up reweighting * ±(2-0.2)% Spin-0 production process * ±(7-3)% Signal modelling * Photon energy resolution  with several generic narrow-width signal shapes overlaid (dotted lines), scaled to 10 times the value of the corresponding expected upper limit at 95% CL on the fiducial cross-section times branching ratio, with pole masses of = 0.4, 1 and 2 TeV. The normalized residuals between the data and the fit are shown in the bottom panel with denoting only the statistical error of the data. There is no data event with > 2.36 TeV. Table 2: Summary of the observed upper limits on the fiducial and the total production cross-section times branching ratio to two photons for the spin-0 and spin-2 models, respectively.   Figure 5: The expected and observed upper limits at 95% CL on (a) the fiducial and (b) the total production cross-section times branching ratio to two photons of (a) a NWA spin-0 resonance as a function of its mass and (b) the lightest KK graviton as a function of its mass * for / Pl = 0.1. For masses greater than 1000 GeV, pseudo-experiments are used to verify the expected and observed limits, and used in place of the asymptotic limit where differences are observed. Cross-section predictions for the RS1 model are shown in (b), where the grey band illustrates the PDF uncertainty.

Conclusion
Searches for new resonances in high-mass diphoton final states with the ATLAS experiment at the LHC are presented. The searches use the full LHC Run 2 proton-proton collision dataset at a centre-of-mass energy of √ = 13 TeV, corresponding to an integrated luminosity of 139 fb −1 . The analyses are optimized to search for spin-0 resonances with masses above 200 GeV and for spin-2 resonances predicted by the Randall-Sundrum model with masses above 500 GeV.
The data are consistent with the Standard Model background expectation. In the spin-0 resonance search, the observed 95% CL upper limits on the fiducial cross-section times branching ratio for a narrow-width signal range from 12.5 fb at 160 GeV to about 0.03 fb at 2800 GeV. In the spin-2 resonance search, the observed limits on the total cross-section times branching ratio for / Pl = 0.1 range from 3.2 fb to about 0.04 fb for a graviton mass between 500 and 5000 GeV, and the RS1 graviton is excluded for masses below 4500 GeV for this coupling hypothesis.