Search for new phenomena in events with at least three photons collected in pp collisions at $\sqrt{s}$ = 8 TeV with the ATLAS detector

Abstract

Results of a search for new phenomena in events with at least three photons are reported. Data from proton–proton collisions at a centre-of-mass energy of 8 TeV, corresponding to an integrated luminosity of 20.3 fb$^{-1}$, were collected with the ATLAS detector at the LHC. The observed data are well described by the Standard Model. Limits at the 95 % confidence level on new phenomena are presented based on the rate of events in an inclusive signal region and a restricted signal region targeting the rare decay $Z\rightarrow 3\gamma$, as well as di-photon and tri-photon resonance searches. For a Standard Model Higgs boson decaying to four photons via a pair of intermediate pseudoscalar particles (a), limits are found to be $\sigma \times {\text{ BR }}(h \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2} < 10^{-3} \sigma _{\text {SM}}$ for 10 GeV $< m_{a}<$ 62 GeV. Limits are also presented for Higgs boson-like scalars (H) for $m_{H} >$ 125 GeV, and for a $Z'$ decaying to three photons via $Z' \rightarrow a+\gamma \rightarrow 3\gamma$. Additionally, the observed limit on the branching ratio of the Z boson decay to three photons is found to be BR$(Z \rightarrow 3\gamma ) < 2.2 \times 10^{-6}$, a result five times stronger than the previous result from LEP.

1 Introduction

Many extensions of the Standard Model (SM) include phenomena that can result in final states consisting of three or more photons. Extensions of the SM scalar sector [1–5], for example, often include pseudoscalar particles (a) with couplings to the Higgs boson [6, 7] (h) and branching ratios into photons that would be visible at the LHC, in addition to scalars (H) with masses different from the SM-like Higgs boson of $m_{h} =$ 125 GeV that can also decay via $H \rightarrow aa \rightarrow 4\gamma$. Other models feature additional vector gauge bosons that can decay to a photon and a new pseudoscalar boson, a, with the subsequent decay of the a into a pair of photons, resulting in a three-photon final state [8]. Moreover, in the SM, the Z boson can decay to three photons via a loop of $W^{\pm }$ bosons or fermions. The decay is heavily suppressed and the branching ratio is predicted to be $\sim$5$\times$10$^{-10}$ [9]. The current most stringent bound on this process comes from the L3 Collaboration, which placed a limit of BR$(Z \rightarrow 3\gamma ) < 10^{-5}$ [10]. The ATLAS detector has collected $\sim$10$^{9}$ Z boson events, and thus an observation of this decay would indicate an enhancement of this decay rate and could be evidence of phenomena not predicted by the SM. Feynman diagrams for some of these beyond-the-Standard Model (BSM) and rare SM scenarios are shown in Fig. 1.

Fig. 1

Feynman diagrams for possible beyond-the-Standard Model (top) and rare Standard Model (bottom) scenarios that result in final states with at least three photons

To ensure sensitivity to these and other possible rare SM and BSM scenarios, an inclusive three-photon search is performed using 20.3 fb$^{-1}$ of LHC proton-proton collisions collected by the ATLAS detector in 2012 at a centre-of-mass energy of 8 TeV. Such a model-independent search is the first of its kind, as are the interpretations for a Higgs boson decaying to four photons via two intermediate pseudoscalar a particles (for a Higgs boson of $m_{h} =$ 125 GeV and for Higgs-like scalars of higher masses) and for three-photon resonances corresponding to a new vector gauge boson.

The dominant backgrounds include the irreducible component with three or more prompt photons, as well as the reducible components consisting of combinations of photons and electrons or hadronic jets misidentified as photons. The contributions from events with jets which are misidentified as photons are calculated from data-driven methods, while simulation is used to estimate the contributions from the irreducible background and the reducible background originating from electroweak processes that lead to electrons which are misidentified as photons in the detector. Collision data is used to derive corrections to the probability obtained from simulation that electrons are misidentified as photons.

2 The ATLAS detector

The ATLAS experiment [11] at the LHC is a multi-purpose particle detector with a forward-backward symmetric cylindrical geometry and a near $4\pi$ coverage in solid angle.¹ It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range $|\eta | < 2.5$. It consists of silicon pixel, silicon micro-strip, and transition radiation tracking detectors. Lead/liquid-argon (LAr) sampling calorimeters provide EM energy measurements with high granularity. A hadronic (iron/scintillator-tile) calorimeter covers the central pseudorapidity range ($|\eta | < 1.7$). The end-cap and forward regions are instrumented with LAr calorimeters for both EM and hadronic energy measurements up to $|\eta | = 4.9$. The muon spectrometer surrounds the calorimeters and is based on three large air-core toroid superconducting magnets with eight coils each. Its bending power ranges from 2.0 to 7.5 T m. It includes a system of precision tracking chambers and fast detectors for triggering. A three-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 reduce the accepted rate to at most 75 kHz. This is followed by two software-based trigger levels that together reduce the accepted event rate to 400 Hz on average depending on the data-taking conditions during 2012.

3 Event and object selection

This search utilises a three-photon trigger that places a minimum requirement on the photon momentum in the plane transverse to the beam axis (transverse momentum, or $p_{\text {T}}$) of 15 GeV, applied on three photon candidates in the EM calorimeter. Each candidate is additionally required to satisfy a set of loose photon identification criteria [12]. Stringent detector and data quality criteria are applied offline. Events are required to contain at least one interaction vertex, with no additional vertex requirements.

