Search for pairs of highly collimated photon-jets in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector

Results of a search for the pair production of photon-jets$-$collimated groupings of photons$-$in the ATLAS detector at the Large Hadron Collider are reported. Highly collimated photon-jets can arise from the decay of new, highly boosted particles that can decay to multiple photons collimated enough to be identified in the electromagnetic calorimeter as a single, photonlike energy cluster. Data from proton-proton collisions at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 36.7 fb$^{-1}$, were collected in 2015 and 2016. Candidate photon-jet pair production events are selected from those containing two reconstructed photons using a set of identification criteria much less stringent than that typically used for the selection of photons, with additional criteria applied to provide improved sensitivity to photon-jets. Narrow excesses in the reconstructed diphoton mass spectra are searched for. The observed mass spectra are consistent with the Standard Model background expectation. The results are interpreted in the context of a model containing a new, high-mass scalar particle with narrow width, $X$, that decays into pairs of photon-jets via new, light particles, $a$. Upper limits are placed on the cross section times the product of branching ratios $\sigma \times \mathcal{B}(X \rightarrow aa) \times \mathcal {B}(a \rightarrow \gamma \gamma)^{2}$ for 200 GeV $<m_{X}<$ 2 TeV and for ranges of $ m_a $ from a lower mass of 100 MeV up to between 2 and 10 GeV, depending upon $ m_X $. Upper limits are also placed on $\sigma \times \mathcal{B}(X \rightarrow aa) \times \mathcal {B}(a \rightarrow 3\pi^{0})^{2}$ for the same range of $ m_X $ and for ranges of $ m_a $ from a lower mass of 500 MeV up to between 2 and 10 GeV.


Introduction
The quest for new particles at the Large Hadron Collider (LHC) at CERN has been greatly rewarded by closely examining collision events that contain photons in the final state. Despite the relatively small branching ratio predicted for the process in the Standard Model (SM), the decay of the Higgs boson into two photons is readily identifiable due to the good energy resolution of the electromagnetic (EM) calorimeters of the ATLAS [1] and CMS [2] detectors and the relatively small backgrounds in final states with only photons. The search for this process was one of the main methods by which the Higgs boson was observed [3,4]. Moreover, the establishment of a wide range of results that so far are consistent with the SM at the LHC at a center-of-mass energy of 13 TeV motivates a renewed focus on searches for new physics that target general experimental signatures, including non-standard photon signatures, rather than specific signal models. In many beyond the Standard Model (BSM) theories [5][6][7][8][9][10][11][12][13], new scalar, pseudoscalar or vector gauge bosons can decay into photon-only final states that lead to collimated groupings of photons ("photon-jets" [14,15]). In some cases, the Lorentz boost of the new particles is large enough to lead to an opening angle between the trajectories of the final-state photons that is smaller than or comparable to the angular size of an energy cluster in the EM calorimeter corresponding to a single photon, resulting in highly collimated photon-jets. Such boosted particles arise, for example, when a high-mass particle produced in the proton-proton collision decays into intermediate particles, with much lower masses, that subsequently decay into photons. Thus, events selected to contain two, well-separated, reconstructed photons can be used to search for pairs of highly collimated photon-jets resulting from BSM particle decays.
A search for highly collimated photon-jets using 36.7 fb −1 of LHC proton-proton collision data collected by the ATLAS detector in 2015 and 2016 at a center-of-mass energy of 13 TeV is presented. Candidate photon-jet pair production events are selected from those containing two reconstructed photons (denoted "γ R "), using a set of identification criteria much less stringent than that typically used for the selection of photons, and with additional criteria applied to provide improved sensitivity to photon-jets. Narrow excesses are searched for in the spectra of the reconstructed diphoton mass m γ R γ R .
The results of the search are interpreted in the context of a benchmark BSM scenario involving a highmass, narrow-width scalar particle, X, with mass m X > 200 GeV, originating from the gluon-gluon fusion process and that can decay into a pair of intermediate particles with spin 0, a, as shown in Figure 1. The a particle can in general decay to several final states, but here is restricted to decay either into a pair of photons, via X → aa → 4γ, or into three neutral pions, via X → aa → 6π 0 → 12γ, yielding events containing a pair of photon-jets of either low or high multiplicity; the result is interpreted for both cases. Because the search is performed using events that contain two calorimeter deposits that are initially loosely identified as individual photons, the search is sensitive to the parameter region in which the a particle is highly boosted, m a < 0.01 × m X .

ATLAS detector
The ATLAS detector [1] is a multipurpose detector with a forward-backward symmetric cylindrical geometry. 1 The detector covers nearly the entire solid angle around the collision point. It consists of an inner tracking detector surrounded by a thin superconducting solenoid, EM and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroid magnets.
The inner-detector system is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range |η| < 2.5. The high-granularity silicon pixel detector covers the vertex region. The innermost layer of the pixel detector, the insertable B-layer [16], was installed between Run 1 and Run 2 of the LHC. The pixel detector typically provides four measurements per track. It is followed by the silicon microstrip tracker that normally provides four two-dimensional measurement points per track. These silicon detectors are complemented by the transition radiation tracker, which enables radially extended track reconstruction up to |η| = 2.0. The transition radiation tracker also provides electron identification information based on the fraction of hits (typically 30 in total) above a higher energy-deposit threshold corresponding to transition radiation.
For |η| < 2.5, the EM calorimeter is composed of three sampling layers in the longitudinal direction of shower depth. The first layer is segmented into high-granularity strips in the η direction, with a typical cell size of ∆η × ∆φ = 0.003 × 0.1 for the ranges |η| < 1.4 and 1.5 < |η| < 2.4, and a coarser cell size of ∆η × ∆φ = 0.025 × 0.1 for other regions. This fine granularity in the η direction allows identification of events with two overlapping showers originating from the decays of neutral hadrons in hadronic jets, mostly π 0 → γγ decays. The second layer has a cell size of ∆η × ∆φ = 0.025 × 0.025. This layer collects most of the energy deposited in the calorimeter by photon and electron showers. The third layer is used to correct for energy leakage beyond the EM calorimeter from high-energy showers. The thicknesses of the first, second, and third layers at η = 0 are 4.3 radiation lengths (X 0 ), 16 X 0 , and 2 X 0 , respectively, and they vary with the pseudorapidity range [1]. Placed in front of these layers, an additional thin LAr presampler layer covering |η| < 1.8 is used to correct for energy loss in material upstream of the 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point in the center of the detector and the z-axis along the beam pipe. The x-axis points from the interaction point to the center of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of ∆R ≡ (∆η) 2 + (∆φ) 2 .
calorimeters. Hadronic calorimetry is provided by the steel/scintillator-tile calorimeter, segmented into three barrel structures within |η| < 1.7, and two copper/LAr hadronic endcap calorimeters. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimized for EM and hadronic measurements respectively.
A two-level trigger system, the first level implemented in custom hardware followed by a software-based level, is used to reduce the event rate to about 1 kHz for offline storage.

