Search in diphoton and dielectron final states for displaced production of Higgs or $Z$ bosons with the ATLAS detector in $\sqrt{s} = 13$ TeV $pp$ collisions

A search is presented for displaced production of Higgs bosons or $Z$ bosons, originating from the decay of a neutral long-lived particle (LLP) and reconstructed in the decay modes $H\rightarrow \gamma\gamma$ and $Z\rightarrow ee$. The analysis uses the full Run 2 data set of proton$-$proton collisions delivered by the LHC at an energy of $\sqrt{s}=13$ TeV between 2015 and 2018 and recorded by the ATLAS detector, corresponding to an integrated luminosity of 139 fb$^{-1}$. Exploiting the capabilities of the ATLAS liquid argon calorimeter to precisely measure the arrival times and trajectories of electromagnetic objects, the analysis searches for the signature of pairs of photons or electrons which arise from a common displaced vertex and which arrive after some delay at the calorimeter. The results are interpreted in a gauge-mediated supersymmetry breaking model with pair-produced higgsinos that decay to LLPs, and each LLP subsequently decays into either a Higgs boson or a $Z$ boson. The final state includes at least two particles that escape direct detection, giving rise to missing transverse momentum. No significant excess is observed above the background expectation. The results are used to set upper limits on the cross section for higgsino pair production, up to a $\tilde\chi^0_1$ mass of 369 (704) GeV for decays with 100% branching ratio of $\tilde\chi^0_1$ to Higgs ($Z$) bosons for a $\tilde\chi^0_1$ lifetime of 2 ns. A model-independent limit is also set on the production of pairs of photons or electrons with a significant delay in arrival at the calorimeter.