Photon candidates must satisfy a pseudorapidity requirement of $|\eta |$ < 2.37, excluding the transition region between the barrel and end-cap of 1.37 $< |\eta |<$ 1.52, and must satisfy requirements on the shape of the energy deposit in the calorimeter. A photon candidate is rejected if the barycentre of its energy deposit is within a cone of $\Delta R \equiv \sqrt{(\Delta \eta )^{2} + (\Delta \phi )^{2}}<$ 0.15 around the barycentre of the energy deposit of a higher $p_{\text {T}}$ photon candidate. Finally, selected photon candidates are required to satisfy a more stringent set of identification criteria, known as tight [12]. Photon isolation is defined by the amount of transverse energy, $E_{\text {T}}^{\text {iso}}$, deposited in the EM calorimeter within a cone of size $\Delta R$ around the photon candidate, excluding the energy of the photon candidate itself. It is a powerful means of distinguishing between photons and hadronic jets misidentified as photons, since the energy clusters deposited by photons in the EM calorimeter tend to be narrower in the transverse direction than those deposited by jets. Because minimum-bias proton–proton interactions in the same or nearby bunch crossings (pileup) can affect the calculated photon isolation energy, a correction is applied based on an event-by-event energy density pileup estimation. This search uses an isolation cone of size $\Delta R<$ 0.4, and a correction to the $E_{\text {T}}^{\text {iso}}$ value of a photon candidate is made when another photon candidate passing the tight identification criteria is found within an annulus of 0.15 $< \Delta R<$ 0.4 around the photon candidate. The correction consists of subtracting the $p_{\text {T}}$ value of the other photon candidate found within the annulus from the $E_{\text {T}}^{\text {iso}}$ value of the photon candidate under consideration. The final isolation criterion is $E_{\text {T}}^{\text {iso}} <$ 4 GeV.

Events in the inclusive signal region are required to have at least three tightly identified and isolated photon candidates, where the two photon candidates with the highest transverse momentum must have $p_{\text {T}}>$ 22 GeV while the third highest must have $p_{\text {T}}>$ 17 GeV. The restricted signal region, targeting the rare decay $Z \rightarrow 3\gamma$, is a subset of these events where an additional criterion of 80 GeV $< m_{3\gamma }<$ 100 GeV is placed on the invariant mass of the three-photon system. The signal regions are supplemented by several control regions, where at least one of the photon candidates fails the isolation requirement.

4 Simulated event samples

Simulated event samples are used to estimate several SM background processes in the search for excesses in the inclusive signal regions, as well as to model signal predictions for both the inclusive searches and the resonance searches for the specific BSM scenarios considered here. The simulated BSM signal samples are also used to define a fiducial region for which the search criteria are largely model-independent.

4.1 Simulated backgrounds

The SM two-photon process is an irreducible background to a three-photon search because a third photon may arise from minimum-bias proton–proton interactions in the same bunch crossing. The SM two-photon background is simulated with Pythia 8 [13], and the three- and four-photon backgrounds are simulated with MadGraph 5 [14], with Pythia 8 used for fragmentation and hadronisation. The production of two, three and four photons in the SM contains large contributions from higher-order Feynman diagrams. Thus, the two-, three- and four-photon simulated event samples calculated at leading order (LO) are multiplied by factors (“K-factors”) determined from studies with generators that include next-to-leading order (NLO) contributions, namely MCFM 6.8 [15] and VBFNLO 2.7.0 [16–18], using the parton distribution function (PDF) sets CTEQ6L1 [19] for the LO cross sections and CT10 [20] and MSTW8NL [21] for the NLO cross sections. These K-factors are 1.9 ± 0.2 for the two-photon process and 3.3 ± 0.5 for the three-photon process. The four-photon process is not included in NLO generators, and since the four-photon background is $\sim$10$^{-3}$ of the total background expectation in the inclusive signal region, the three-photon K-factor is applied to the four-photon background sample. The uncertainties for these K-factors are determined by multiplying the renormalisation and factorisation scales independently by 2 and 0.5 and taking the largest deviations from the nominal value of the K-factor.

The reducible backgrounds where electrons are misidentified as photons originate from multiple sources. Processes where a Z boson decays to an $e^{+}e^{-}$ pair, accompanied by a photon not from the matrix element, are modelled with Powheg-Box 1.0 [22], using Pythia 8 for fragmentation and hadronisation, and Z+$\gamma$ production is modelled with Sherpa 1.4.1 [23]. Backgrounds from processes involving the leptonic decay of the W boson in association with photons and/or hadronic jets are simulated with Alpgen [24] and Herwig [25, 26]. Possible mis-measurement of the rate of electrons misidentified as photons in simulation is addressed by comparing to electrons misidentified as photons from $Z \rightarrow e^{+}e^{-}$ events in data.

To obtain estimates of the rates at which true photons populate the regions of kinematic phase space assumed to be dominated by jets (used in the calculation of the systematic uncertainty for the data-driven estimate of jet backgrounds), a sample containing events with one hard-process quark or gluon and one prompt photon is simulated using Pythia and the CTEQ6L1 PDF set.

The PDF sets for the simulated event samples of background processes used for the final background estimate in the inclusive search, for the MadGraph, Pythia and Alpgen + Herwig samples, are taken from CTEQ6L1, while for the Powheg-Box and Sherpa samples the PDF sets are taken from CT10.

4.2 Simulated signal processes

The $Z\gamma \gamma \gamma$ effective vertex has been implemented with FeynRules [27, 28] and then used in a customised MadGraph 5 model which is employed to simulate events, using the CTEQ6L1 PDF set and Pythia 8 for fragmentation and hadronisation. Each of the two non-trivial, independent, lowest-order effective Lagrangians for this process [29] contains a dimensionful coupling constant, and the values of these constants have been calculated [30] using the SM expected $Z \rightarrow 3\gamma$ branching ratio of 5.41$\times 10^{-10}$. These SM values are used in the simulation. The BSM process of a Higgs boson produced via gluon fusion and decaying to four photons via a pair of intermediate a particles is simulated with Powheg-Box and Pythia 8 (using the CT10 NLO PDF set). The BSM process of a new vector gauge boson decaying to three photons via $Z' \rightarrow a+\gamma \rightarrow 3\gamma$ is simulated with Pythia 8 (using the MSTW2008LO [21] leading-order PDF set).

4.3 Minimum-bias interactions and the ATLAS detector simulation