Photon-jet signal characteristics
Photon-jets, defined as collimated groupings of photons, can arise from decays of particles that are highly boosted as a result of themselves being the decay products of higher-mass particles. For the benchmark BSM scenario considered here, the extent to which photons from decays of a particles are collimated depends on the ratio of the masses of the X and a, particularly in the case where the X particle is produced with a momentum significantly less than its mass.
For large values of the ratio m a /m X , the boost of the a is small enough to yield more than two individual photons, well separated and isolated, that can be identified in the detector. In this regime, a general search for new phenomena in events with at least three isolated photons, using a three-photon trigger, was performed by ATLAS at 8 TeV [17]. This search was sensitive to cases where the angular separation between photons was large, for ∆R γγ 0.3, which corresponds to m a /m X 0.08 for the benchmark signal scenario. For slightly smaller values of the ratio m a /m X , the individual final-state photons appear too close together in the detector and fail isolation criteria, limiting the sensitivity of the 8 TeV ATLAS search in this regime.
For very small values of the ratio m a /m X , the boost of the a is large enough to lead to angular separations between the final-state photons of ∆R γγ 0.04, which is approximately the same size as a standard singlephoton energy cluster in the ATLAS EM calorimeter. In this case, existing triggers cannot distinguish a calorimeter energy deposit resulting from highly collimated photons from that of a single photon. Thus, diphoton-like events can be used as a starting point for a search for highly collimated photon-jets, and the sensitivity to this region of the photon-jet parameter space can be increased by placing criteria on the shape of the shower in the EM calorimeter in addition to those applied in the trigger. This analysis presents a search for highly collimated photon-jets that is sensitive to a wide mass range for the parent X particle, m X > 200 GeV, and for m a /m X < 0.01 in the benchmark signal scenario.
For this benchmark scenario, for the process X → aa → 4γ, the distribution of ∆R γγ is shown in Figure 2. Due to the kinematics of boosted particles, ∆R γγ has a maximum at a value of 2/γ, where γ is the Lorentz factor of the a particle, γ = E a /m a . When the X particle is produced nearly at rest, since the energy of the a particle has a median value of E a ∼ m X /2, the distribution of ∆R γγ has a maximum at ∼ 4 × m a /m X . The approximate proportionality of the angular spread of photons within the photon-jet to m a /m X holds for photon-jets in general, including those with larger photon multiplicity resulting from processes such as X → aa → 6π 0 . Since the two different final states of the benchmark scenario are similar, some parts of the descriptions in the following sections are only mentioned for the X → aa → 4γ decay to avoid repetition, although they apply to the X → aa → 6π 0 decay as well.
For values of the ratio m a /m X greater than 0.01, the final-state photons are separated enough to lead to a relatively large cluster of energy in the calorimeter, and such events do not satisfy the isolation criteria or the initial loose identification of photons at the trigger level. The signal selection efficiency for the present  Figure 2: The distribution of ∆R γγ , the angular separation between two photons that are reconstructed as a single photon-jet in the ATLAS detector, for the benchmark signal scenario for the process X → aa → 4γ, using simulated signal samples at generator level. The distribution has a peak at ∼ 4 × m a /m X and a long tail on the right side. For the values of m X and m a presented in the figure, ∆R γγ is smaller than or comparable in size to an EM cluster. analysis in this m a /m X > 0.01 region is lower than 4%, and so no attempt is made to search in this regime. There is therefore an intermediate region, 0.01 < m a /m X 0.08, which is covered by neither this search nor the previous search for three-photon final states at 8 TeV.

Event samples
The data sample used for this search corresponds to an integrated luminosity of 36.7 fb −1 (after applying data-quality requirements), collected under normal data-taking conditions for pp collisions during 2015 and 2016 at a center-of-mass energy of √ s = 13 TeV. The data were selected using an unprescaled trigger that filters events with two energy deposits in the EM calorimeter that satisfy trigger-level loose photon identification criteria with transverse energy values of E T,1 > 35 GeV and E T,2 > 25 GeV.
Samples of the benchmark signal scenario with two different final states, X → aa → 4γ and X → aa → 6π 0 , were simulated. For the production of the X via gluon-gluon fusion, M G 5_aMC@NLO [18] Version 2.3.3, at next-to-leading order (NLO) in quantum chromodynamics (QCD) with the NNPDF30NLO parton distribution function (PDF) set [19], was used. For the subsequent decay of the X into aa and into the photon-jet final states, P 8 [20] Version 8.210, with the A14 set of tuned parameters [21], was used, as well as for the parton-shower and hadronization simulation of initial state radiation jets. The samples were produced using a narrow-width approximation (NWA) approach with the resonance widths of the X and a set to 4 MeV and 1 MeV, respectively. Samples were simulated for mass ranges of 200 GeV < m X < 2000 GeV and 0.1 GeV < m a < 0.01 × m X .
The non-resonant production of diphoton events in the SM is the dominant background source for this analysis, and these events were simulated with S 2.1.1 [22]. Matrix elements were calculated with up to two additional partons at leading order (LO) in QCD and merged with the S partonshower simulation [23] using the ME+PS@LO prescription [24]. The CT10 PDF set [25] was used in conjunction with a dedicated parton-shower tune of S . These samples are used to validate the background modeling based on analytic functions (described in Section 6.2). Simulated samples of the reducible SM background consisting of one photon and one hadronic jet from the hard process were also generated with S 2.1.1-using the same PDF set, parton-shower tune, and merging prescription as for the diphoton sample-with matrix elements calculated at LO with up to four additional partons. These samples are used for optimizing the search strategy described in Section 5.
Additional interactions in the same or nearby bunch crossings (pileup) were simulated using Pythia 8.186 [20] using the A2 set of tuned parameters [26] and the MSTW2008LO PDF [27] set and overlayed on the simulated signal and SM background events. All simulated event samples were produced using the ATLAS simulation infrastructure [28], using the full G 4 [29] simulation of the ATLAS detector. Simulated events were then reconstructed with the same software as used for the data.