Introduction
The Standard Model (SM) of particle physics is a renormalizable quantum field theory that provides a framework for understanding fundamental particles and their interactions. Predictions of the SM have been substantiated by experimental results over decades, one highlight being the 2012 discovery of the Higgs boson by the ATLAS and CMS experiments [1,2] at the Large Hadron Collider (LHC) at CERN. However, the SM does not describe gravity, does not contain a dark matter candidate, or provide a solution to the hierarchy problem, pointing to the need for new fundamental physics.
Supersymmetry (SUSY) [3][4][5][6][7][8][9] is a well-motivated theoretical extension to the SM that offers possible answers to many of these questions. The theory predicts the existence of SUSY partners (sparticles) for particles in the SM. Each sparticle has identical quantum numbers to its SM partner, differing only by half a unit of spin. A new symmetry called R-parity acts on supersymmetric fields, and assigns a quantum number of +1 to all SM particles and -1 to sparticles. In R-parity-conserving SUSY models [10][11][12][13][14], the lightest SUSY particle (LSP) is stable. At colliders, sparticles would be produced in pairs, which would then decay in cascades involving other sparticles and SM particles until two final state LSPs are produced.
The weak eigenstates of the SUSY partners of the Higgs and gauge bosons mix to form mass eigenstates called electroweakinos that can be electrically neutral or charged fermions. These are respectively referred to as neutralinos (˜0 1 ,˜0 2 ,˜0 3 ,˜0 4 ) and charginos (˜± 1 ,˜± 2 ), with the subscripts indicating increasing mass. These mass eigenstates are model-dependent combinations of the individual electroweakino degrees of freedom. Of the many new particles predicted in SUSY, existing exclusions from the LHC are generally weaker for these electroweak sparticles as compared to sparticles produced via strong interactions, due in part to the low electroweak process cross sections in collisions [15][16][17][18][19][20].
In gauge-mediated SUSY breaking (GMSB) models [21][22][23][24][25][26], the superpartner of the graviton called the gravitino (˜) is the LSP for typical model parameter values. The weak coupling of the next-to-lightest SUSY particle (NLSP) to the gravitino LSP could generate a lifetime of the NSLP that is non-negligible on the detector length scale, leading to displaced NLSP decay vertices [25]. In GMSB models, the lightest neutralino˜0 1 is often the NLSP. If the combination of weak eigenstates in the˜0 1 mass eigenstate is mostly composed of the supersymmetric Higgs eigenstate (higgsino), then the most likely decay modes will bẽ 0 1 → / +˜.
This search considers GMSB models with the˜± 1 ,˜0 2 and˜0 1 forming an almost degenerate triplet of SUSY partners of the SM electroweak bosons, and with the heavier charginos and neutralinos sufficiently massive to be essentially decoupled. As shown in the example Feynman diagrams in Figure 1, the search focuses on direct pair production of members of the nearly degenerate triplet. Each˜0 2 and˜± 1 that is produced decays to the NLSP plus two light SM fermions (denoted by ), followed by the subsequent NLSP decay via˜0 1 → / +˜. The small mass splitting between electroweakinos results in final state SM fermions with very low , making them challenging to reconstruct and thus not useful for event selection. The Feynman diagrams shown include both → and → final states.
The analysis exploits the precision spatial and timing capabilities of the ATLAS liquid argon (LAr) electromagnetic (EM) calorimeter to achieve sensitivity to the displaced production of the SM Higgs or boson by reconstructing the resultant → or → decays. The main characteristic of such events is the presence of two electromagnetic (EM) objects, either photons or electrons, that originate from the decay of the same LLP parent. The EM objects are reconstructed using only EM calorimeter information, so no attempt is made to separate between the diphoton and dielectron final states. These EM objects are  Figure 1: Feynman diagrams of the signal process considered, targeting pair production of electroweak sparticles decaying to two light SM fermions (including all leptons and quarks except top and bottom) and a˜0 1 particle, which then decays to either a or Higgs boson along with a˜. Each of the˜0 1 particles is required to decay to a Higgs (Z) boson as shown on the left (right) which decays to a diphoton (dielectron) final state. The other˜0 1 is not used in the analysis, and the Higgs/ boson decays with its Standard Model branching ratio.
detected with some delay compared to prompt objects in the final state. Given the size of the ATLAS detector and the LAr timing resolution, the requirement of LAr calorimeter measurements of the delayed EM objects restricts the sensitivity of the analysis to NLSP lifetimes of O(ns). In addition, due to the opening angle between the / boson and the gravitino LSP produced in the NLSP decay, the EM objects would have flight paths that are inconsistent with originating from the primary vertex (PV), and hence are called non-pointing. Hereafter, the daughter particles will be referred to as photons, and the signature as a displaced diphoton vertex (DDV), covering both electron and photon final states. Precise LAr timing and spatial information is used to enhance sensitivity to this signature, and this result represents the first application of a new method to localize a displaced decay vertex position using only LAr measurements.
This analysis utilizes the full Run 2 ATLAS dataset of 13 TeV collisions and is the first LHC search optimized for the DDV signature. Previous ATLAS analyses searched for non-pointing and delayed photons produced in long-lived NLSP decays in the datasets of collisions collected at center-of-mass energies of 13 TeV [27] during Run 2 of the LHC, and at both 7 [28] and 8 TeV [29] during Run 1. No search found an excess above the SM background expectation, and results were interpreted in the context of a particular set of GMSB SUSY models. A recent Run 2 CMS result searching for such models also found data agreeing with the SM prediction [30]. These previous searches provide generic sensitivity to events containing one or more non-prompt photons that do not necessarily originate from a common vertex, making this result the first to exploit the correlation between e/ measurements expected from the GMSB delayed higgsino signal.
The analysis considers a simplified model for signals, where the mass and lifetime of the NLSP are treated as independent parameters. A branching ratio (BR) of unity is considered for the combination of the two NLSP decay modes considered, namely˜0 1 → / +˜, though the relative probability of the two modes is considered a free parameter. The and are assumed to have their Standard Model branching ratios. Both → and → signals are reconstructed as final state photons, as electrons and photons have the same EM shower shape and thus the same EM calorimeter reconstruction. Furthermore, both prompt and displaced track reconstruction efficiencies are low for highly displaced electrons, and the application of calorimeter-based selection can recover sensitivity to these signatures. The photons should have a delay compared to prompt objects, and as = 0 is defined as the expected value for a prompt photon from the PV of the hard collision, the photons are required to have > 0. Measurements of the trajectories of the two photons, as determined by their EM shower shapes, are used to determine a common origin. The separation distance between this secondary vertex candidate and the PV is calculated in the ( , ) plane, and used to categorize the events according to the degree of displacement. This procedure is shown diagrammatically in Figure 2 and described in detail in Section 4.2.
180 6 El ect r omagnet i c bar r el cal or i met er and pr esampl er  Figure 2: Illustration of the two-dimensional calo-vertexing procedure to calculate the R and discriminating variables used in the analysis, with on the y-axis and on the x-axis. The three layers of the LAr calorimeter are highlighted, along with the energy deposits left by the passage of the two daughter photons 1 and 2 of the LLP. The location of the secondary vertex (SV) is determined by the pointing measurements of the two photons. R is defined as the distance in R from the SV to the beamline, and as the distance in from the SV to the PV.
The signal region (SR) is defined by the presence of at least two photons, as well as a high value of missing transverse momentum ( miss T ) resulting from the escaping gravitinos. In addition, a few optimized kinematic selections, designed to enhance the signal-to-background ratio, are applied. The events in the SR were kept blinded until the analysis design was finalized. The background estimation procedure assesses the contribution of SM processes that populate the SR selection and is fully data-driven due to non-Gaussian tails in the timing and vertexing LAr measurements. The predicted background is determined using a control region (CR) with low values of miss T , which is validated using two different orthogonal validation regions (VR). The first VR requires intermediate values of miss T and is denoted VR( miss T ), while the second, denoted VR( ), imposes the same miss T requirement as the SR but reverses the photon timing requirements by requiring < 0 for each photon. By construction the various analysis regions are orthogonal, and any signal contamination in the CR, VR( miss T ), or VR( ) was found to be negligible. The development of the background modeling using the CR, and its validation using the VRs, were finalized before analyzing the data in the SR. A simultaneous likelihood fit of the background model to SR data is performed using the average photon timing distribution in categories of the photon secondary vertex displacement.
The ATLAS trigger and data acquisition system [34] consists of a hardware-based first-level (L1) trigger followed by a software-based high-level trigger (HLT) that reduces the rate of events selected for offline storage to 1 kHz. An extensive software suite [35] is used in the reconstruction and analysis of real and simulated data, monitoring the detector during operation, and the trigger and data acquisition systems of the experiment.