Minimum-bias proton–proton interactions in the same or nearby bunch crossings (pileup) are modelled with Pythia 8, using the MSTW2008LO PDF set. These pileup events are overlaid onto the hard-scattering process for all simulated signal and background samples to reproduce the distribution of the average number of interactions per bunch crossing observed over the course of data-taking in 2012.

All signal and background samples are processed with the full ATLAS detector simulation [31] based on Geant 4 [32] and reconstructed using the same software as that used for collision data.

5 Background composition estimate

The backgrounds in the search for excesses in the inclusive signal regions are estimated from a combination of simulated samples (detailed in the previous section) and methods employing collision data. The dominant backgrounds in the inclusive signal region are the irreducible SM two-, three- and four-photon processes, while for the $Z \rightarrow 3\gamma$ search channel, backgrounds involving electrons misidentified as photons are dominant.

5.1 Backgrounds estimated from simulation

The irreducible SM two-, three- and four-photon backgrounds, as well as backgrounds from processes involving electrons in the final state originating from Z decays and those involving the leptonic decay of the W boson in association with photons and/or hadronic jets, are estimated via simulation. The third photon for the SM two-photon background process typically arises from pileup interactions, but can occasionally be a quark- or gluon-initiated jet radiated from the incoming partons which is misidentified as a photon. Possible double-counting with the $2\gamma$ + 1-jet final state (estimated via a data-driven method described in the following section) is avoided by omitting from consideration events in the SM two-photon simulated sample where one of the three photon candidates is a jet, using generator-level information. Possible mis-measurement of the rate of electrons misidentified as photons in simulation is addressed by comparing $Z \rightarrow e^{+}e^{-}$ processes in simulation and in data. The per-electron scale factor is the ratio of the misidentification rate determined in data to that determined in simulated samples. This scale factor is independent of $p_{\text {T}}$ and $\eta$ for the ranges considered here, and is found to be 1.03 ± 0.04.

5.2 Data-driven estimates of $2\gamma$ + 1-jet, $1\gamma$ + 2-jet, and 3-jet backgrounds

A crucial aspect of the analysis is the data-driven estimate of the backgrounds where hadronic jets are misidentified as photons (hereafter called “jet fakes”), i.e., SM processes that can produce $2\gamma$ + 1-jet, $1\gamma$ + 2-jet, and 3-jet events. Collision data are used to derive efficiencies for photons passing the isolation criterion ($\epsilon _{\gamma }$) and rates at which jets are misidentified as isolated photons ($f_{\text {jet}}$). These values of $\epsilon _{\gamma }$ and $f_{\text {jet}}$ are then used in a likelihood matrix method (described below) to estimate the jet backgrounds.

A sample of photon candidates consisting mainly of jet fakes is defined in the following way. The standard tight and loose photon identification categories are augmented with a medium definition [33], intermediate between tight and loose. The medium photon is defined by relaxing some EM shower shape requirements that provide high levels of rejection of jet fakes. When the medium definition is combined with a further requirement that the photon candidates fail tight (the combination hereafter called non-tight), the result is a sample of photon candidates that is primarily composed of jet fakes. This method presupposes that the $E_{\text {T}}^{\text {iso}}$ distribution of the non-tight sample is composed primarily of jet fakes, and that the subset of tight photons with higher values of $E_{\text {T}}^{\text {iso}}$ (the “tail”, here for $E_{\text {T}}^{\text {iso}} >$ 7 GeV) is dominated by jet fakes. Under these assumptions, the tail of the non-tight distribution is scaled to match that of the tight distribution, thus providing a determination of the contribution of jet fakes to the signal region, i.e., the collection of photons that pass both the tight and the isolation criteria. The scaled non-tight distribution is then subtracted from the tight distribution. Photon isolation efficiency, $\epsilon _{\gamma }$, is then calculated as the ratio of the number of isolated photons (those that satisfy $E_{\text {T}}^{\text {iso}}<$ 4 GeV) to the total number of photons in the tight distribution after this subtraction has been performed. The rate at which jet fakes are identified as photons, $f_{\text {jet}}$, is the ratio of the number of isolated photons to the total number of photons in the non-tight distribution.

The assumptions described above are validated using simulated samples of events containing photons and jets, described in Sect. 4. Any collection of photon candidates consists of some combination of actual photons, which can be defined as “true”, and other objects that are misidentified as photons. The non-zero true photon contamination in the set of non-tight photons, and the set of tight photons that fail the isolation criterion, is taken from the simulated samples and is used to derive a systematic uncertainty (described in Sect. 6) on the jet background estimate procedure.

The procedure is performed separately for three kinematic regions as follows. Photons are ordered by $p_{\text {T}}$, highest to lowest. Three regions in the $p_{\text {T}}$–$\eta$ plane are defined as (1) 15 GeV $< p_{\text {T}}<$ 40 GeV and $|\eta | < 1.37$, (2) $p_{\text {T}} >$ 40 GeV and $|\eta | < 1.37$, and (3) $1.52< |\eta | < 2.37$. The separation into lower and higher $p_{\text {T}}$ bins around 40 GeV is chosen because this is the value at which $\epsilon _{\gamma }$ and $f_{\text {jet}}$ are changing rapidly, and the three regions were chosen to maintain a large number of events in each bin. The values of $\epsilon _{\gamma }$ and $f_{\text {jet}}$ are then calculated for each of the three regions, and the results are shown in Table 1.

Table 1 Photon isolation efficiencies ($\epsilon _{\gamma }$) and rates of jets misidentified as photons ($f_{\text {jet}}$) from collision data for the three kinematic regions, used for the jet background estimate. The isolation criterion is $E_{\text {T}}^{\text {iso}}<$ 4 GeV. The three regions were chosen to maintain a large number of events in each bin. The first uncertainty is statistical while the second is systematic