Object and event selection
This analysis selects events containing at least two reconstructed photons, obtained from a diphoton trigger, and then searches for pair-produced photon-jets. This is accomplished by applying additional selection criteria and scanning for deviations from the expected background in the m γ R γ R spectrum, defined as the distribution of the mass values of the two reconstructed photons, which would correspond to the mass of the high-mass particle m X in the case of a signal event. No attempt is made to reconstruct the mass of the a in the process X → aa → photon-jets (although specifics of the a are taken into account in several parameters of the signal modeling, which is detailed in Section 6.1).

Initial event selection with two loose photons
Reconstructed photons are obtained from clusters of energy deposited in the EM calorimeter [30]. In the barrel section a cluster size of 3 × 7 cells in the middle layer is used (equivalent to an area of size ∆η × ∆φ = 0.075 × 0.175), while a cluster of 5 × 5 cells in the middle layer is used in the endcap sections (equivalent to an area of ∆η × ∆φ = 0.125 × 0.125). Reconstructed photons are required to match photon objects calculated at the trigger level, within the separation of ∆R < 0.07, and may have associated tracks and conversion vertices reconstructed in the inner detector.
The two leading reconstructed photons are required to be within the fiducial calorimeter region of |η| < 2.37, excluding the transition region at 1.37 < |η| < 1.52 between the barrel and endcap calorimeters. The criterion E T,1 > 0.4 × m γ R γ R is applied to the leading reconstructed photon, and E T,2 > 0.3 × m γ R γ R to the subleading reconstructed photon. These criteria increase the sensitivity to photon-jet pairs from a scalar resonance, since such candidate signal events tend to contain photons with larger E T /m γ R γ R ratios compared with those from background events dominated by t-channel processes [31]. Only events with m γ R γ R > 175 GeV are selected for further analysis.
The two leading reconstructed photons are required to be isolated from other calorimeter energy deposits and from nearby tracks not associated with the photon. The calorimeter isolation variable E iso T is defined as the sum of energy deposits in the calorimeter in a cone of size ∆R = 0.4 around the barycenter of the photon cluster (excluding the energy associated with the photon cluster) minus 0.022 × E T . This cone energy is corrected for the leakage of the photon energy from the photon cluster and for the effects of pileup [32]. The calorimeter isolation variable is required to satisfy E iso T < 2.45 GeV. The track isolation variable p iso T is defined as the scalar sum of the transverse momenta of tracks not associated with the photon in a cone of size ∆R = 0.2 around the barycenter of the photon cluster. It is required to satisfy p iso T < 0.05 × E T .

Optimized photon selection for photon-jet signatures
Photon identification in ATLAS [30] is based on a set of requirements placed on several discriminating variables that characterize the shower development in the calorimeter ("shower shapes"), defined to reject the background from hadronic jets misidentified as photons. Nine discriminating variables are defined, and they are described in detail in Table 1 of Ref. [30]. One variable quantifies the shower leakage fraction in the hadronic calorimeter, and three variables quantify the lateral shower development in the EM calorimeter second layer. The other five variables quantify the lateral shower development in the finely segmented strips of the first layer, and two of them are utilized to identify photon candidates with two separate local energy maxima in the fine strips, which are characteristic of neutral hadron decays in hadronic jets, primarily from π 0 → γγ.
Several reference selections are defined, including those referred to as "Loose" and "Tight". The Loose selection is based only on shower shapes in the second layer of the EM calorimeter and on the leakage in the hadronic calorimeter, and is used by the photon triggers, including the diphoton trigger used for the collection of the data sample for this search. The Tight selection is based on all nine variables and is used for the standard photon identification in ATLAS, but is not used in this search. The criteria for the Tight selection change as a function of the η values of the reconstructed photons, to account for the calorimeter geometry and effects from the material upstream of the calorimeter, and are separately optimized for reconstructed photons with and without an associated conversion vertex to increase the photon identification efficiency.
In this search, both reconstructed photons are required to fulfill the "Loose " selection. This selection is defined by removing requirements on all five variables quantifying the shower development in the finely segmented strip layer of the calorimeter (w s 3 , w s tot , F side , ∆E, and E ratio , defined in Table 1 of Ref. [30]), with respect to the standard Tight selection. The requirements on the other four variables (R had , R η , w η2 , and R φ ) remain the same as for the standard Tight selection. By definition, the Loose is an intermediate selection between Loose and Tight. Based on simulated samples of signal and SM background processes, this Loose selection provides good sensitivity to photon-jet signals. This is explained by the fact that energy clusters of photon-jets exhibit multiple local energy maxima in the fine strip layer, since the angular separation of photons constituting the photon-jet can be larger than the segmentation of the strips, depending on the mass parameters m X and m a of the benchmark signal scenario, as seen in Figure 2.
For signal mass values 0.003 < m a /m X < 0.006 and m X > 200 GeV, the total selection efficiency is less than 5% when the standard Tight selection is applied, in addition to the selection criteria described in Section 5.1, and this increases to 20%-50% with the Loose selection. Comparing the two selection criteria, an increase in the overall event yield of roughly 30% is observed with the Loose selection. Thus, the analysis sensitivity to photon-jet signals is increased by the use of the Loose selection, rather than the standard Tight selection.
Additionally, the choice of Loose allows the definition of a set of "not Loose " criteria (i.e., where at least one of the two reconstructed photons fails the Loose selection) that is used to define the control regions for the evaluation of the background composition, as described in Section 6.2.