Data and Monte Carlo samples
This search was performed with the full Run 2 LHC dataset, collected by the ATLAS detector between 2015 and 2018. After the application of data quality requirements [36] that ensure good working condition of all detector components, the dataset corresponds to an integrated luminosity of 139.0 ± 2.4 fb −1 [37,38]. The recording of events used for this search was triggered by the presence of two high-T photons, where the full-rate trigger with the lowest available T threshold is used across data-taking years [34,39]. For 2015 and 2016, the trigger used requires the two photons to pass Loose identification (ID) selection, defined in Section 5.1. A trigger based on Medium ID photons became available in 2017 and 2018, and is thus used here for the data collected in those years due to its lower T threshold. Kinematic selections were imposed on the photons to ensure that the selected events lie in the fully efficient regime of the trigger, namely for the leading (subleading) photon T to be greater than 40 (30) GeV.
Monte Carlo event generators were used to simulate the signal targeted by this search. Signal matrix elements were generated at leading order (LO) using MadGraph 2.7.3 [40] with showering and hadronization performed by Pythia 8 [41]. The A14 tune [42] and the NNPDF2.3lo PDF set [43] were used in the event generation. Each event has two long-lived NLSPs and two final state / bosons. Events were filtered such that only one / is required to decay to the desired di-/e resonance, and the other takes its SM branching ratios. Therefore, four separate processes are simulated: → + → SM, → + → SM, → + → SM, and → + → SM. The generated signals are parameterized by the mass of the NLSP, which ranges from 100 GeV to 725 GeV, and its lifetime. The mass of the gravitino LSP is set at 1 MeV, and the lightest chargino and next-to-lightest neutralino have the same mass which is set to be 1 GeV heavier than the lightest neutralino.
For each NLSP mass value, at least two different NLSP lifetimes, typically 2 ns and 10 ns, were simulated. Since the distribution of particle decays follows an exponential decay curve, it is possible to reweight the shape of that curve and thus generalize to other lifetime values. Each event was assigned a weight according to a source signal lifetime, target signal lifetime, and the decay of the event in question. Weights for target lifetimes between 0.25 and 1000 ns were calculated using the generated signal point with the closest lifetime as the source distribution. The signal model MC events were passed through a Geant4 [44] simulation of the ATLAS detector [45] and reconstructed with the same software [35] as used for the data. The generation of the simulated event samples includes the effect of multiple interactions in the same or neighboring bunch crossings (pileup). The effect is assessed with the inclusion of overlaid minimum-bias events, as well as the effect on the detector response due to interactions from bunch crossings before or after the one containing the hard interaction. Events in the simulation were weighted in order to reproduce the amount of pileup observed in the Run 2 data-taking period.
As the background estimation is fully data-driven, no simulation is required of the background processes. Prompt SM → Monte Carlo was used to study the modeling of the specialized e/ variables described in Section 4. These samples were generated using Powheg_v1 [46] interfaced to Pythia 8 with the AZNLOCTEQ6L1 PDF/tune [47,48].

Photon variables
The capability of the ATLAS LAr calorimeter to provide precision spatial and timing information for EM objects is essential for the reconstruction of and sensitivity to DDV events. In addition to standard four-vector information, two key variables are used to characterize photons: the timing of the photon signal, and the pointing of its trajectory back to the beamline. These measurements are almost completely uncorrelated for prompt backgrounds, but the signal is expected to populate higher values of both quantities than the SM expectation, making them excellent variables to discriminate between signal and prompt backgrounds. Since the targeted final state has two photons that share a common secondary vertex, the two pointing measurements are algorithmically combined into two novel vertexing variables, which describe the position of the diphoton vertex in the two-dimensional ( , ) plane. Details of the timing and vertexing calculations and their use in the analysis are provided below. The photons and electrons used in this section are selected with the criteria enumerated in Section 5.1.

Timing
Photons from long-lived NLSP decays reach the LAr calorimeter with a slight delay compared to prompt photons produced directly in the hard scattering. This delay results mostly from the flight time of the heavy NLSP, which would have a relativistic speed ( = / ) that is less than 1. In addition, the opening angle in the NLSP decay, which causes the photon to be non-pointing, results in the geometrical path to the calorimeter being longer than that for a prompt photon from the PV.
The LAr calorimeter has an accordion geometry and excellent timing resolution, with a readout which incorporates fast shaping, and a clock jitter on the readout board that is less than 20 ps. The energy deposited by an EM object in the LAr calorimeter is measured using samples read out from the calorimeter channel at 25 ns intervals. The arrival time is determined using the energy deposit from the second-layer calorimeter cell with the maximum energy deposit among cells in the associated EM cluster ( cell ). For the EM shower of an electron or photon with an energy in the range of interest, this cell typically contains about 20%-50% of the total energy deposited in the EM shower. The energy and timing for each cell are reconstructed by applying the optimal filtering coefficient (OFC) algorithm [49] to four samples of the signal shape.
The time resolution ( ) improves as the cell increases, reaching a lower plateau at O(10) GeV above which the resolution is flat at approximately 200 ps. More specifically, it follows the form ( ) = 0 / cell ⊕ 1 , where cell is the cell energy, ⊕ denotes addition in quadrature, and parameters 0 and 1 are the coefficients of the so-called noise term and constant term, respectively, determined by a fit to data. The time measurements are observed to include a correlated contribution of ≈ 190 ps, which agrees well with the expected spread of times from a single vertex due to their spatial distribution along the beamline.
An offline calibration procedure is necessary to obtain the best possible resolution on each timing measurement. Offline corrections are determined for each interval of validity of the LAr online calibrations, of which there are 13 in the Run 2 dataset. Timing calibration corrections are determined with a dedicated procedure that uses a large sample of → data. This procedure defines a time of zero as the expected measurement from a prompt photon originating at the PV, and includes corrections for offsets between channels, energy dependence, electronic crosstalk, and the position of the PV. The calibration is validated over an independent sample of → events, which also provide a measurement of the expected resolution that is obtained by performing Gaussian fits to the time distributions in bins of cell energy. Application of the calibration constants offline achieves a final resolution of O(100) ps in response to high-energy e/ objects [27].
A source of early and delayed photons in data emerges through satellite bunches of protons that, due to the radio-frequency structure of the LHC accelerator and injection complex, are present in the LHC beams but separated from the main bunches by multiples of ±5 ns. These contribute as a background process to the signal region of interest, while also allowing for an assessment of the OFC reconstruction method in data that more closely matches the expected timing distribution of the signal. The typical population of a satellite bunch is about a factor of one thousand lower than that of the nearby nominal bunch, so collisions between two satellite bunches are suppressed by roughly a factor of a million. Despite their low rate, such satellite collisions are nonetheless observable in the ATLAS data. Figure 3 shows the timing distribution of the electrons with the highest ( 1 ) and second highest ( 2 ) T in data events that are subject to a selection, defined in Section 5.2, that isolates the prompt → process. In addition to the bulk of events clustered around times of zero, satellite collisions are seen with both electron times around ±5 ns and also at +10 ns. Known features of the LHC bunch structure cause a slight asymmetry between the positive and negative populations. The shape of the signal timing distribution is constructed using the timing variable of the simulated signal samples. The photon time resolution is not modeled precisely in the MC simulation, which underestimates the time resolution observed in data. Additional resolution uncertainties are assigned to cover the observed discrepancy between radiative → ℓℓ events in data and those from MC simulation. In addition to a contribution applied independently to each photon, the additional smearing includes a correlated event-level contribution to account for the impact of the spread of the actual time of the collision, which results from the longitudinal profiles of the proton bunches along the LHC beamline. The combined smearing contributions are tuned to match the time performance observed in data using electrons since, due to their similar EM shower developments, electrons have the same timing performance as photons. The correlated and uncorrelated contributions to the time measurement are de-convolved by studying the times of electron-positron pairs in events from the → selection.
In order to exploit the correlation of the timing between two photons produced in the same parent decay, the final analysis variable is the average of the times of the two leading photons ( 1 and 2 ), defined as The distribution of avg that is expected in the signal region as defined in Section 5, for data and several representative simulated signals, is shown in Figure 4. The decay channels have longer timing tails due to the window discussed in Section 5.2. Because the calculation of is dependent on the opening angle between the photons and thus highly displaced signals have underestimated values of the reconstructed , the signals with a lower truth are more likely to be cut away by the lower bound.