The data-derived $\epsilon _{\gamma }$ and $f_{\text {jet}}$ values are applied to events with three photon candidates to estimate the SM $2\gamma$ + 1-jet, $1\gamma$ + 2-jet, and 3-jet backgrounds. This is done using a likelihood-based version of a standard matrix method (here called the “likelihood matrix method”). In standard matrix methods [33], a matrix of efficiencies relates an observed event that falls into a particular event category (based on some discriminating variable or variables) to the true, unknown final states to which the event has the possibility of corresponding, and the matrix is inverted to determine probabilities that a given observed event corresponds to one of these true final states. When summed over a large number of events, these per-event estimators average to the overall estimate of the number of events in each true final state.

In the likelihood matrix method, by contrast (and with respect to the present three-photon search), the expected yield for each three-object final state consisting of jets plus photons or all jets is the result of fitting a likelihood function to data. For the event sample where all three photons, ordered from highest to lowest $p_{\text {T}}$, have satisfied the tight requirements, events are placed into 160 orthogonal categories designated by six criteria. These are defined first by the three regions in the $p_{\text {T}}$–$\eta$ plane to which each photon candidate belongs. These are the same three regions that are used to categorise photons and to calculate $\epsilon _{\gamma }$ and $f_{\text {jet}}$, described and labeled previously as regions 1–3. The remaining three criteria by which each event is categorized are three boolean variables, one for each photon candidate, indicating whether it passed or failed the isolation criterion. Since each of the three photon candidates either passes (P) or fails (F) isolation, there are $2^{3} = 8$ possible isolation combinations for three photons: PPP, PPF, PFP, FPP, PFF, FPF, FFP, and FFF. The three photons are ordered by $p_{\text {T}}$, from highest to lowest, and, since one of the three kinematic regions defined above depends only upon $\eta$, there are twenty possible $p_{\text {T}}$–$\eta$ bin combinations: 333, 332, 323, 322, 331, 313, 311, 321, 222, 223, 232, 233, 221, 211, 231, 213, 111, 113, 131, 133. This results in $8\times 20 = 160$ categories, denoted PPP_333 for those events where the three photon candidates all passed isolation and had $p_{\text {T}}$–$\eta$ values placing them in the “3” kinematic region, PPF_321 for those events where the leading and subleading photons passed isolation and the sub-subleading photon failed isolation, and the $p_{\text {T}}$–$\eta$ value combinations placed them successively in the “3”, “2”, and “1” regions, etc.

Each of the 160 categories corresponds to a Poisson function where the observed number of events is the number of events seen in data for that category and the expected number of events is a sum of terms corresponding to each of the possible true (unknown) final states consisting of photons and jets or only jets for a particular $p_{\text {T}}$–$\eta$ combination. Each term in a given sum is multiplied by the appropriate values of $\epsilon _{\gamma }$ and $f_{\text {jet}}$. A likelihood is then constructed consisting of a product of the 160 Poisson functions. The expectations for each true final state are the maximum likelihood estimators that result from fitting this likelihood function to the data. That is, the true unknown expectations are allowed to float in the fit and are constrained to be positive and, hence, physical. The estimated number of events of a given final state in a particular signal or control region – defined by whether the photons passed or failed isolation – is determined by summing the resulting expectations from the fit times the appropriate $\epsilon _{\gamma }$ and $f_{\text {jet}}$ values.

6 Systematic uncertainties

6.1 Data-driven uncertainties

For the data-driven jet background estimate, systematic uncertainties arise in the calculation of the rate of photons passing the isolation criterion, $\epsilon _{\gamma }$, and the rate of jets misidentified as isolated photons (“jet fakes”), $f_{\text {jet}}$. This calculation relies upon the assumptions that both the tail ($E_{\text {T}}^{\text {iso}} >$ 7 GeV) of the $E_{\text {T}}^{\text {iso}}$ distribution of tight photons and the entirety of the $E_{\text {T}}^{\text {iso}}$ distribution of non-tight photons are primarily composed of jet fakes. Tests on simulations of photons and jets indicate that the true photon contamination in these jet-dominated regions is between 5 and 15 %, depending on the region. These values are used to calculate different values of $\epsilon _{\gamma }$ and $f_{\text {jet}}$ (where the number of photon candidates in a given region is altered by the corresponding percentage) which are then used in the jet background estimate. The deviations from the nominal signal region yield (assuming no true photon contamination) are calculated separately for the three final states of $2\gamma$ + 1 jet, $1\gamma$ + 2 jets, and 3 jets, and these values (4, 10 and 21 %, respectively) are taken as systematic uncertainties on the estimates of these backgrounds in the signal region.

An additional uncertainty associated with the data-driven methods employed arises from the choice of kinematic variable used to categorise photons and then to calculate and apply $\epsilon _{\gamma }$ and $f_{\text {jet}}$. The baseline analysis uses three bins in the $p_{\text {T}}$–$\eta$ plane, described in Sect. 5, and separate analyses using either $p_{\text {T}}$-dependence only or $\eta$-dependence only are conducted as well. The largest deviation of the two different methods from the nominal method is 13 %, which is taken as a systematic uncertainty.

6.2 Simulation uncertainties

The uncertainty on the integrated luminosity for the data sample is 2.8 %, derived using the same methodology as that detailed in Ref. [34].

The photon identification efficiency has been directly measured in data using photons from $Z \rightarrow e^{+}e^{-}/\mu ^{+}\mu ^{-}$ radiative decays [12]. The systematic uncertainties on the signal region yield due to the uncertainty on this efficiency measurement are found to range from <1 to 6 % for simulated backgrounds, and from 3 to 7 % for simulated signal processes, depending on the sample.