Categorization of events by the shower shape variable ∆E
After the preselection of events with two leading reconstructed photons satisfying the isolation and Loose identification criteria described in the previous sections, the final signal region is defined by dividing the events into two orthogonal categories based on the value of the calorimeter variable ∆E for the reconstructed photons. The quantity ∆E corresponds to a shower shape variable based on information in the first layer of the EM calorimeter, and quantifies the relative size of multiple, individual energy deposits that may be contained within a single energy cluster.
It is defined as min where E S1 2nd max is the energy of the strip cell with the second-largest energy, and E S1 min is the energy in the strip cell with the lowest energy located between the strips with the largest and the second-largest energy. If the strip cells with the largest and the second-largest energy are located next to each other, or if there is no second-largest energy strip, then ∆E = 0. This variable is useful for identifying the π 0 → γγ process, prevalent in hadronic jets, which leaves a characteristic signature in the first layer of the EM calorimeter that often yields two peaks in the η direction, resulting in large ∆E values. When the photon-jet signals from decays such as a → γγ and a → 3π 0 → 6γ have angular separation of photons larger than the segmentation of the first layer of the EM calorimeter, they leave signatures in the calorimeter similar to π 0 → γγ events. Thus, the variable ∆E is used to effectively select photon-jet signals.
The categorization by ∆E is as follows: • low-∆E category: Both reconstructed photons are required to have values of ∆E below given thresholds. This requirement corresponds to reconstructed photons with a signature in the fine strip layer similar to that of single photons.
• high-∆E category: At least one of the two leading reconstructed photons is required to have a value of ∆E above a given threshold. This requirement corresponds to events containing reconstructed photons which have a π 0 -like signature.
The thresholds for the value of ∆E used to determine whether an event appears in either the low-or high-∆E category are the same as those used in the standard Tight photon selection. These thresholds range from 100 MeV to 500 MeV, depending on the photon η and whether there are associated tracks or conversion vertices.
The high-∆E category is found to have a significantly better signal-to-background ratio compared with the low-∆E category, since reconstructed photons with large ∆E values typically correspond to photonjets with a larger angular spread among the constituent photons. The high-∆E criterion also effectively reduces the contribution of single photons, which tend to have small ∆E values. Hadronic jets from SM processes containing π 0 → γγ decays are likely to fall into the high-∆E category, but the contribution of these events is small due to the isolation requirements. This leads to lower expected background event yields in the high-∆E category, resulting in better signal-to-background ratios compared with the low-∆E category. The number of events observed in each category for different ranges of m γ R γ R is shown in Table 1. Although overall the ratio of signal-to-background is lower for the low-∆E category, it still provides increased sensitivity to photon-jet signals with smaller angular separation, and so both categories are used in this search.

Summary of the selection
The overall efficiency, ε, of selecting signal events after applying all criteria, including kinematic acceptance and excluding the categorization by ∆E, is shown in Figure 3 (a), and the fraction, f , of signal events that appear in the low-∆E category is shown in Figure 3 (b), both as a function of m a and m X for the decay X → aa → 4γ. The selection efficiency is low for small values of m X and large values of m a , and almost all events are in the low-∆E category for large values of m X and small values of m a . For smaller m a and larger m X , f increases because of the small angular spread of photons inside the photon-jet, which leads to a calorimeter signature similar to that of a single photon. Additionally, for larger m a and smaller m X , f also increases because individual photons are reconstructed separately due to the large angular separation, resulting in events containing more than two reconstructed photons, each of which more resembles a single photon. The results for the decay X → aa → 6π 0 are similar to those of the decay X → aa → 4γ. [GeV]  Table 2 displays the number of events in data that satisfy each selection criterion. The fraction of events with both of the two leading reconstructed photons found in |η| < 1.37 (i.e. the barrel section) is 59% (63%) for the low-∆E (high-∆E) category.

Signal and background modeling
The reconstructed signal mass shape is modeled with a double-sided Crystal Ball (DSCB) function. The backgrounds are determined by fitting functions to the observed mass spectra of two reconstructed photons, Table 2: Number of events in collision data that satisfy the successive selection criteria, as well as the cumulative and relative fraction of events remaining after applying each criterion. The values in the last two lines of the "Relative" column are the fraction of events relative to the "m γ R γ R > 175 GeV" line. The values in the "Preselection" line include the offline Loose photon selection, E T > 25 GeV, |η| < 2.37, excluding the transition region between the barrel and endcap calorimeters, and the matching of the reconstructed photon to the photon trigger object applied to the two leading reconstructed photons. The label "Relative E T " denotes the requirements on E T /m γ R γ R for the reconstructed photons, described in Section 5.1.

N observed
Fraction of events Cumulative Relative All triggered events 6.4 × 10 9 --