Displaced Vertex Finding
The precise spatial resolution and segmentation of the LAr calorimeter provides geometrical information about the origin and direction of travel of EM objects. In contrast to standard displaced vertex finding methods which rely on tracking information from charged particles passing through the inner detector, this analysis uses a novel method called calo-vertexing, where the di-e/ production vertex is localized using only information from the LAr calorimeter. Calo-vertexing is the only way that unconverted photons can be vertexed, while also providing enhanced acceptance for highly displaced electrons that do not have associated tracks.
The calo-vertexing method developed for this search uses well-established pointing variables to localize the DDV in two dimensions. The granularity of the LAr calorimeter in is insufficient for unambiguous localization in this dimension, so the algorithm first rotates all photon measurements onto =0. It then finds each photon's pointing value, defined as the position on the beamline ( = =0) with respect to the PV that its trajectory "points" back to. Pointing is determined using the signals in the first two calorimeter layers, each of which has an associated distance from the beamline and . A line can be drawn through the measurements from the first two calorimeter layers, ( 1 , 1 ) and ( 2 , 2 ), and extrapolated back to the beamline to determine the origin of the photon in . The resolution on the pointing measurement is ≈15 mm for a photon with energy of ≈50-100 GeV in the barrel, with good agreement between simulation TeV, 139 fb s Figure 5: Resolution of the photon pointing variable | origin | as a function of | origin | in units of mm. Shown are prompt → data and simulation in a selection of events with at least two electrons that have 68 GeV < m < 108 GeV, and |Δ ( 1 , 2 )| > 0.1. Overlaid for comparison are representative signals in the SR selection, labeled by the˜0 1 mass in GeV, the˜0 1 lifetime in ns, and the decay channel to or . and data over the observable pointing distribution [27]. Figure 5 shows the resolution of the pointing variable | origin | as a function of | origin | , comparing → -selected data, → simulation, and signal simulation for several benchmark points. The resolution of the predicted signal, → MC, and → data agrees well over the range that is well-populated with events from prompt data, confirming that data and simulation are similarly modeled. Further, the agreement of the signal and → MC across higher pointing values confirms the use of the pointing variable to describe both signal and background processes.
The intersection of the two photon paths as determined by the pointing procedure then defines the location of the reconstructed secondary vertex candidate. The variables R and refer to the location of the vertex along the R and axis respectively, measured with respect to the PV of the event. Figure 2 diagrammatically illustrates an LLP that is produced at the PV and decays after some travel distance into two photons, whose pointing measurements are used to reconstruct the SV location and thus the two vertexing variables.
As the pointing value will tend to be larger for photons that are produced in the decay of an LLP, and the pointing for two photons that share a common vertex will be correlated, the R and measurements are highly discriminating against prompt background. Furthermore, the correlation between R and can be exploited by summing them in quadrature, generating the variable = √︃ 2 R + 2 that is used in the analysis. The distribution of that is expected in the signal region for data and several representative simulated signals is shown in Figure 6. 2 The reconstruction of the invariant mass of the di-e/ system is biased by displacement, as the calculation of the opening angle between the two objects assumes they originate at the beamline. Decays with larger displacement leads to a greater underestimation of the value, leading to a population of events reconstructed below the truth boson mass. As the boson mass is lower than the Higgs boson mass, more mis-reconstructed boson invariant masses are removed by the 60 GeV lower bound.

Event selection
Events are selected based on object quality requirements, event-level features and kinematics, as well as the timing and displacement of the two photon objects. The selection criteria are defined based on optimization procedures that maximize signal sensitivity.