As mentioned in Sect. 3, the analysis supplements the isolation prescription – $E_{\text {T}}^{\text {iso}}<$ 4 GeV, with a cone size of $\Delta R<$ 0.4 – with an isolation energy correction that is applied to photons with overlapping isolation cones. This procedure improves sensitivity to lower-mass two-photon resonances where the photon pairs are close together in $\Delta R$. To account for possible over- or under-correction due to a photon being near the edge of the isolation annulus, an additional systematic uncertainty is assessed. The $p_{\text {T}}$ values of all isolated photons in simulated samples are calibrated to yield agreement with the values observed in data [35]. Since the calibration factors for isolated photons deviate from one by typically less than 5 %, a value of 5 % is a conservative estimate of the uncertainty on photon $p_{\text {T}}$. To assess the systematic uncertainty on the isolation energy correction, the measured value of the $p_{\text {T}}$ of the other tight photon in the isolation cone is varied by 5 %, the correction procedure is applied, and the effects are propagated to the final event selection in the signal region. For example, using simulated samples of a Higgs boson decaying to four photons via $H \rightarrow aa \rightarrow 4\gamma$, the systematic uncertainty due to this effect is smaller for higher ratios of $m_{a}/m_{H}$ (as large as 6 % when the $p_{\text {T}}$ is varied by $-5$ and <1 % when the $p_{\text {T}}$ is varied by $+5\,\%$, for $m_{H} = 900$ GeV and $m_{a} = 440$ GeV), and the uncertainty is larger for smaller ratios of $m_{a}/m_{H}$ (as large as 69 % when the $p_{\text {T}}$ is varied by $-5$ and 12 % when the $p_{\text {T}}$ is varied by $+5\,\%$, for $m_{H} = 900$ GeV and $m_{a} = 50$ GeV), as the photons tend to overlap within the isolation cone more frequently.

The uncertainties on the event yields due to systematic uncertainties in the photon energy scale and resolution [35] are found to range from <1 to 4 % for the simulated signal and background samples. The uncertainties on the event yields due to systematic uncertainties in the scale factors used to yield agreement between photon identification efficiencies calculated in data and simulated samples [35] are found to range from <1 to 8 % for simulated backgrounds, and from 1 to 4 % for simulated signal processes. The systematic uncertainty on the scaling factor for electrons misidentified as photons in simulated samples is taken to be the statistical uncertainty arising from the calculation, i.e., 4 %, since this is as large as or larger than the systematic uncertainties due to the photon energy scale and resolution, above. The efficiency and uncertainties of the three-photon trigger chain have been determined to be 98.5 ± 0.1 (stat.) ±0.2 % (syst.). The trigger efficiency is calculated using single photons (with $p_{\text {T}}$ values corresponding to the values used for the analysis event selection) from Z boson radiative decays and then, under the assumption that the per-event performance of the photon trigger for one photon is uncorrelated to that for another photon in the same event, multiplying these values to obtain the overall trigger efficiency.

Uncertainties on calculated cross sections for simulated background processes due to QCD renormalisation and factorisation scales and due to the choice of PDF set and value of $\alpha _{\text{ s }}$ used in simulation are addressed via the recommendations of PDF4LHC [36]. The resulting combined uncertainties are found to range from 7 to 16 %, depending on the simulated background process. The total theoretical uncertainty on the SM 3$\gamma$ background process due to the uncertainty on the LO to NLO correction, combined with the uncertainties due to choice of PDF set and renormalisation and factorisation scales, is found to be $\sim$30 %.

Uncertainties exist for the measured or calculated production cross sections for SM particles for which BSM decays are considered as signal scenarios and are accounted for. For the BSM Higgs boson scenario of $h \rightarrow aa \rightarrow 4\gamma$ the gluon fusion production cross section for the SM Higgs boson with $m_{h} =$ 125 GeV [37, 38], $\sigma _{h,\text {SM}} =$ 19.27 pb is used, with an uncertainty of $\pm 10.4$ % due to choice of PDF set and renormalisation and factorisation scales. For the rare decay $Z \rightarrow 3\gamma$, the measured $pp \rightarrow Z$ cross section of $(2.79 \pm 0.02 \pm 0.11)\times 10^{4}$ pb [39] is used. An additional uncertainty of $\pm 12.3$ % – determined by varying the QCD renormalisation and factorisation scales, PDF set, and value of $\alpha _{\text{ s }}$ – is also assessed, to account for variations in the simulation of the kinematics of the final-state photons and, hence, the acceptance in the signal region.

These systematic uncertainties are summarized in Table 2.

Table 2 Systematic uncertainties (%) on the expected event yields in the signal region. The values given for data-driven backgrounds correspond to the three jet backgrounds described in the text. For simulated samples, when a range is given it corresponds to the smallest and largest uncertainties for all simulated backgrounds or signals

7 Results and interpretations

This section presents results based on the number of events in a broad inclusive region and a restricted region focusing on the rare decay $Z\rightarrow 3\gamma$, as well as results from the search for resonances in the di-photon and tri-photon invariant mass spectra.

7.1 Inclusive and $Z\rightarrow 3\gamma$ regions

The number of SM background events expected in the signal region is 1370 ± 140 (combined statistical and systematic uncertainties) and the observed number of events is 1290. The observation is in agreement with the SM expectation. Additionally, while the event selection is optimised for a search for physics beyond the SM as opposed to a measurement of the $3\gamma$ inclusive cross section, the results are nevertheless compatible with the irreducible all-photon process expectations from the SM. The expected and observed yields in the signal region are presented in Table 3, and the expected and observed yields in signal and control regions where all three photons have passed the tight identification criterion are shown in Fig. 2. In the figure, the red hatched band, in the signal region bin, is the combination of statistical and systematic uncertainties on all background sources, while the black hatched band, in the control regions, is the combination of statistical uncertainties of the data-driven jet background estimate and the expected yields from simulated samples of SM background processes. For the inclusive signal region, this corresponds to a model-independent observed (expected $\pm 1\sigma$) upper limit, at the 95 % confidence level (CL), on the number of signal events of 240 (273$^{+83}_{-66}$), and to the model-dependent upper limits on the inclusive fiducial cross section in the aforementioned acceptance for the signal scenarios of the BSM Higgs boson and Higgs boson-like decays and the $Z'$ decays shown in Table 4, where hypothesis testing and limit setting are calculated using the profile likelihood ratio as the test statistic for the $CL_{s}$ technique [40] in the asymptotic approximation [41, 42]. The fiducial efficiencies for each signal scenario are determined with respect to a generator-level kinematic region with the same requirements applied to three-photon events as those used for the analysis signal region. This fiducial region is defined as the set of events that contain three photons where (1) each photon satisfies a pseudorapidity requirement of $|\eta |$ < 2.37, excluding the transition region between the barrel and end-cap of 1.37 $< |\eta |<$ 1.52, (2) the three photons satisfy $p_{\text {T}}>$ 22, 22, and 17 GeV, and (3) each photon satisfies $E_{\text {T}}^{\text {truth iso}}<$ 4 GeV, where $E_{\text {T}}^{\text {truth iso}}$ is a generator-level definition of the photon isolation criterion equivalent to that used for event selection on reconstructed events. The fiducial efficiencies are similar for the considered signal scenarios for mass points where the distributions of photon $p_{\text {T}}$, for all photons, tend to peak higher than $p_{\text {T}}>$ 50 GeV. This is because the overall photon identification efficiency decreases for photons with $p_{\text {T}}<$ 50 GeV [12]. Since the $p_{\text {T}}$ distribution for at least one of the photons for signal scenarios with lower-mass resonances tends to peak at lower values, the fiducial efficiencies are lower.