Signal modeling
The DSCB function has been shown to be effective in modeling new-particle resonances expected to have a Gaussian core surrounded by asymmetric and non-Gaussian low-and high-mass tails, and is described in detail elsewhere [31]. In this analysis the DSCB is a function of the mass of two reconstructed photons (photon-jets in simulated signal samples), with parameters to account for the peak position and width of the Gaussian part, as well as for the upper and lower tails where the resonance shape meets the smoothly falling two-photon mass background. For the benchmark signal scenario investigated here, since the reconstructed photons are photon-jets (e.g., a → γγ and a → 3π 0 → 6γ), the reconstructed m γ R γ R corresponds to the mass of two a particles, i.e., the mass of the parent particle, X.
For the benchmark signal scenario, for very small values of m a /m X , the behavior of the DSCB as a function of the mass of two photon-jets is nearly identical to that of the BSM process X → 2γ. The position of the fitted peak of the DSCB is slightly lower than the mass input to the generator. With the NWA approach, the width of the Gaussian core σ CB is dominated by detector resolution, and it increases linearly with m X , from 2 GeV for m X = 200 GeV to 14 GeV for m X = 2 TeV. For larger values of m a /m X , the wider opening angle between the photons inside a photon-jet leads to a greater fraction of the energy of the shower leaking out of the window defined by the cells of the EM calorimeter to collect energy for the reconstruction of photons, leading to a further increase in the mass shift and width of the DSCB. For instance, for m X = 600 GeV and m a = 5 GeV, the width is σ CB = 8 GeV for the process X → aa → 4γ, and σ CB = 9 GeV for the process X → aa → 6π 0 . For a given m X and m a , the same signal mass shape modeling results are used for the analysis of the two orthogonal event categories (the low-∆E category and the high-∆E category), since only a small dependence of the signal mass distributions on ∆E is observed.
To validate the mass shape modeling results, injection tests are performed, where a fixed number of signal events are inserted into a pseudo-dataset reproducing a background-only m γ R γ R spectrum of one of the two event categories, and the number of events inserted is then compared with the number determined by fitting the DSCB. The pseudo-datasets are generated from background probability density functions (represented by Eq. (1), described in Section 6.2) with the parameters determined from a fit to the observed m γ R γ R spectra in collision data. For each simulated sample of the benchmark scenario, with different values of m a and m X , separate tests are performed for an increasing number of injected signal events. The average of the number of events determined from the fit to multiple pseudo-datasets and the number inserted should be identical in an ideal case, and the difference between these two numbers is taken as a systematic uncertainty in the signal mass shape modeling.
The fraction, f , of signal events that appear in the low-∆E category is parameterized as a function of the mass parameters m X and m a of the benchmark signal scenario, to have a continuous model for all the masses considered in the results. The values of f are taken from simulation and a third-order spline interpolation is performed as a function of m a /m X .
Similarly, the total signal selection efficiency, ε, is calculated from the individual signal mass points generated, and is parameterized as a function of m X and m a . This serves as an input to the calculation of the cross-section times branching ratios for the benchmark signal scenario.
The modeling of signal mass shape, f , and ε as functions of (m X , m a ) is performed separately for the two different final states of the benchmark signal scenario, X → aa → 4γ and X → aa → 6π 0 . In general, the results are similar for the two decay scenarios. The main distinction is in the different trend in f with respect to m X and m a , especially the threshold in m a /m X at which the values of f transition from f > 0.5 to f < 0.5. This threshold is found to be at m a /m X 0.0015 for X → aa → 4γ, and at m a /m X 0.0020 for X → aa → 6π 0 .

Background modeling
The backgrounds in this search mainly consist of the SM production of events containing either two prompt photons; one prompt photon and one hadronic jet; or two hadronic jets. Prompt photons are defined as photons not originating from hadron decays. Hadronic jets can be misreconstructed as a photon. The three background components are denoted γγ, γ j or jγ, and j j, respectively, with the first symbol indicating the one with a higher value of E T . The m γ R γ R distribution of the sum of these background components is described by an analytic function, separately for each of the two ∆E categories. The parameters of the two analytic functions are determined from fits to the m γ R γ R distributions in the analysis signal region of collision data from a lower edge of m γ R γ R = 175 GeV.
Based on simulated samples, the contribution from Drell-Yan processes, where the two isolated electrons are misreconstructed as photons, is expected to be at the sub-percent level in the analysis signal region. The shape of the m γ R γ R distribution of the Drell-Yan contribution in the mass range m γ R γ R > 175 GeV is expected to be similar to that of the γγ component, and it is therefore absorbed into the analytic function fit for the continuum background components.
The choice of the functional form describing the background distribution is based on studies of background templates. A variety of functional forms are considered for the background parameterization to achieve a good compromise between limiting the size of a potential bias toward the identification of a signal when none is present (the spurious signal) and retaining good statistical power. The size of the spurious signal for a given functional form is estimated by performing a maximum-likelihood fit to the background templates using the sum of signal and background parameterizations.
To determine the overall shape of the background mass spectra, background templates are determined using both the simulation and collision data, separately for each of the two ∆E categories. A simulated sample of prompt diphoton events is used to model the shape of the contribution from γγ events. Subsets of collision data that are similar but orthogonal to the signal region are used to determine the shapes of the γ j, jγ and j j components, where the subleading reconstructed photon, leading reconstructed photon, or both reconstructed photons, respectively, are required to fail the default isolation criterion but satisfy a looser one. This looser criterion is defined by loosening the requirement for the calorimeter isolation variable to E iso T < 7 GeV. The resulting samples of γγ, γ j, jγ, and j j are summed to derive the background templates, scaled with the background composition fractions determined from the matrix method described below.
The background composition of a given mass spectrum of two reconstructed photons is estimated using a matrix method [33], where events are categorized into four subsets by whether both, only the leading, only the subleading, or neither of the two leading reconstructed photons satisfy the calorimeter isolation requirement. The method relies on external estimates of the efficiency for prompt photons satisfying calorimeter isolation and the rate at which hadronic jets can mimic a photon satisfying calorimeter isolation (the "fake rate"). Photon isolation efficiency is estimated with simulated samples of prompt photons. The isolation variables of photons in simulated samples are adjusted by applying correction factors obtained from small differences observed between photon-enriched control samples of collision data and simulation. An uncertainty is assessed for the photon isolation efficiency by comparing the nominal efficiencies with those derived without applying the corrections to the isolation variable in simulated samples. Fake rates are determined using subsets of collision data with selection criteria imposed so that they are similar but orthogonal to the analysis signal region ("control regions"). These control regions are defined by requiring reconstructed photons to fail the baseline Loose photon selection but satisfy another, looser photon selection. This looser photon selection, with respect to the Loose selection, is defined by removing requirements on two additional shower shape variables that quantify the lateral shower development in the EM calorimeter second layer (w η2 and R φ , described in Table 1 of Ref. [30]). A difference of approximately 1 GeV is found between the isolation energy spectra in the signal and control regions. This is accounted for by shifting the threshold of the isolation selection criteria by ± 1 GeV, determining the resulting change in the calculated fake rates, and assigning the difference as a systematic uncertainty in these values.
An additional uncertainty is assessed by altering the definition of the control regions. To accomplish this, a looser photon selection, with respect to the Loose selection, is defined by removing the requirement on one shower shape variable (w η2 ) instead of two and comparing the difference between the resulting fake rates.
The resulting background compositions are shown in Table 3. Good agreement is seen between the observed isolation spectrum and the expected spectrum based on the matrix method results, within uncertainties, as shown in Figure 4.
The background templates are derived with the summation of the γγ, γ j, jγ and j j components scaled by the background composition fractions, separately for each of the two ∆E categories, as described above. The resulting background templates are presented in Figure 5.
To evaluate the size of the spurious signal, a test is performed using these background templates and the signal modeling described in Section 6.1. The background templates are normalized to the integrated  , is chosen to describe the shape of the m γ R γ R distribution:     The variable x is defined as x = m γ R γ R / √ s. The parameters a and b j are free parameters and N is the normalization factor. The spurious signal tests are then performed using a maximum-likelihood fit of the sum of the signal and background parameterizations to each of the two background templates. The spurious signal is allowed to be negative as well as positive. The final functional form used to model the background when performing the search for resonances is one where the estimated spurious signal is required to be smaller than 30% of the statistical uncertainty in the fitted signal yield across the full mass spectrum. The cutoff of 30% is chosen to ensure that the contribution of this systematic uncertainty to the total uncertainty, including all statistical and systematic uncertainties, is subdominant and smaller than 5%.
The method is validated by checking that similar results are obtained when the test is performed using variations of the background templates, for which the background compositions are shifted within the uncertainties presented in Table 3. When the fraction of the γγ component is shifted up and those for γ j, jγ, j j are shifted down, or vice versa, the size of the resulting spurious signals are consistent within the statistical uncertainty of the background templates.
The resulting functional form used for the background mass spectrum evaluation of the two categories is shown in Eq. (1) with k = 1. Figure 5 shows the level of agreement between this functional form and the background templates. The resulting background model and its associated systematic uncertainties are used when searching for resonances in the mass spectra of the signal region.