Object Selection
Electrons and photons are reconstructed from calorimeter signals using a dynamical, topological cell clustering-based algorithm [50]. While the signal is reconstructed only with photons, prompt electron objects are used to define the prompt → selection.
Photons are required to satisfy T > 10 GeV and | | < 2.37, excluding the region 1.37 < | | < 1.52 which corresponds to the transition region between the EM barrel and endcap calorimeters. To reduce the background of jets misreconstructed as photons, they are also subject to track-and calorimeter-based isolation requirements. The calorimeter isolation variable is defined as the energy of calorimeter clusters around the photon candidate in the EM calorimeter in a radius of Δ < 0.2. Additional corrections based on the leakage of photon energy outside this window, pileup, and the underlying event contribution are applied [50]. The track-based isolation variable is defined as the scalar sum of the T of all tracks with T > 1 GeV within Δ < 0.2 of the photon candidate. The calorimeter (track) isolation is required to be less than 6.5% (5%) of the photon transverse energy.
A set of photon ID selections uses shower shape variables that describe the energy profiles of the EM showers in the calorimeter, and further enhance the photon efficiency while providing rejection against the background. Three working points are employed in this analysis. Loose ID uses only variables pertaining to the second layer of the EM calorimeter and leakage into the hadronic calorimeter, which minimizes its bias against the identification of non-pointing photons. It is applied to all baseline photon objects, including those used to calculate the miss T , and in an overlap removal procedure that prohibits the double-counting of overlapping objects in an event. The second working point, Medium ID, shares all Loose requirements with the addition of an -dependent selection on a first-layer shower shape variable ratio . While this selection decreases signal efficiency by a few percent across the expected range of displacement values, it also offers a lower T threshold by way of the Medium ID-based diphoton triggers, which was found to enhance the signal-to-background ratio overall. Therefore, photons that are considered as potential signal objects are required to pass Medium ID. The discrepancy between data and simulation related to the Medium identification selection is found to be negligible across displacement values, thus no systematic uncertainty on the efficiency is added. Finally, Tight ID is used to define the photon-enriched background template described in Section 6. It includes a variety of additional selections to ensure good rejection of fake photons. All ID working points are defined in Ref. [50].
Electron candidates must pass the same isolation and identification criteria as the photons (where the electron track is excluded from the track isolation calculation.) The reconstructed track associated to the electron candidate must be consistent with originating at the PV, in that its longitudinal impact parameter 0 and transverse impact parameter 0 must satisfy | 0 · sin | < 0.5 mm and | 0 |/ 0 < 5, respectively. Finally, electron objects are required to have T > 10 GeV and | | < 2.47, excluding the crack region (1.37 < | | < 1.52).
Muons and jets only enter the analysis via overlap removal and contributions to the calculation of the missing transverse momentum miss T . Muons are reconstructed by matching tracks from the inner detector and MS subsystems. Muons without an inner detector track in the range 2.5 < | | < 2.7 but with a MS track that is compatible with the interaction point are also considered. Muon candidates are required to have T > 10 GeV and | | < 2.7, and must satisfy Medium muon identification requirements [51]. Muons are further required to satisfy calorimeter-and track-based isolation requirements [51] that are 95%-97% efficient for muons with T ∈ [10, 60] GeV and 99% efficient for T > 60 GeV. Finally, muon tracks must satisfy | 0 · sin | < 0.5 mm and | 0 |/ 0 < 3. Jets are reconstructed with FastJet [52] using a particle flow algorithm [53] from noise-suppressed positive-energy topological clusters [54,55] in the calorimeter using the anti-algorithm [56] with a radius parameter R = 0.4. Energy deposited in the calorimeter by charged particles is subtracted and replaced by the momenta of tracks that are matched to those topological clusters. The jet four-momentum is corrected for the non-compensating calorimeter response, signal losses due to noise threshold effects, energy lost in non-instrumented regions, and contributions from the pileup. Jets are required to have T > 25 GeV and rapidity | | < 4.4. A jet-vertex-tagger (JVT) multivariate discriminant [57] is applied to jets with T < 60 GeV and | | < 2.4, to suppress jets from the pileup.
An overlap removal procedure is performed in order to avoid double-counting objects, with photons given the highest priority. The procedure is as follows: remove electrons overlapping with photons (Δ < 0.4); remove jets overlapping with photons (Δ < 0.4) and those closely overlapping with electrons (Δ < 0.2); remove electrons "close to" the remaining jets (0.2 < Δ < 0.4); and remove muons overlapping with photons or jets (Δ < 0.4).
The magnitude of the miss T is the absolute value of the negative vector sum of the transverse momenta of the entire event and is calculated via the track-based soft term (TST) approach [58]. It uses selected photon, electron, muon, and jet objects surviving the overlap removal procedure as well as the remaining "soft" tracking terms that were not assigned to any of the remaining physics objects.