Fig. 2

Observed and expected yields in signal and control regions for the full mass range (left) and the restricted range of 80 GeV $< m_{3\gamma }<$ 100 GeV (right), for events where all three photon candidates satisfy the tight photon identification criteria. The bins along the horizontal axis correspond to orthogonal subsets of events where each subset is categorised by whether the three photons – ordered from largest to smallest values of $p_{\text {T}}$– passed (“P”) or failed (“F”) the isolation criterion. The leftmost bin is the signal region, composed of events satisfying PPP, and the other bins are the different control regions, where at least one of the photon candidates failed the isolation criterion. The red hatched band, in the signal region bin, is the combination of statistical and systematic uncertainties, while the black hatched bands represent statistical uncertainties. As a result of the data-driven jet background estimate, the statistical uncertainty in each bin is partially correlated with the uncertainty on the data in that bin

Table 3 Expected and observed event yields in the inclusive signal region and for the signal region with a further requirement of 80 GeV $< m_{3\gamma }<$ 100 GeV. Background expectations estimated via simulations are marked sim., whereas data-driven calculations are denoted as D–D. The uncertainties for each row are the combination of statistical and systematic uncertainties for a given background process, and the overall uncertainties in the second to last row are the combined uncertainties for the total background expectations for each signal region

Table 4 Top row observed and expected model-independent upper limits on event yields for new physics processes for the inclusive signal region. Also shown are the efficiencies for the fiducial kinematic region defined in the text for some example mass points for the signal scenarios explicitly considered here, and the corresponding observed and expected ($\pm 1 \sigma$) upper limits on the fiducial cross section within the acceptance. Total statistical uncertainties are quoted for the fiducial efficiencies, and the uncertainties for the upper limits correspond to the uncertainties arising from the $\pm 1 \sigma$ upper limits calculated via hypothesis testing using the combination of statistical and systematic uncertainties

Using the same data-driven and simulation-based methodology restricted to the region 80 GeV $< m_{3\gamma } < 100$ GeV provides a test for the rare decay of the Z boson to three photons. The SM branching ratio for the process is predicted to be $\sim$5$\times$10$^{-10}$ [9], but it has yet to be observed. Table 3 (right) and Fig. 2 (right) summarise the observed counts as well as background expectations in this restricted region. The data are consistent with the SM expectation: 244 events are observed and 233 ± 28 events are expected, while the signal expectation from simulation, for BR$(Z \rightarrow 3\gamma ) = 10^{-5}$ (corresponding to the previous limit from LEP [10]), is 418 ± 9 events. Using the same hypothesis-testing and limit-setting procedure described above, but taking the signal expectation from the simulated sample described in Sect. 4, the observed (expected) limit, at the 95 % CL, on the branching ratio of the Z boson decay to three photons is found to be BR$(Z \rightarrow 3\gamma )<$ 2.2 (2.0) $\times 10^{-6}$, a result five times stronger than the previous LEP limit of $10^{-5}$.

7.2 The $2\gamma$ and $3\gamma$ resonance searches

In addition to the tests based on the number of events in the inclusive signal regions, searches are performed for resonances in the two-photon and three-photon invariant mass ($m_{2\gamma }$ and $m_{3\gamma }$) distributions for events in the inclusive signal region. For these resonance searches, the background contribution is estimated from a fit to the $m_{2\gamma }$ or $m_{3\gamma }$ sideband regions, and thus does not rely upon simulated samples for the background estimate. The sideband is modelled as a fourth-order polynomial, and the size of the sideband is mass-dependent, symmetric around the hypothesised resonance mass, following a local-spectrum approach. The range of the observed mass spectrum that is used for the sideband fit is a local, truncated subset of the full spectrum. For the $m_{2\gamma }$ ($m_{3\gamma }$) resonance search, the sideband is 20 (25) GeV in each direction for $m_{2\gamma }$ ($m_{3\gamma }$) < 90 (230) GeV, where the event counts change rapidly as a function of $m_{2\gamma }$ ($m_{3\gamma }$), and rises to a sideband size of 80 (100) GeV in each direction for $m_{2\gamma }$ ($m_{3\gamma }$) > 195 (425) GeV, increasing roughly linearly with mass as the spectrum becomes smoother. The $m_{2\gamma }$ ($m_{3\gamma }$) resonance search begins at a mass hypothesis of 10 (100) GeV, and proceeds in steps of 0.5 GeV. The signal component of the resonance search is a Gaussian function with a fixed width that varies with particle mass, and the widths are determined from simulated signal samples. Since the simulated signal samples are generated with a narrow-width approximation for both the pseudoscalar a and the $Z'$ in all cases, the $2\gamma$ and $3\gamma$ mass resolutions for this search are equivalent to Gaussian functions that account for detector resolution, and are determined via fits to the simulated signal samples. Hypothesis testing and limit setting are performed using the profile likelihood ratio as the test statistic for the $CL_{s}$ technique in the asymptotic approximation.