Systematic uncertainties
Several sources of systematic uncertainties that affect the determination of the signal yield are taken into account. In most cases, systematic uncertainties are smaller than statistical errors.
The uncertainty in the combined 2015+2016 integrated luminosity is 2.1%. It is derived, following a methodology similar to that detailed in Ref. [35], from a calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016.
The impact of the photon energy resolution on signal modeling is evaluated. It mainly affects the mass shape width, σ CB , of the Crystal Ball function used to model the signal mass shape. The photon energy resolution is adjusted by one standard deviation from the nominal value in both positive and negative directions, and the resulting change in the fitted signal width is determined. The relative difference in the fitted value of σ CB ranges from as small as a few percent to as large as 37%, increasing with larger m X and dependent slightly on m a , and is taken as a systematic uncertainty.
Systematic uncertainties in the extracted signal yield due to signal mass shape modeling are evaluated via injection tests, described in Section 6.1. The final fitted values of the number of signal events deviate from the injected values by less than 1% almost everywhere, rising to a maximum of 5% for some signal mass values at the edge of the analysis search region of m a = 0.1 GeV for the high-∆E category and m a = 0.01 × m X for the low-∆E category. This is taken as the estimate of the systematic uncertainty in the signal yield.
Uncertainties in the modeling of the category fraction, f , are evaluated by an envelope to cover the deviations of the values of f from simulation and the parameterization. The absolute value of the change in f varies as a function of m a /m X , from 3% at m a /m X = 0, increasing to 12%-14% at around m a /m X = 0.002, and decreasing to 6%-10% at 0.002 < m a /m X < 0.01. This is taken as the estimate of the systematic uncertainty in f .
Other systematic uncertainties in the extracted signal yield and the migration of signal events between the two orthogonal ∆E categories are evaluated by comparisons between nominal and systematically varied versions of various experimental uncertainty sources, such as the photon energy scale and resolution, isolation selection efficiency, shower shape modeling, and pileup. The systematic uncertainties due to the photon energy scale and resolution are adapted from results determined during LHC Run 1 [32], with minor updates derived from data-driven corrections determined using Run 2 data. Uncertainties related to the Loose photon identification scheme are evaluated with the systematic variations for the shower shape modeling, without the correction factors applied to simulation derived from small differences observed between photon-enriched control samples of collision data and simulation [36]. The uncertainty in the photon calorimeter isolation efficiency is calculated from changes due to applying and not applying corrections derived from small differences observed between photon-enriched control samples of collision data and simulation. The uncertainties of the efficiency correction factors using photon-enriched control samples of collision data are used to derive the uncertainty in the photon track isolation efficiency. The pileup uncertainty is taken into account by propagating it through the event selection. The uncertainties in ε and f due to these sources for the mass regions considered for the benchmark signal scenario are calculated. The uncertainties are less than 1% in almost all cases, rising to ∼4% for some isolation and shower shape uncertainties for larger values of m a /m X at the edge of the analysis sensitivity.
Additional systematic uncertainties in the loose diphoton trigger efficiency are not assessed. The E T requirements for reconstructed photons are much larger than the value at which the diphoton trigger utilized becomes nearly 100% efficient, and any additional uncertainties in signal efficiency due to mismodeling of the trigger-level shower shape variables are accounted for when calculating uncertainties in offline Loose identification, because the loose photon identification definitions at the trigger and offline levels are strongly correlated.
The uncertainty in the signal kinematic acceptance, which is included in the definition of the total signal selection efficiency, is evaluated for the choice of PDF set used for the simulation of the signal samples. It is less than 1% in most cases, rising to ∼ 4% for large m X around m X ∼ 2 TeV.
The systematic uncertainties related to the evaluation of the background mass spectrum are determined from the spurious signal method, described in Section 6.2. The spurious signal as a function of m X and m a is parameterized so that the modeling between mass points is continuous. This parameterization is performed in such a way that it can slightly overestimate the size of spurious signals, especially at the lower end of the m X range, m X = 200 GeV. The size of the parameterized spurious signal decreases for larger m X and depends slightly on the m a value, ranging from 85 to 6 × 10 −3 events for the low-∆E category, and from 32 to 1 × 10 −2 events for the high-∆E category.
The systematic uncertainties are generally smaller than the statistical errors, with the systematic uncertainty in the background evaluated from the spurious signal being the largest contribution. This is because the parameterization of the size of spurious signals slightly overestimates the values at the lower end of the m X range, as described above. The impact of the systematic uncertainties on the expected limit decreases with the resonance mass m X from 51% at most for m X = 200 GeV to 5% at most for m X > 800 GeV. The impact of the systematic uncertainties on the signal yield obtained from the fit is summarized in Table 4. Table 4: Breakdown of the relative contributions to the total uncertainty in the signal yield obtained from the fit. For each source of uncertainty σ source , the fraction σ source /σ total is presented, where σ total is the total uncertainty that includes the statistical uncertainty. The sum in quadrature of the individual components differs from 100% due to small correlations between the components. The values here are for the signal process X → aa → 4γ. The mass points (m X , m a ) = (200 GeV, 0.3 GeV), (600 GeV, 0.9 GeV) correspond to those values for which the systematic uncertainty of the category fraction f is the highest. Similar results are found for the decay X → aa → 6π 0 .