Analysis Regions
Selected events are required to have a candidate PV for the hard scatter, reconstructed from at least two charged tracks, each with T > 500 MeV. In case of multiple reconstructed vertices, the PV is selected as the one with the largest sum of the 2 T values of the associated tracks. No photon pointing information is taken into consideration for PV assignment.
As the primary feature of the signal events is the presence of the DDV, at least two photons are required in all analysis events. Therefore, all events are subject to the offline selection of at least two photons that are matched to the triggers described in Section 3. The leading and subleading photon T must be greater than 40 and 30 GeV respectively, and their invariant mass must be > 60 GeV, to ensure the trigger is used in its efficiency plateau. At least one photon must be in the barrel (| | < 1.37), while the other must satisfy | | < 1.37 or 1.52 < | | < 2.37. An event cleaning procedure is applied to reject events from calorimeter noise bursts or other non-collision background, which has a negligible impact on the signal efficiency.
The SR is designed to select events that are consistent with the presence of a DDV. As the DDV selections are sufficient to provide essentially background-free signal regions, the other SM boson produced by the NLSP decay is not used in the analysis in order to mitigate model-dependence of the result. The DDV is assumed to come from the two photons that have the highest momenta in the event, whose pointing values are thus used to compute R and . The two leading photons must both have > 0, for consistency with a delayed signal, as well as < 12 ns, to avoid contamination from adjacent bunch crossings. They are additionally required to have an cell measurement that is read out on High or Medium LAr gain [59], to ensure good performance of the offline calibration. Further, the cell must be at least 5 GeV, a selection that balances signal acceptance with the rejection of very low-energy deposits that contaminate the timing resolution. Selections on 0< R <1500 mm and | |<3740 mm ensure that the DDV is produced within the boundaries of the inner detector. To ensure good resolution on the R and reconstruction, |Δ | between the two photons is required to be greater than 0.1. As the resonance producing the DDV is assumed to be a Higgs or Z boson, < 135 GeV is required. The presence of gravitinos in the final state motivates a selection of miss T > 30 GeV. Finally, two additional selections are imposed in the SR only, on kinematic variables that were found to help isolate the signal, specifically > 70 GeV and the azimuthal angle difference Δ ( 1 , 2 ) < 2.4.
The data-driven background estimation, described in Section 6, is derived from a CR that is orthogonal to the SR by requiring miss T < 20 GeV. Due to correlations between the photon times in the data, the timing shapes in the regions where both photons have opposite-sign timing values are slightly narrower than those in the regions with photons of same-sign times. To ensure that the timing correlations and their impact on the analysis variable shapes match what is expected in the SR, the two photons in CR events are required to have same-sign timing values, though both "positive-positive" and "negative-negative" combinations are allowed in order to increase the CR sample size by a factor of two. The selection of > 135 GeV in the CR is inverted with respect to the SR requirement. Signal contamination in the control region is calculated to be less than 0.1% for all simulated signal points.
Two VRs, which are defined to be orthogonal to both the SR and CR as well as to each other, are used to validate the background prediction. VR( miss T ) is defined in the intermediate region of the miss T spectrum, namely 20 GeV < miss T < 30 GeV, providing a region of kinematic phase space that is close to that of the SR. The inverted selection used in the CR is kept for this region. VR( ) uses the same miss T and selection as the SR, but inverts the timing selection such that only events that have two photons with negative times are allowed. An additional prompt → -enriched region is defined to study the timing and vertexing variables, where events have at least two electrons with an invariant mass that satisfies 68 GeV < m < 108 GeV, and |Δ ( 1 , 2 )| > 0.1. As this region is defined with electron objects and the overlap removal procedure prioritizes photons, this selection is orthogonal to all other analysis regions. Signal contamination in the CR is less than 0.1% across all generated signal mass and lifetime hypotheses, and less than 0.2% in both VRs and the → selection.
As both signal and background have a distinct shape in avg and , a shape fit is performed to exploit the entire spectrum of both variables. The final likelihood fit is performed over the binned avg distribution, for several categories of . Dedicated optimization studies were used to determine the binning in both variables that maximizes signal exclusion significance, while ensuring sufficient background statistics for all bins. The projected significance of a signal hypothesis incorporating all bins was calculated for each binning considered, and optimization was performed separately for each signal point to ensure that all SR phase space is considered. The final binning was chosen based on the discovery reach across the signal grid and is given in Table 1, along with a summary of all region selections discussed above.

Background estimation
SM processes do not produce a DDV with a significant invariant mass. Background events in the SR are therefore a result of processes with either real photons that are misreconstructed to pass the DDV selection criteria (including those from satellite collisions), or other objects that fake photons. Due to non-Gaussian tails in the LAr timing distribution, Monte Carlo generators do not give sufficiently good modeling of the data for this variable. Therefore, a fully data-driven background estimation is used to predict the size and shape of these two background sources in the SR.
The low-miss T CR is used to extract templates of the avg shape from data. It is known [27,29] that the measured pointing and timing distributions of genuine photons are narrower than those of other physics objects (specifically jets) that can be misreconstructed as photons. To capture this shape difference, two templates are defined, one enriched in real photons and the other enriched in fake photons. The real-enhanced photon template is defined by the CR selection given in Table 1, in addition to a requirement that both photons satisfy the Tight identification criteria as discussed in Section 5.1. Similarly, the fake-enhanced photon template is defined as the set of events where at least one photon fails Tight ID.
The background in the high-miss T SR is predicted by transforming and mixing these two templates in the following way. The two templates are scaled to match the observed purity in the SR, derived from the value obtained from data events in the SR. The real photon CR template is multiplied by the fraction of SR events where both photons pass Tight ( ), and the fake photon template is weighted by 1 − . The observed in the SR ranges from approximately 0.8 in the lowest bin to approximately 0.5 in the highest bin. This is done separately in each category. Prior to purity scaling, the largest deviation in purity across the analysis regions is on the order of 10%, and after scaling the distributions of all regions agree within error. Due to the dependence of calorimeter performance on the magnitude of the energy deposit, the photon timing resolution is correlated with its kinematics, specifically its energy. A reweighting procedure is next applied to the templates that match each photon's cell distribution to that of the SR, which mitigates this correlation and ensures the validity of extrapolating over miss T in the background estimation. cell -reweighting is done inclusively in timing and exclusively in categories, separately for each template. The cell shape differences between the CR and SR before reweighting are on the order of 20%, and after reweighting the distributions agree within error. Finally, a mean correction is imposed on all data distributions to correct for small non-zero means (O(ps)) introduced by the application of timing corrections derived from electrons to photon objects. The correction is derived and applied separately in each cell -reweighting region. This procedure requires two features of the SR data events to be unblinded, namely the fraction and the cell distribution. As neither of these variables discriminates between signal and background in the region of interest and the signal contamination in the CR is below 0.1%, this procedure doesn't affect the potential detection of signal events in the SR. The final background template after these transformations is shown in Figures 4 and 6 as a function of avg and , respectively.