The resonance search is performed separately for the three two-photon mass spectra defined by the three possible photon pairings for three photons in the inclusive signal region. As mentioned previously, the photons are ordered by $p_{\text {T}}$, from highest to lowest, and so the three two-photon mass spectra are denoted $m_{12}$, $m_{13}$ and $m_{23}$, where the 1, 2, and 3 refer to the $p_{\text {T}}$-ordered photons. The observed $m_{2\gamma }$ and $m_{3\gamma }$ spectra in the inclusive signal region are shown in Fig. 3. Also shown in Fig. 3, for visualisation purposes only, is the background expectation per bin, determined from the sideband fit to data as a part of the resonance search. The resonance search is performed with a step size of 0.5 GeV and so the final results shown in Figs. 4 and 5 demonstrate sensitivity to resonances with widths appropriate to the BSM models considered here. The widths of the bins in Fig. 3 do not correspond to the mass resolution for the signal scenarios in question. The background estimates and significances shown in Fig. 3 provide a complementary comparison of the local agreement between data and expectation. The lower panels show the significance, in units of standard deviations of a Gaussian function, of the observation in each bin, taking into account the fractional uncertainty on the background as a result of the sideband fit. The significances shown in the lower panels in Fig. 3 are derived from the p value for the background-only hypothesis for each bin, calculated using a frequentist binomial parameter test [43–45]. For regions beyond the sensitivity of the search, no background estimate is shown.

For the H / h $\rightarrow aa \rightarrow 4\gamma$ BSM signal scenario, the photon pairing (among the three $p_{\text {T}}$-ordered photons) that most often corresponds to the photons arising from the decay of the same a particle is (2, 3). As a result, the resonance search in the $m_{23}$ spectrum provides the best sensitivity to this model. The widths of the Gaussian signal component – corresponding to the detector resolution – are taken from simulated samples of a Higgs boson decaying to four photons via a pair of intermediate pseudoscalar a particles, and vary from 0.6 GeV $< \sigma _{\text {Gauss}}<$ 3.2 GeV for 10 GeV $< m_{a}<$ 440 GeV, and are largely independent of $m_{h/H}$.

Fig. 3

Observed spectra of $m_{12}$, $m_{13}$, and $m_{23}$, where the 1, 2, and 3 refer to three $p_{\text {T}}$-ordered photons, as well as $m_{3\gamma }$. For illustration purposes only, also shown is the expected background per bin, determined via unbinned sideband fits to the data as a part of the resonance search, for a hypothesised resonance mass defined by the centre of the bin, as well as the signal expectation for a few mass points for the BSM scenarios considered here. The lower panels show the significance, in units of standard deviations of a Gaussian function, of the observation in each bin, taking into account the fractional uncertainty on the background as a result of the sideband fit. This significance is derived from the p value for the background-only hypothesis for each bin, calculated using a frequentist binomial parameter test [43–45]. The signal distributions used for the $m_{2\gamma }$ resonance searches have two components, a narrow Gaussian core for correctly paired two-photon combinations and a wide distribution for incorrectly paired combinations that is well described by the polynomial used to simultaneously model the background shape for the resonance search described in Sect. 7.2

No excess above background is detected, and upper limits, for a SM-like Higgs boson of $m_{h} = 125$ GeV (and assuming kinematics associated only with gluon fusion SM Higgs boson production), are calculated. Additionally, limits are set for Higgs boson-like scalars with masses larger than 125 GeV. The results of these resonance searches are shown in Fig. 4 for the SM-like Higgs boson of $m_{h} = 125$ GeV (in the top row) and, as an example of a higher scalar mass, for $m_{H} = 600$ GeV (in the bottom row). The resonance search limits for higher values of $m_{H}$ are limited by the small number of events in the mass spectra at higher values of $m_{2\gamma }$. As shown in Fig. 4, the limits vary as a function of two-photon invariant mass, but an overall upper bound on the limits is determined to be $\sigma \times {\text{ BR }}(h \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}$ $1 \times 10^{-3} \sigma _{\text {SM}}$, for 10 GeV $< m_{a}<$ 62 GeV for the SM-like Higgs boson of $m_{h} = 125$ GeV and, for the higher scalar mass case, $\sigma _{H} \times {\text{ BR }}(H \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}<$ 0.02 pb for lower $m_{a}$ values in the range 10 GeV $< m_{a}<$ 90 GeV and < 0.001 pb for higher $m_{a}$ (up to 245 GeV for the resonance search in the $m_{23}$ spectrum for $m_{H} = 600$ GeV shown in Fig. 4). Additionally, using the expected signal yields from simulated samples, inclusive limits are calculated for 300 GeV $< m_{H}<$ 900 GeV and for a range of $m_{a} < m_{H}/2$, including values beyond the range of the mass spectra used for the resonance search. These inclusive limits are shown in Table 5.

Table 5 Expected and observed 95 % CL upper limits on $\sigma _{H} \times {\text{ BR }}(H \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}$. The uncertainties for the expected limits are the $\pm 1 \sigma$ uncertainties resulting from the hypothesis tests for each mass point, taking into account statistical and systematic uncertainties

Fig. 4

Left local p values for the background-only hypothesis as a result of a resonance search with respect to the BSM process $h/H \rightarrow aa \rightarrow 4\gamma$, for $m_{h}$ = 125 GeV (top row) and $m_{H}$ = 600 GeV (bottom row), as a function of $m_{a}$, determined via a search for local excesses in the $m_{23}$ spectrum. Right upper limits, at the 95 % CL, on $(\sigma / \sigma _{\text {SM}}) \times {\text{ BR }}(h \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}$ (top row) and $\sigma _{H} \times {\text{ BR }}(H \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}$ (bottom row). Also shown are the $\pm 1$ and $2 \sigma$ uncertainty bands resulting from the resonance search hypothesis tests, taking into account the statistical and systematic uncertainties from simulated signal samples which are used to determine signal efficiency and Gaussian resonance width due to detector resolution for each mass hypothesis