Statistical procedure
For a given fixed signal mass hypothesis, a mass spectrum fit including both the background and signal components is performed to the full mass spectrum of m γ R γ R > 175 GeV, using an unbinned maximumlikelihood approach, simultaneously for the two categories (low-∆E and high-∆E categories). A constraint is placed on the ratio of the two separate normalization factors of the signal component for the two categories, evaluated from the category fraction f , which depends on the signal masses m X and m a . Deviations from the background-only hypothesis are searched for starting from m X = 200 GeV, and the entire m γ R γ R range is used for the background component for each hypothesis test. The p-values are calculated with the profile likelihood ratio as the basis for the test statistic and utilizing an asymptotic approximation [37].
Systematic uncertainties (described in Section 7) are treated as nuisance parameters in the likelihood function, where each is a floating parameter constrained by either a Gaussian function (for spurious signal and uncertainties related to the migration of events between the ∆E categories) or a log-normal function (for all other uncertainties). Two nuisance parameters are introduced for the extracted signal yield due to signal mass shape modeling uncertainties, one for each ∆E category, and they are multiplied by the signal normalizations of each category. One nuisance parameter is introduced for the impact of the photon energy resolution on the mass shape width, and it is multiplied by the signal mass shape width σ CB . One nuisance parameter is introduced for the modeling of the category fraction, f , which is added to f to shift its value. Several nuisance parameters are introduced for experimental uncertainty sources and PDF uncertainty that affect the extracted signal yield, total signal selection efficiency, ε, and category fraction, f . Two nuisance parameters are introduced for the spurious signals, one for each ∆E category. For a given signal mass (m X , m a ) hypothesis, the spurious signals are given the same m γ R γ R shape as the signal component, and normalized by the size of the spurious signals.
The calculation of p-values for the background-only hypothesis (p 0 ) is performed for a narrow resonance from m X = 200 GeV to m X = 2.7 TeV, with a scan step of 1 GeV. Since the samples for the benchmark signal scenario were simulated for the m X values in the range 200 GeV < m X < 2 TeV, the results of signal mass shape modeling, modeling of category fraction f , and systematic uncertainties are extrapolated for the m X values in the range 2 TeV < m X < 2.7 TeV.
Expected and observed upper limits, at the 95% confidence level (CL), on the production cross-section times the product of branching ratios are calculated as a function of the mass parameters of the benchmark signal scenario, m X and m a , following the CL s modified frequentist prescription [38]. Upper limits are determined separately for the two final states of the benchmark signal scenario where the a particle decays into either a pair of photons or three neutral pions.
The assumptions inherent in the use of the asymptotic approximation are validated by sampling distributions of the test-statistic using pseudo-experiments, for a few signal mass points. The asymptotic approximation yields median values of the expected upper limits within 5% of those calculated with a large number of pseudo-experiments for most of the values of m X and m a tested. Due to the small number of events in data in the region m γ R γ R > 1 TeV in the high-∆E category, larger deviations are observed for m X > 1 TeV and large m a . The deviation is smaller than 5% at (m X , m a ) = (1 TeV, 10 GeV), but the expected upper limits obtained from the asymptotic approximation are smaller than those from pseudoexperiments by 20% for (m X , m a ) = (1.5 TeV, 10 GeV), and 30% for (m X , m a ) = (2 TeV, 10 GeV).