Systematic uncertainties
While the sensitivity of this search is dominated by the statistical uncertainties in the dataset, several sources of systematic uncertainty are considered. These are modeled via dedicated studies and accounted for as nuisance parameters in the final statistical treatment, as described in Section 8. The systematic uncertainties can be either evenly distributed variations that affect only the signal normalization, or uncertainties that impact the distribution shape of the signal and/or the background. Systematic uncertainties are not needed on the background normalization in each category, since these normalizations are implemented as nuisance parameters that are free to float in the fit. As expected given the limited statistical precision of the analysis region, the expected sensitivity with and without systematic uncertainties agrees to within 1% for all mass and lifetime signal points considered.

Background Shape Uncertainties
Background shape systematic uncertainties are related to the data-driven background estimation, and are assessed through variations on the estimation procedure.
An uncertainty exists on the purity fraction used to weight the real and fake templates to construct an SR prediction. Since the selection criteria of an analysis region require both photons to either pass or fail Tight ID selection, a binomial error is assigned to the photon purity fraction. This photon purity fraction error is then used to define maximum deviations in the up and down directions of the photon purity fraction in bins of . These range from <1% for small avg and values, and are < 10% at their largest values in the higher bins. A similar uncertainty exists on the cell -reweighting. The reweights are calculated from the ratio of the two-dimensional histograms between the target and the source photon cell spectra. Therefore, the error on the reweighting is the statistical error of the ratio of the two histograms, which define the envelope of variation considered. This effect ranges from O (1)% in the smallest avg and bins to O (10)% for the highest.
Finally, a systematic uncertainty is needed to account for the extrapolation from a low-miss T to a high-miss T pileup profile. This is necessary as the timing resolution for photon objects degrades with increasing pileup, due to the increased contribution of jets from pileup interactions. To estimate this effect, the up and down variations are constructed from CR photons based on whether the average number of simultaneous interactions is less than or greater than a specified -dependent threshold value. These threshold values are all approximately 33, which is the average value of pileup for Run 2, and they agree within 10% across the categories. Given the limited background statistics, a smoothing procedure is applied to the pileup variation in the highest category by merging neighboring bins to reduce statistical fluctuations. The size of this uncertainty is ∼10% across all five categories.

Signal Uncertainties
Both normalization and shape uncertainties on the signal prediction are considered. These arise from experimental conditions, theoretical modeling, or analysis methodologies.
The instrumental systematic uncertainties which affect the signal yield are uncertainties on the photon reconstruction efficiencies, reconstruction of miss T , trigger object matching, pileup reweighting, and the integrated luminosity measurement. The uncertainty in the integrated luminosity is ±1.7%, evaluated using the methodology described in Ref. [37]. The total uncertainty on the signal normalization comes from these instrumental effects and is calculated as the quadrature sum of the detector uncertainties, yielding a value of 2.5-3% depending on the signal model. The impact of these normalization systematic uncertainties on the CL s value, with respect to the fit with all background shape systematic uncertainties, is less than 0.2%.
Theoretical uncertainties exist on the choice of the strong coupling constant ( ), the renormalization and the factorization scale variations, as well as the choice of PDFs. Uncertainties are evaluated by comparing the avg distribution after varying each of the parameters by a factor of two and taking an envelope of the resulting distributions. An overall ±12.6% normalization uncertainty is applied to all signals as a result of these theoretical modeling effects.
Additional uncertainties are considered that affect the distribution of relative event yields across avg bins, namely on the timing resolutions in data as compared to simulation. The impact is measured by implementing an alternate timing smearing, adopting the same procedure of previous non-pointing photon analyses [27]. This procedure has a negligible effect on the results.

Statistical analysis
The compatibility of the data with the background-only hypothesis is evaluated with a binned profile likelihood fit of the avg distribution of the background prediction to that observed in the data, performed simultaneously across all five categories. Similarly, to evaluate the compatibility of the data under the signal presence assumption, a fit of both signal and background templates to data is performed.
The likelihood function is constructed as a product of Poisson probabilities P for each bin of the average timing distribution for each category. It depends on the signal strength , as well as the set of all nuisance parameters (NPs) with Gaussian constraints. The parameter multiplies the expected signal cross-section and is fully correlated across categories. The normalization factor for the background in each category is modeled as a floating NP that is uncorrelated to other bins.
A profile likelihood ratio was used to perform frequentist hypothesis tests. The 0 -value, defined as the probability of statistical fluctuations making the background distributions appear to contain at least as much signal as the data distributions, was calculated for different signal models. The statistical fit procedure, as well as the background estimation as derived from the CR, was validated through background-only fits to the data in the two VRs, which each have negligible signal contamination. Figure 7 shows the post-fit avg distribution of data and background in the five categories for VR( ), showing good agreement between data and background estimation.
The validity of the fit setup and background model was additionally assessed with generated pseudo-datasets. 500 pseudo-datasets were generated by drawing events from the CR, and applying the background estimation transformations described in Section 6 to match the pseudo-dataset timing shape to that of the targeted VR. The number of events for each pseudo-dataset was fixed to the expected yield in the signal region. This process was performed separately for both VR( ) and VR( miss T ), and a background-only fit was performed for each of the 500 pseudo-datasets to obtain a distribution of 0 values for the background-only hypothesis. The mean of the obtained 0 distribution was 0.53 and 0.54 for VR( ) and VR( miss T ) respectively, reflecting good behavior of the background model in each fit.
Further validation was achieved via signal-plus-background fits for spurious signal and signal injection tests. The spurious signal was studied via fits of the background expectation to data in the signal-depleted VRs, ensuring that the fit does not find a signal if none is present. A conservative upper bound on the spurious signal was obtained using the signal point in the grid to which the analysis is most sensitive, which is˜0 1 (135 GeV, 2 ns). The maximum fitted significance of the aforementioned signal point across all three VRs is 0.8 . In the signal injection tests, the fit was performed to a pseudo-dataset consisting of CR events transformed to match the SR timing distribution shape, to which a signal template was added with varying significance. The fitted signal is compared to the injected and they are found to be the same within statistical error across the signal grid. The relatively similar shape of the various signal points in both avg and leads to similar signal-plus-background results across the grid, which also serves as validation that the systematic uncertainties are appropriate for the fit.   Figure 8 shows the avg distributions of SR data and predicted background as determined by a backgroundonly fit with all systematic uncertainties. The observed data are generally found to be in good agreement with the expected background. The largest deviation occurs in the highest category and highest avg bin, where a single event is observed with = 560 mm, a leading photon time of 5.82 ns, and a subleading photon time of 0.45 ns. Such an event is highly incompatible with the signal hypothesis of this search, which should have similarly delayed times for both photons.