Moreover, upper limits on the $Z'$ production cross section times the product of branching ratios, $\sigma _{Z'} \times {\text{ BR }}(Z' \rightarrow a+\gamma ) \times {\text{ BR }}(a \rightarrow \gamma \gamma )$, are found to be in the range of 0.04 pb to 0.3 pb, depending upon $m_{Z'}$ and $m_{a}$. Upper limits, at the 95 % CL, on $\sigma _{Z'} \times {\text{ BR }}(Z' \rightarrow a+\gamma ) \times {\text{ BR }}(a \rightarrow \gamma \gamma )$ are shown in Table 6, as a function of $m_{a}$, using the expected signal yields from simulated samples. Additionally, using a narrow-width approximation to the $Z'$ resonance width, local excesses corresponding to Gaussian resonances due to detector resolution are searched for in the $m_{3\gamma }$ spectrum. The Gaussian widths are determined via fits to the $Z'$ simulated signal samples, described in Sect. 4. For the range of $m_{Z'}$ for which the resonance search is possible, the $Z'$ width exhibits a small dependence on $m_{a}$. For each $Z'$ mass point, three different samples are simulated, with different values of $m_{a}$. The average of the three measured $Z'$ widths for each of the $m_{a}$ points simulated is taken as the width for a given $m_{Z'}$, and these values are used for the three-photon resonance search, interpolating for $m_{Z'}$ points between those for which samples are simulated. These values range from 1.5 GeV $< \sigma _{\text {Gauss}}<$ 2.4 GeV for 100 GeV $< m_{Z'}<$ 500 GeV. The results, along with the local p values for the background-only hypothesis, are shown in Fig. 5. The smallest local p value is found to be 0.0003 ($3.4 \sigma$ local significance), at $m_{3\gamma } =$ 212 GeV which, after adjusting for a trials factor [46], corresponds to a global p value of 0.087 ($1.4 \sigma$ global significance).

Table 6 Expected and observed 95 % CL upper limits on $\sigma _{Z'} \times {\text{ BR }}(Z' \rightarrow a+\gamma ) \times {\text{ BR }}(a \rightarrow \gamma \gamma )$. The uncertainties for the expected limits are the $\pm 1 \sigma$ uncertainties resulting from the hypothesis tests for each mass point, taking into account statistical and systematic uncertainties

Fig. 5

Left local p values for the background-only hypothesis as a result of a resonance search with respect to the production of a new vector gauge boson $Z'$ as a function of $m_{Z'}$, determined via a search for local excesses in the $m_{3\gamma }$ spectrum, using a narrow-width approximation to the $Z'$ resonance width. The smallest local p value is found to be 0.0003 ($3.4 \sigma$) which corresponds to a global p value of 0.087 ($1.4 \sigma$). Right upper limits, at the 95 % CL, on $\sigma _{Z'} \times {\text{ BR }}(Z' \rightarrow a+\gamma ) \times {\text{ BR }}(a \rightarrow \gamma \gamma )$. Also shown are the $\pm 1$ and $2 \sigma$ uncertainty bands resulting from the resonance search hypothesis tests, taking into account the statistical and systematic uncertainties from simulated signal samples which are used to determine signal efficiency and Gaussian resonance width due to detector resolution for each mass hypothesis

8 Conclusion

A search for new phenomena in events with at least three photons has been performed using 20.3 fb$^{-1}$ of LHC pp collision data at $\sqrt{s} = 8$ TeV collected with the ATLAS detector at CERN. The SM background expectation is in agreement with the data, and is determined to be 1370 ± 140 events while 1290 events are observed. The model-independent observed (expected) 95 % CL upper limit on the number of signal events is found to be 240 (273$^{+83}_{-66}$). Upper limits at the 95 % CL are calculated on the fiducial cross section $\sigma _{\text {fid}}$ for events from non-SM processes for several signal scenarios. The observed (expected) limit on the branching ratio of the Z boson decay to three photons is found to be BR$(Z \rightarrow 3\gamma )<$ 2.2 (2.0) $\times 10^{-6}$, a result five times stronger than the previous result from LEP.

In addition, a search for local excesses in the two-photon and three-photon invariant mass distributions is conducted. For the two-photon mass spectra, no significant excesses are detected, and the 95 % CL upper limit on $(\sigma / \sigma _{\text {SM}}) \times {\text{ BR }}(h \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}$ (assuming kinematics associated only with gluon fusion SM Higgs boson production) is calculated to vary from $\sim$3$\times$10$^{-4}$ to $\sim$4$\times$10$^{-4}$ for 10 GeV $< m_{a}<$ 62 GeV for a SM-like Higgs boson with a mass of $m_{h} =$ 125 GeV. Limits are set for Higgs boson-like scalars H with masses up to $m_{H} = 900$ GeV and are found to be $\sigma _{H} \times {\text{ BR }}(H \rightarrow aa) \times {\text{ BR }}(a \rightarrow \gamma \gamma )^{2}<$ 0.02–0.001 pb, depending upon $m_{H}$ and $m_{a}$. For the three-photon mass spectrum, the resonance search is conducted in the context of a $Z'$ decaying to three photons. The smallest local p value is found to be 0.0003 ($3.4 \sigma$ local significance), at $m_{Z'} =$ 212 GeV which, after adjusting for a trials factor, corresponds to a global p value of 0.09 ($1.4 \sigma$ global significance). Upper limits at the 95 % CL on the $Z'$ production cross section times the product of branching ratios, $\sigma _{Z'} \times {\text{ BR }}(Z' \rightarrow a+\gamma ) \times {\text{ BR }}(a \rightarrow \gamma \gamma )$, are found to be in the range of 0.04–0.3 pb, depending upon $m_{Z'}$.

These model-independent results are the first of their kind, as are the interpretations for a Higgs boson decaying to four photons via two intermediate pseudoscalar a particles (for a SM-like Higgs boson of $m_{h} =$ 125 GeV and for Higgs-like scalars of higher masses) and for three-photon resonances corresponding to a new vector gauge boson.