Results
The observed m γ R γ R spectra in the signal region are shown in Figure 6. The results of the two-dimensional scan of p 0 , equivalently expressed in terms of the local significance-the number of standard deviations away from the mean of a normal distribution-are shown in Figure 7. Two different regimes can be seen in this plot, above and below the threshold at m a ∼ 0.0015 × m X . These are a result of the categorization of events based on the ∆E variable. For m a 0.0015 × m X , a larger fraction of signal events is expected in the low-∆E category, and for m a 0.0015 × m X , a larger fraction of signal events is expected in the high-∆E category. The largest local deviation from the background-only hypothesis is found to be 2.7σ, corresponding to m X = 729 GeV and m a = 0.1 GeV for the decay X → aa → 4γ. The width of the signal mass shape for m X = 729 GeV and m a = 0.1 GeV is 6 GeV, and thus this deviation appears as a small area in Figure 7. A small excess of events is also observed centered around m X = 1.1 TeV and m a = 7 GeV, which corresponds to a local deviation of 2.2σ. The observed maximum local deviation is less significant than the median of the largest deviation obtained in background-only pseudo-experiments, calculated in the search region defined by m X values from 200 GeV to 2.7 TeV and m a values from 0.1 GeV to 0.01×m X . The m γ R γ R mass distribution is found to be consistent with the background-only hypothesis.
The 95% CL observed and expected upper limits on the cross-section for the production via gluon-gluon fusion of a high-mass scalar particle, X, with narrow width times the branching ratios into a pair of a particles and the subsequent decay of each a into a pair of photons, σ X × B(X → aa) × B(a → γγ) 2 , are shown in Figure 8, separately for different values of m a . The same result is presented in Figure 9, with the ratio m a /m X shown on the horizontal axis. This plot illustrates the two features of this search. First, when the ratio m a /m X is larger than a threshold of roughly 0.0015, more signal events are expected in the high-∆E category, which has a significantly better signal-to-background ratio compared with the low-∆E category, thus leading to stronger upper limits. Second, for larger values of m a /m X , the decrease in the signal selection efficiency leads to weaker upper limits.
The 95% CL observed and expected upper limits on the cross-section times product of branching ratios for the decay of the a into three neutral pions, σ X × B(X → aa) × B(a → 3π 0 ) 2 , are shown in Figure 10, separately for different values of m a . This result shows features similar to that shown in Figure 8, with slight differences arising mainly from the different trend of the category fraction, f , with respect to the mass values m X and m a .  Figure 6: Observed distributions of the mass of two reconstructed photons passing all analysis selections, m γ R γ R , for the two signal region categories. The background-only fit result is superimposed. The ±1σ uncertainty originating from the uncertainties in the fit function parameter values is shown as a shaded band around the fit. The lower panel of each plot displays the significance associated with the observed event yield in each bin, calculated before considering systematic uncertainties. The calculation assumes that the event yield in each bin is Poisson-distributed with a mean given by the background-only fit. The computation is performed with a one-sided test based on the positive or negative tail of the Poisson distribution, depending on the sign of the difference between the event yield and the fit estimate, with negative significance values quoted for negative differences [39].  Figure 8: The observed and expected upper limits on the production cross-section times the product of branching ratios for the benchmark signal scenario involving a scalar particle X with narrow width decaying via X → aa → 4γ, σ X × B(X → aa) × B(a → γγ) 2 . The limits are calculated using the asymptotic approximation. This leads to an underestimate of the limits, especially for m X > 1 TeV and large m a . The limits for m a = 5 GeV and 10 GeV do not cover as large a range as the other mass points, since the region of interest is limited to m a < 0.01 × m X .   Figure 9: The observed and expected upper limits on the production cross-section times the product of branching ratios for the benchmark signal scenario involving a scalar particle X with narrow width decaying via X → aa → 4γ, σ X × B(X → aa) × B(a → γγ) 2 . They are evaluated as a function of m a /m X for fixed values of m X . The limits are calculated using the asymptotic approximation. This leads to an underestimate of the limits, especially for m X > 1 TeV and large m a . The results for the X → aa → 6π 0 case are qualitatively similar.  Figure 10: The observed and expected upper limits on the production cross-section times the product of branching ratios for the benchmark signal scenario involving a scalar particle X with narrow width decaying via X → aa → 6π 0 , σ X × B(X → aa) × B(a → 3π 0 ) 2 . The limits are calculated using the asymptotic approximation. This leads to an underestimate of the limits, especially for m X > 1 TeV and large m a . The limits for m a = 5 GeV and 10 GeV do not cover as large a range as the other mass points, since the region of interest is limited to m a < 0.01 × m X .

Conclusion
A search for pairs of highly collimated groupings of photons-photon-jets-that are identified as single, photon-like energy clusters in the EM calorimeter of the ATLAS detector at the LHC is presented. Data from proton-proton collisions at a center-of-mass energy of 13 TeV collected in 2015 and 2016, corresponding to an integrated luminosity of 36.7 fb −1 , are used. Pairs of photon-jets can arise, for example, as the final-state decay products of a new high-mass resonance decaying via new light resonances into highly collimated groupings of photons. Candidate photon-jet events are initially selected with a loose diphoton trigger and then potential photon-jets are selected using a combination of variables that model EM shower development. Sensitivity to photon-jets is then increased by categorizing reconstructed photons by one of those shower shape variables and narrow resonances are searched for in the resulting mass distributions of two reconstructed photons. The observed mass spectra are consistent with the SM background expectation.
The results are interpreted in the context of a BSM scenario containing a high-mass scalar particle with narrow width, X, that decays into photon-jets via low-mass intermediate particles with spin 0, a. For the range of m X investigated, from 200 GeV to 2 TeV, upper limits on σ × B(X → aa) × B(a → γγ) 2 are found to range from 0.2 fb to 1 fb over most of the range of m X , for 100 MeV < m a < 2 GeV, rising to 10-100 fb for values of m X at the low end of the range, depending upon m a . Similarly, upper limits on σ × B(X → aa) × B(a → 3π 0 ) 2 are found to range from 0.2 fb to 1 fb over most of the range of m X , for 500 MeV < m a < 2 GeV, rising to 10-100 fb for values of m X at the low end of the range. These limits are calculated using an asymptotic approximation. In addition to the calculated upper limits for this benchmark signal scenario, the results, including the evaluation of the observed upper limits, are provided in HepData [40] in a largely model-independent way, to enable reinterpretation in the context of other signal models containing highly collimated photon-jets of low or high photon multiplicity.