Results
The timing of the leading photon around 5 ns and the subleading photon closer to 0 ns makes this event a likely satellite collision candidate, with the leading photon coming from a satellite collision and the second photon arising from an overlaid in-time collision. As a known source of background, events with photons arising from satellite collisions are covered by the data-driven background estimation described in Section 6, and are thus included in the prediction of the event yield in the SR. A simple estimate for the specific contribution of satellite collisions to the total predicted background in the SR is obtained by determining the number of satellite events present in the CR, defined by requiring photon times to be between 4.5 and 6.5 ns inclusively in . Since the presence of satellite collisions is uncorrelated to photon kinematics, a prediction for the SR can be calculated by scaling this number by the ratio of events in the SR to CR. This procedure predicts that 0.5 ± 0.3 satellite events should be observed in the -inclusive SR, where no systematic uncertainties are considered, thus providing further context for the likelihood of this single event observation.
As no significant excess above the background prediction is observed in the SR data, signal-plus-background fits are performed to set upper limits at 95% confidence on the signal production cross-section via the CL s technique [60,61] under the asymptotic approximation [62]. Limits are presented as a function of the two parameters of the signal grid, the˜0 1 mass and lifetime, under the assumption of 100% BR (B) to either Higgs or Z bosons, as well as the BR of the˜0 1 to the SM Higgs boson for specific mass and lifetime hypotheses. For the signal points with insufficient statistical precision where the asymptotic approximation breaks down, no limit is provided. Limits on the cross-section ( →˜0 1˜0 1 ) in fb are shown in Figure 9, as a function of both˜0 1 mass and lifetime, for both and decay modes. Figure 10 shows the 95% CL exclusion for the signal hypothesis in the two-dimensional (m˜0 1 ,˜0 1 ) plane. In these results, the limits for → decaying signals are more stringent than the → decays, as the BR of the boson to the dielectron final state is higher than that of the Higgs boson to two photons, so more events are able to pass the analysis selection criteria. Limits as a function of the˜0 1 BR can also be computed by combining samples with both Higgs and boson decays of the˜0 1 , as well as the decay of the other˜0 1 in the event that does not produce the DDV. In these interpretations, the sum of the˜0 1 BR to the Higgs and bosons is assumed to be 1. Figure 11 shows an example of these limits as a function of B(˜0 1 → +˜) for several mass hypotheses with a fixed 0 1 lifetime of 2 ns.
In addition to the signal cross-section limits, an additional test is performed using only the final timing bin ( avg > 0.9 ns) and no vertexing categorization. This region enables a less model-dependent search for generic DDV signatures in data. The results of this model-agnostic test are summarized in Table 2  To provide some indication of the variations in signal yield and shape, three signal models are shown for each of thẽ 0 1 decay modes, namely˜0 1 → +˜and˜0 1 → +˜. The models shown include a rather low˜0 1 mass value of 135 GeV for lifetimes of either 2 ns or 10 ns, and a higher˜0 1 mass value which is near the 95% CL exclusion limit for each decay mode for a lifetime of 2 ns. Each signal model is shown with the signal normalization corresponding to a BR value of unity for the decay mode in question.  events under the background-only hypothesis ( 0 -value) and its significance (Z) are calculated. Based on the observed and expected number of events, a 95% CL upper limit on the number of excluded events (N excl ) under the background-only hypothesis is calculated.

Conclusion
The first search for displaced diphoton vertices originating from the decay of a massive LLP is presented. The data set used was recorded by the ATLAS detector at the LHC and corresponds to an integrated luminosity of 139 fb −1 of collisions at a center-of-mass energy of √ = 13 TeV. Precise measurements from the ATLAS LAr calorimeter are used to select these events based on the delayed timing and displaced vertices of the final state photons. No significant deviations are observed in the data with respect to the predicted background. One event is observed in the highest timing bin and highest vertexing category, whose features are consistent with a source of expected background from satellite collisions. Results Observed Expected Expected ± 1 Figure 10: The 95% CL exclusion limits on the target signal hypothesis, for˜0 1 lifetime in ns as a function of˜0 1 mass in GeV. The overlaid curves correspond to different decay hypotheses, where the assumed cross-section is for higgsino production, and the˜0 1 decays to +˜or +˜such that B ( +˜) + B ( +˜) = 100%. The curve shown in red represents the decay hypothesis where the˜0 1 decays to +˜with 100% branching ratio. The curve shown in blue represents the decay hypothesis where the˜0 1 decays to +˜with 100% branching ratio.   [55] ATLAS Collaboration, Properties of jets and inputs to jet reconstruction and calibration with the ATLAS detector using proton-proton collisions at √ = 13 TeV, ATL-PHYS-PUB-2015-036, 2015, url: https://cds.cern.ch/record/2044564.    The ATLAS Collaboration