Searches for the $Z\gamma$ decay mode of the Higgs boson and for new high-mass resonances in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

This article presents searches for the $Z\gamma$ decay of the Higgs boson and for narrow high-mass resonances decaying to $Z\gamma$, exploiting $Z$ boson decays to pairs of electrons or muons. The data analysis uses 36.1 fb$^{-1}$ of $pp$ collisions at $\sqrt{s} = 13$ TeV recorded by the ATLAS detector at the CERN Large Hadron Collider. The data are found to be consistent with the expected Standard Model background. The observed (expected - assuming Standard Model $pp\to H\to Z\gamma$ production and decay) upper limit on the production cross section times the branching ratio for $pp\to H\to Z\gamma$ is 6.6 (5.2) times the Standard Model prediction at the 95% confidence level for a Higgs boson mass of 125.09 GeV. In addition, upper limits are set on the production cross section times the branching ratio as a function of the mass of a narrow resonance between 250 GeV and 2.4 TeV, assuming spin-0 resonances produced via gluon-gluon fusion, and spin-2 resonances produced via gluon-gluon or quark-antiquark initial states. For high-mass spin-0 resonances, the observed (expected) limits vary between 88 fb (61 fb) and 2.8 fb (2.7 fb) for the mass range from 250 GeV to 2.4 TeV at the 95% confidence level.

The data analysis uses 36.1 fb −1 of pp collisions at √ s = 13 TeV recorded by the ATLAS detector at the CERN Large Hadron Collider. The data are found to be consistent with the expected Standard Model background. The observed (expected -assuming Standard Model pp → H → Zγ production and decay) upper limit on the production cross section times the branching ratio for pp → H → Zγ is 6.6 (5.2) times the Standard Model prediction at the 95% confidence level for a Higgs boson mass of 125.09 GeV. In addition, upper limits are set on the production cross section times the branching ratio as a function of the mass of a narrow resonance between 250 GeV and 2.4 TeV, assuming spin-0 resonances produced via gluon-gluon fusion, and spin-2 resonances produced via gluon-gluon or quark-antiquark initial states. For high-mass spin-0 resonances, the observed (expected) limits vary between 88 fb (61 fb) and 2.8 fb (2.7 fb) for the mass range from 250 GeV to 2.4 TeV at the 95% confidence level.

Introduction
Since the observation of a Higgs boson (H) by the ATLAS and CMS collaborations [1, 2], its properties have been measured and presented in subsequent publications. Its mass was determined to be m H = 125.09 ± 0.21(stat) ± 0.11(syst) GeV [3]. The couplings to Standard Model (SM) elementary particles were measured and confirmed to be consistent with the predictions for a SM Higgs boson within the present uncertainties [4-6], and alternative spin and CP hypotheses were rejected in favour of the SM hypothesis [7-10].
In the SM, the Zγ decay proceeds through loop diagrams similar to the H → γγ decay. The branching ratio for the Higgs boson decay to Zγ is predicted by the SM to be B(H → Zγ) = (1.54 ± 0.09) × 10 −3 at m H = 125.09 GeV, which is comparable with the branching ratio of the Higgs boson decay to a photon pair, B(H → γγ) = (2.27 ± 0.05) × 10 −3 . However, the number of reconstructed events is significantly smaller if Z boson decays into electron or muon pairs are considered (B(Z → ee) = (3.363 ± 0.004)% and B(Z → µµ) = (3.366 ± 0.007)% [11]). A H → Zγ branching ratio different from the SM prediction would be expected if H were a neutral scalar of different origin [12,13], or a composite state [14], or in models with additional colourless charged scalars, leptons or vector bosons coupled to the Higgs boson and exchanged in the H → Zγ loop [15][16][17].
The H → Z(→ )γ decay ( = e or µ) has been searched for by the ATLAS and CMS collaborations using the full data sets collected at √ s = 7 and 8 TeV [18,19]. No significant excess over the expected background was observed. The ATLAS Collaboration reported an observed (expected) upper limit at the 95% confidence level (CL) of 11 (9) times the SM prediction for a Higgs boson mass of m H = 125.5 GeV. The observed (expected) limit from the CMS Collaboration is 9.5 (10) times the SM prediction for a 125 GeV Higgs boson.
Many models of physics beyond the SM introduce new heavy bosons (X) through extensions of the Higgs sector or as additional gauge fields. Searches for heavy Zγ resonances were carried out by the D0 Collaboration at the Tevatron and by the ATLAS and CMS collaborations. The D0 Collaboration set upper limits on σ(pp → X) · B(X → Zγ) ranging from 2.5 (3.1) pb for a scalar (vector) X mass of 140 GeV to 0.19 (0.20) pb for a mass of 600 GeV [20] using about 1 fb −1 of pp collision data recorded at √ s = 1.96 TeV. The ATLAS Collaboration excluded technimesons decaying to Zγ for technimeson masses between 200 and 700 GeV and between 750 and 890 GeV, and composite spin-0 resonances decaying to Zγ for most of the resonance mass range between 200 GeV and 1.18 TeV for certain model parameters using pp collision data recorded at √ s = 8 TeV [21]. Using 3.2 fb −1 of pp collision data recorded at √ s = 13 TeV, the ATLAS Collaboration set limits on σ(pp → X) · B(X → Zγ) between 295 fb and 8.2 fb for spin-0 resonances produced in gluon-gluon fusion for a X mass range from 250 GeV to 2.75 TeV using leptonic and hadronic Z boson decays [22]. The CMS Collaboration set upper limits on σ(pp → X) · B(X → Zγ) between about 300 fb and about 2.5 fb for spin-0 resonances produced in gluon-gluon fusion for a mass range of 200 GeV to 3 TeV using leptonic and hadronic Z decays and pp collision data taken at √ s = 8 and 13 TeV [23].
This article presents improved searches for decays of the Higgs boson to Zγ as well as for narrow highmass resonances decaying to Zγ using Z boson decays to electron or muon pairs. The Z(→ )γ final state can be reconstructed completely and with high efficiency, good invariant mass resolution, and relatively small backgrounds. It is therefore a powerful experimental signature. The main background is the non-resonant production of Z bosons in conjunction with photons, which is modelled using samples of simulated events. A smaller contribution arises from Z boson production together with hadronic jets when a jet is misidentified as a photon. This background is studied with a dedicated selection of the photon candidate. Both searches are based on pp collision data recorded at √ s = 13 TeV with the ATLAS detector at the LHC during 2015 and 2016, corresponding to a total integrated luminosity of 36.1 fb −1 . The search for decays of the Higgs boson to Zγ benefits from the increased Higgs production cross section due to the increased centre-of-mass energy and also from the larger data set compared to the previous search [18]. In addition, the event categorisation now includes a category sensitive to Higgs boson production via vector-boson fusion. The search is optimised based on the expected production processes and kinematics for a SM Higgs boson. The search for high-mass resonances uses a significantly larger data set than the previous search [22] and is extended to spin-2 resonance production.

ATLAS detector and data sample
The ATLAS detector [24] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4π coverage in solid angle. 1 It consists of an inner tracking detector, electromagnetic and hadronic calorimeters, and a muon spectrometer. The inner detector (ID), immersed in a 2 T axial magnetic field provided by a thin superconducting solenoid, includes silicon pixel and microstrip detectors, which provide precision tracking in the pseudorapidity range |η| < 2.5, and a transition-radiation tracker (TRT) providing additional tracking and information for electron identification for |η| < 2.0. For the √ s = 13 TeV data-taking period, the ID was upgraded with a silicon-pixel insertable B-layer [25], providing additional tracking information from a new layer closest to the interaction point. The solenoid is surrounded by sampling calorimeters: a lead/liquid-argon (LAr) electromagnetic calorimeter covering the region |η| < 3.2, a hadronic calorimeter with a steel/scintillator-tile barrel section for |η| < 1.7 and two copper/LAr endcaps for 1.5 < |η| < 3.2. The forward region is covered by additional coarser-granularity LAr calorimeters up to |η| = 4.9. The calorimeter is surrounded by the muon spectrometer (MS) consisting of three large superconducting toroids each containing eight coils. Precise momentum measurements for muons with pseudorapidity up to |η| = 2.7 are provided by three layers of tracking chambers. The muon spectrometer also includes separate trigger chambers covering |η| up to 2.4.
A two-level trigger system [26] was used during the √ s = 13 TeV data-taking period. The first-level trigger (L1) is implemented in hardware and uses a subset of the detector information. This is followed by a software-based level which runs algorithms similar to the offline reconstruction software, reducing the event rate to approximately 1 kHz from the maximum L1 rate of 100 kHz.
The pp data collected by ATLAS in 2015 and 2016 were taken at a centre-of-mass energy of √ s = 13 TeV and with a bunch spacing of 25 ns. After requiring that the full detector was operational during data-taking and application of requirements on the data quality, the integrated luminosity corresponds to 36.1 fb −1 , of which 3.2 and 32.9 fb −1 were taken during 2015 and 2016, respectively. The average number of pp interactions per bunch crossing (pile-up) ranged from about 13 in 2015 to about 25 in 2016, with a peak instantaneous luminosity of 1.37 × 10 34 cm −2 s −1 achieved in 2016.
The events were collected with triggers requiring either one or two electrons or muons in the event. The single-muon trigger has a transverse momentum (p T ) threshold of 26 GeV and applies a requirement on the muon track isolation. The track isolation is defined as the scalar sum of the transverse momenta of the ID tracks in a cone of ∆R = (∆η) 2 + (∆φ) 2 = 0.3 around the muon. For the trigger used during 2016, the cone size was modified to be ∆R = 10/(p T /GeV) for muons with p T > 33.3 GeV. The track isolation is computed from ID tracks with p T > 1 GeV and with a longitudinal impact parameter z 0 within 6 mm of the z 0 of the muon track, excluding the muon track itself. The track isolation is required to be less than 6% (7%) of the muon's transverse momentum in the 2015 (2016) data set. A second single-muon trigger with a p T threshold of 50 GeV has no requirement on the track isolation. The dimuon trigger has p T thresholds of 22 GeV and 8 GeV and does not apply track isolation criteria. Single-electron triggers with p T thresholds at 24 GeV (26 GeV), 60 GeV, and 120 GeV (140 GeV) are used, as well as a dielectron trigger with a p T threshold of 12 GeV (17 GeV) during the 2015 (2016) data taking. In the 2016 data-taking period, the lowest-threshold single-electron trigger required the track isolation in a cone of ∆R = 0.2 for p T < 50 GeV and ∆R = 10/(p T /GeV) for p T > 50 GeV to be less than 10% of the electron's transverse energy. For all electron triggers, electron candidates must satisfy identification criteria based on the properties of the energy cluster in the electromagnetic calorimeter and its associated track. The singleelectron triggers with lower thresholds use tighter criteria for the electron identification. For H → Zγ events that pass the analysis preselection (see Section 4), the efficiency to pass the trigger selection is z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). The transverse energy is defined as E T = E sin(θ). 92.9% (96.9%) for Z boson decays to muon (electron) pairs. For a high-mass resonance at 1 TeV, the corresponding efficiencies are 94.3% and 99.8% for muon and electron final states, respectively.

Simulation samples
Samples of simulated Monte Carlo (MC) events are used to optimise the search strategy, evaluate the selection efficiency and to study the different background contributions. The generated event samples were processed with the detailed ATLAS detector simulation [27] based on Geant4 [28] (one exception is noted below).
For the H → Zγ search, the mass of the Higgs boson is chosen to be m H = 125 GeV and the corresponding width is Γ H = 4.1 MeV [29]. The SM Higgs boson production was simulated with Powheg Box v2 [30][31][32] using the combined parton distribution function (PDF) set PDF4LHC following the recommendations [33] based on the CT14 [34], MMHT14 [35] and NNPDF3.0 [36] PDF sets and the Hessian reduction method [37][38][39]. The techniques used for the simulation of gluon-gluon fusion (ggF) production, vector-boson fusion (VBF) production, and production in association with a vector boson (WH and ZH, together referred to as V H) and the perturbative order achieved are summarised in Table 1 and detailed below. Higgs boson production in association with a tt pair and other production processes are not considered as their contributions to the total Higgs production cross section are of the order of 0.1% or less. Table 1: Higgs boson production processes produced with Powheg Box with the techniques used and their precision in α s for the event generation (gen.). The total cross section is known with higher precision in QCD and electroweak (norm.) than available in the event generation. The events were reweighted to reproduce the more precise total cross section. Higgs boson production via ggF was simulated with Powheg Box, using the MiNLO approach [40], which achieves next-to-leading-order (NLO) precision for both the inclusive and the H+1-jet process in quantum chromodynamics (QCD). In addition, the NNLOPS approach [41] was used to improve the precision for inclusive observables to next-to-next-to-leading-order (NNLO) in QCD: the Higgs transverse momentum spectrum achieved by this technique was found to be in agreement with the result obtained using QCD resummation with next-to-next-to-leading logarithmic (NNLL)+NNLO precision from the HqT calculation [42,43]. Top and bottom quark mass effects are included up to NLO precision in QCD. The central scale choice for the nominal factorisation (µ F ) and renormalisation scales (µ R ) is µ F = µ R = m H /2. The events were reweighted to reproduce the inclusive cross section at next-to-next-to-next-to-leading-order (NNNLO) precision in QCD and NLO precision in electroweak corrections [29,33,[44][45][46][47].

Process
Higgs boson production via VBF was simulated with Powheg Box at NLO precision in QCD [48]. The events were reweighted to reproduce the inclusive cross section with approximate-NNLO precision in QCD and NLO precision in electroweak corrections [29,33,[49][50][51].
Higgs boson production in association with a vector boson via quark-antiquark initial states was simulated at NLO precision in QCD for inclusive events and H + 1-jet events using the MiNLO technique [52].
The events were reweighted to reproduce the total V H production cross section, including also production via gluon-gluon initial states, at NNLO precision in QCD with NLO electroweak corrections [29,33,[53][54][55]. Three additional event samples of gluon-gluon fusion production are used for the studies of theoretical uncertainties. The first is an event sample generated with MadGraph5_aMC@NLO version 5.2.3.3 with FxFx multijet merging [61,62] of H + 0-jet and H + 1-jet at NLO precision in QCD, using the NNPDF 3.0 PDF set. The decay of the Higgs boson and the parton showering, hadronisation and MPI were provided by Pythia 8.186 using the A14 set of parameters [63] and the NNPDF2.3 PDF set [64]. Two more samples were simulated with Powheg Box v1 [30][31][32]65] with the CT10 PDF set [66] and Pythia 8.186 for parton showering and hadronisation using the AZNLO set of parameters and the CTEQ6L1 PDF set. The two samples were produced with and without MPI to study the uncertainties in the signal acceptance related to the modelling of non-perturbative effects.
Production of CP-even, high-mass spin-0 resonances X in the mass range m X ∈ [300-2500] GeV was simulated for the gluon-gluon fusion and vector-boson fusion production processes and for an intrinsic resonance width of 4 MeV, which is much smaller than the experimental resolution (see Section 5) and referred to as narrow width assumption (NWA). Due to the assumed narrow width of the resonance, the interference between the resonant signal and the non-resonant background is neglected. The ggF (VBF) process was simulated for m X = 300, 500, 700, 750, 800, 1000, 1500, 2000 and 2500 GeV (m X = 300, 500, 1000 and 2500 GeV). Both the ggF and VBF processes were produced with Powheg Box v1 with the CT10 PDF set.
The signal shape and the reconstruction and selection efficiency of the studied high-mass resonances are parameterised as a function of m X . The parameterisation allows the extraction of the signal shape and efficiency for any mass point at which no simulation sample is available.
The background mainly originates from non-resonant production of a a Z boson and a prompt photon (Z+γ), with a smaller contribution from production of Z bosons in association with jets (Z+jets), with one jet misidentified as a photon. Z+γ production within the SM is primarily due to radiation of photons from final-state leptons (FSR) or initial-state quarks (ISR). Both SM processes were simulated using the Sherpa generator [68] (version 2.1.1 for Z+γ and version 2.2.0 for Z+jets), and the matrix elements were calculated using the Comix [69] and OpenLoops [70] generators, where Z+γ production was calculated for real emission of up to two partons at leading order (LO) in QCD and merged with the Sherpa parton shower [71] using the ME+PS@LO prescription [72]. The process of Z+jets was calculated for up to two partons at next-to-leading-order (NLO) and four partons at LO and merged with the parton shower using the ME+PS@NLO prescription [73]. For the Z+γ (Z+jets) samples, the CT10 (NNPDF3.0) PDF set was used in conjunction with dedicated parton shower tuning developed by the Sherpa authors. To study the background model in detail, a large sample of Z+γ events was simulated using fast simulation of the calorimeter response [74].
For all event samples, the additional inelastic pp collisions per bunch crossing were simulated with Pythia 8.186 using the A2 set of tuned parameters [75] and the MRSTW2008LO PDF set [76]. The MC events were reweighted to reproduce the distribution of the average number of interactions per bunch crossing observed in the data.
Corrections derived from trigger, identification, reconstruction, and isolation efficiency measurements for electrons and muons, from identification and isolation efficiency measurements for photons, and from selection efficiency measurements for jets are applied to the simulated events to improve the description of the data. Similarly, energy scale and resolution corrections for all simulated objects are also taken into account.

Event selection and categorisation 4.1 Event preselection
Events are required to have at least one primary vertex candidate, determined using the tracks with transverse momentum p T > 400 MeV reconstructed in the ID. The primary vertex candidate with the largest sum of the squared transverse momenta of the associated tracks ( p 2 T ) is considered to be the primary vertex of the interaction of interest. For H/X → Zγ signal events, the selected primary vertex is within 0.3 mm of the true primary interaction vertex for more than 99% of the events.
The H/X → Z(→ )γ candidate events are selected by requiring two same-flavour opposite-charge leptons ( = e, µ) to form a Z boson candidate and at least one photon candidate.
Muon candidates with |η| < 2.5 are reconstructed by combining tracks in the ID with tracks in the MS. To extend the acceptance beyond that of the ID, muon candidates are reconstructed from tracks reconstructed only in the MS up to |η| = 2.7 [77]. Muon candidates are required to satisfy the medium criterion and have p T > 10 GeV. In order to ensure good track quality, the ID tracks associated with muons in |η| < 2.5 are required to have at least one hit in the silicon pixel detector and at least five hits in the silicon microstrip detector, as well as to extend into the TRT for 0.1 < |η| < 1.9. The muon candidates in 2.5 < |η| < 2.7 are required to have hits in each of the three layers of MS tracking chambers.
Electron candidates are reconstructed from a cluster of energy deposits in neighbouring cells of the EM calorimeter and a track, matched to the cluster, in the ID. They are required to have p T > 10 GeV and be within the fiducial region |η| < 2.47 excluding the candidates in the transition region between the barrel and endcap EM calorimeters, 1.37 < |η| < 1.52. The electrons are identified with the medium likelihoodbased criterion [78] built from the shower shapes of the clusters, the number of hits associated with the track in the ID and the quality of the track-cluster matching.
Both the muon and electron candidates are required to be associated with the primary vertex by requiring the longitudinal impact parameter, ∆z 0 , computed with respect to the primary vertex position along the beam-line, to satisfy |∆z 0 · sin θ| < 0.5 mm, where θ is the polar angle of the track. In addition the significance of the transverse impact parameter d 0 calculated with respect to the measured beam-line position must satisfy |d 0 |/σ d 0 < 3 (5) for muons (electrons) where σ d 0 is the uncertainty in d 0 .
[77]), while the efficiency of the electron identification ranges from about 80% for electrons with p T = 10 GeV to higher than 90% for electrons with p T > 50 GeV (similar to Ref. [78]). The efficiency is typically about 5% higher in the barrel region of the detector than in the endcaps.
The lepton candidates are further required to satisfy additional criteria for track isolation, which is defined similarly to the track isolation used in the trigger (see Section 2), but uses a different track selection and a different cone size in some cases. The track isolation is computed as the scalar sum of the transverse momenta of all tracks in a cone around the lepton candidate with p T > 1 GeV which satisfy loose trackquality criteria and originate from the selected primary vertex, excluding the track associated with the lepton candidate. For muon candidates, the cone size is chosen to be ∆R = 0.3 for p T < 33.3 GeV and ∆R = 10/(p T /GeV) for p T > 33.3 GeV. For electron candidates, the cone size is chosen to be ∆R = 0.2 for p T < 50 GeV and ∆R = 10/(p T /GeV) for p T > 50 GeV. The requirement on the track isolation is chosen such that it is 99% efficient over the full lepton p T range.
An overlap removal procedure is applied to the selected lepton candidates. If two electrons share the same track, or the two electron clusters satisfy |∆η| < 0.075 and |∆φ| < 0.125, then only the highest-p T electron is retained. Electron candidates that are within ∆R = 0.02 of a selected muon candidate are also discarded.
Photon candidates are reconstructed from energy clusters in the electromagnetic calorimeter. Clusters matched to a conversion vertex, reconstructed from either two tracks consistent with a vertex originating from a photon conversion or one track that does not have any hits in the innermost pixel layer and has an electron-like response in the TRT, are reconstructed as converted photon candidates. Clusters without any matching track (clusters with a matching track are reconstructed as electrons as described above) or conversion vertex are reconstructed as unconverted photon candidates [79]. Photon candidates are required to have p T > 10 GeV and |η| < 1.37 or 1.52 < |η| < 2.37. The identification of photon candidates is based on the lateral and longitudinal shape of the electromagnetic shower [79,80]. A loose identification is used for preselection and for background studies.
In order to suppress the events arising from FSR processes and H → * → γ decays, photon candidates within a ∆R = 0.3 cone around a selected electron or muon candidate are rejected.
The selection criteria described in the preceding paragraphs define the event preselection for the leptons and photons included in the reconstruction of the invariant mass of the and γ systems.
The event categorisation described in Section 4.3 used in the search for decays of the Higgs boson to Zγ makes use of hadronic jets produced in association with the Higgs boson candidate. Jets are reconstructed using the anti-k t algorithm [81] with a radius parameter of 0.4 with three-dimensional topological clusters as input [82]. Jets are corrected on an event-by-event basis for soft energy deposits originating from pileup interactions [83] and calibrated using a combination of simulation-and data-driven correction factors accounting for the non-compensating response of the calorimeter and energy loss in inactive regions [84]. Jets are required to have a transverse momentum larger than 25 GeV and |η| < 4.4. To reduce the contamination from jets produced in pile-up interactions, jets with transverse momentum smaller than 60 GeV and contained within the inner detector's acceptance (|η| < 2.4) are required to pass a selection based on the jet vertex tagging algorithm [85], which is 92% efficient for jets originating from the hard interaction.
The jet vertex tagging algorithm is based on the tracks associated with the jet which are consistent with originating from the selected primary vertex. Jet-lepton and jet-photon overlap removal is performed where the jet is removed if the lepton or photon is within a cone of size ∆R = 0.2.

Reconstruction of Z candidates and H/X candidates and final selection
The Z boson candidates are reconstructed from two same-flavour opposite-sign leptons satisfying the preselection criteria and with an invariant mass m larger than 45 GeV. Leptons are required to be consistent with the objects that triggered the event. Trigger efficiency turn-on effects are mitigated by transverse momentum requirements on the leptons that fired the single-lepton or dilepton trigger. If the event was triggered by a single-lepton trigger, the transverse momentum is required to be at least 27 GeV for the leading lepton, and at least 1 GeV higher than the respective trigger threshold in cases where the event was triggered by one of the higher-threshold triggers. If the event was triggered by a dilepton trigger, the transverse momentum is required to be at least 24 GeV (18 GeV) for the leading muon (electron) and 10 GeV (18 GeV) for the subleading muon (electron). For Z → µµ candidates with an invariant mass between 66 and 89 GeV, the invariant mass resolution of the Z boson candidate is improved by correcting the muon momenta for collinear FSR by including any reconstructed electromagnetic cluster with p T above 1.5 GeV lying close to a muon track (with ∆R < 0.15) if the corrected invariant mass is below 100 GeV [86]. A constrained kinematic fit is applied to recompute the four-momenta of the dilepton pair [87] for Z → µµ and Z → ee candidates. The fit models the lepton energy and momentum response as a Gaussian distribution for each lepton, and the Gaussian width is given by the expected resolution. The Z lineshape is used as a constraint with the approximation of the leptons being massless. The lineshape is modelled by a Breit-Wigner distribution. After the application of the FSR corrections and the kinematic fit, Z boson candidates are required to have an invariant mass within 15 GeV of the Z boson mass, m Z = 91.2 GeV [11]. If multiple Z boson candidates pass all requirements, the candidate with the mass closest to the Z boson mass is chosen. About 0.2% (0.5%) of events that pass the final H → Zγ (X → Zγ) selection have more than one Z boson candidate within the 15 GeV mass window.
Higgs boson and X candidates are reconstructed from the Z boson candidate and the highest-p T photon candidate after the preselection.
For the main analyses with the exception of the background studies, the photon candidate used for the reconstruction of the H/X candidate is required to pass the tight identification [79]. The efficiency of the tight identification ranges from 67% (60%) to 90% (95%) for unconverted (converted) isolated photons from p T of 15 GeV to 50 GeV and larger.
Photon candidates are furthermore required to be isolated from additional activity in the detector. A combined requirement on the isolation energy in the calorimeter and the inner detector is used. The calorimeter isolation is computed as the sum of transverse energies of positive-energy topological clusters [82] in the calorimeter within a cone of ∆R = 0.2 centred around the photon shower barycentre. The transverse energy of the photon candidate is removed and the contributions of the underlying event and pile-up are subtracted based on the method suggested in Ref. [88]. The track isolation for a cone size of ∆R = 0.2 is used and for converted photons the tracks associated with the conversion are removed. The calorimeter (track) isolation is required to be less than 6.5% (5%) of the photon p T . The efficiency of the isolation requirement for photons satisfying the tight identification criteria ranges from approximately 60% for p T of 15 GeV to more than 90% for p T of 40 GeV and larger.
For the H → Zγ (X → Zγ) search, the photon transverse momentum requirement is tightened to 15 GeV The invariant mass of the final-state particles, m Zγ , is required to satisfy 115 GeV < m Zγ < 170 GeV for the H → Zγ search and 200 GeV < m Zγ < 2500 GeV for the high-mass resonance search. Figure 1 shows the invariant mass distribution for simulated H → Zγ candidates after the final selection with and without the lepton momentum corrections from the FSR recovery and the kinematic fit. Improvements of the m µµγ resolution of 3% are observed for m H = 125 GeV from the FSR recovery. The kinematic fit improves the m µµγ (m eeγ ) resolution by 7% (13%) at the same mass. For high invariant masses, the m µµγ resolution improvement varies from 10% at m X = 300 GeV to about 50% for m X > 1.5 TeV, while the m eeγ resolution is improved by 9% at m X = 300 GeV and by 3% or less above m X = 500 GeV. The constrained kinematic fit is particularly effective at large m X for the Z → µµ final state due to the decreasing precision of the momentum measurement for increasing muon p T .

Categorisation
Events are split into mutually exclusive event categories that are optimised to improve the sensitivity of both the H → Zγ and X → Zγ searches. The event categories separate events on the basis of the expected signal-to-background ratio and of the expected three-body invariant mass resolution. Different categories are used in the search for decays of the Higgs boson to Zγ and the search for high-mass resonances.
The H → Zγ search uses six exclusive event categories and events are assigned to the categories in the following order: • VBF-enriched : Events are required to have at least two jets. A boosted decision tree (BDT) that was trained to separate VBF events from other Higgs boson production modes and non-Higgs back- grounds is applied. It uses six kinematic variables as input, computed from the Z boson candidate, the photon candidate and the two jets with the largest transverse momenta: -The invariant mass of the two jets (m j j ), -The separation of the jets in pseudorapidity (∆η j j ), -The azimuthal separation of the Zγ and the dijet systems (∆φ Zγ, j j ), -The component of the transverse momentum of the Zγ system that is perpendicular to the difference of the 3-momenta of the Z boson and the photon candidate The smallest ∆R separation between the Z boson or photon candidate and the two jets (∆R min Z/γ, j ), -The difference between the pseudorapidity of the Zγ system and the average pseudorapidity of the two jets (|η Zγ − (η j1 + η j2 )/2|).
The variable p Tt is strongly correlated with the transverse momentum of the Zγ system, but has better experimental resolution [89,90]. Any requirement on ∆φ Zγ, j j effectively vetoes additional jets in the event by restricting the phase space for additional emissions and, to avoid uncontrolled theoretical uncertainties, the BDT does not use shape information for events with ∆φ Zγ, j j > 2.94 by merging these events into one bin. A minimum value of the BDT output (BDT > 0.82) is required. The expected and observed distributions for two input variables, m j j as a typical variable to select events with VBF topology and ∆φ Zγ, j j , which serves as an implicit third-jet veto, are shown in Figure 2 for selected events with at least two jets.
• High relative p T : Events are required to have a high p T photon, p γ T /m Zγ > 0.4. • ee high p Tt : Events are required to have high p Tt (p Tt > 40 GeV) and a Z boson candidate decay to electrons.
• ee low p Tt : Events are required to have low p Tt (p Tt < 40 GeV) and a Z boson candidate decay to electrons.
• µµ high p Tt : Events are required to have high p Tt (p Tt > 40 GeV) and a Z boson candidate decay to muons.
• µµ low p Tt : Events are required to have low p Tt (p Tt < 40 GeV) and a Z boson candidate decay to muons.
For SM H → Z(→ )γ events, the reconstruction and selection efficiency (including kinematic acceptance) is 21.5%. Table 2 shows the expected signal efficiency times acceptance for each of the different SM Higgs boson production processes in each category, as well as the expected relative contribution of a given production process to each category. The VBF-enriched category is expected to be about 68% pure in VBF events. The high relative p T and high p Tt categories are expected to be slightly enriched in VBF and V H events. Overall, about 40 H → Zγ events are expected to be selected.  Figure 2: Kinematic variables used in the BDT used to define the VBF-enriched category: (a) the invariant mass of the two jets with the highest transverse momenta, m j j and (b) the azimuthal separation of the Zγ and the dijet system, ∆φ Zγ, j j for events with at least two jets and 115 GeV < m Zγ < 170 GeV. The observed distribution (normalised to unity) is shown as data points. The contributions from Z + γ events (obtained from simulation) and the contribution from Z+jets (obtained from data control regions described in the text) are shown as stacked histograms. The corresponding expected distributions for Higgs bosons produced via gluon-gluon fusion and vector-boson fusion production for m H = 125 GeV are shown as open histograms. The ∆φ Zγ, j j distribution is shown before the suppression of the shape information for ∆φ Zγ, j j > 2.94. Table 2: The expected signal efficiency times acceptance, denoted by , per production mode for each category after the full event selection, as well as the expected fraction f of each production process relative to the total signal yield, for simulated SM Higgs boson production assuming m H = 125 GeV. The expected number of signal events per production process is also given.  The search for high-mass resonances uses two categories, one for each Z boson candidate decay mode, Z → ee and Z → µµ, to benefit from both the better invariant mass resolution in the electron channel at large m Zγ and the differences in the systematic uncertainties between electrons and muons. The invariant mass resolution, measured by the Gaussian width of the signal model (see Section 5), ranges from 2.8 GeV (3.1 GeV) at m X = 250 GeV to 16 GeV (36 GeV) at m X = 2.4 TeV for Z → ee (Z → µµ).

Signal and background modelling
The signal and background yields are extracted from the data m Zγ distribution by assuming analytical models. The parameters that describe the shape of the signal are obtained from simulated signal samples. The analytical models used for the background shape are chosen using simulated background samples and the values of their free parameters are determined from the fit to data.

Signal modelling
The signal mass distribution in the searches for both the Higgs boson and the high-mass resonance decay to Zγ is well modelled with a double-sided Crystal Ball the simultaneous fit. The parameterisation is done separately for each of the three models considered, a spin-0 resonance and a spin-2 resonance produced via either gluon-gluon or quark-antiquark initial states. Figure 3 shows the MC-simulated m Zγ distribution at m H = 125 GeV for the low p Tt categories and at m X = 1000 GeV. Similar fit qualities are obtained for all the categories in both searches.
Additionally, the signal efficiency defined as the number of events satisfying all the selection criteria (as given in Section 4) normalised to the total number of events is needed to extract σ · B(H/X → Zγ) from the measured yield. For the H → Zγ search, the signal efficiency times the acceptance in each category are shown in Table 2. 2 For the search for high-mass resonances, the signal efficiency is parameterised as a function of the resonance mass with an exponentiated second-order polynomial. Figure 4(a) shows the reconstruction and selection efficiency for X → Z(→ )γ events for a spin-0 resonance produced in gluon-gluon fusion, separately for Z → ee and Z → µµ. The efficiencies range from about 30% to about 46% in the mass range from 250 GeV to 2.4 TeV. For a spin-0 resonance produced via vector-boson fusion, the efficiency is larger by up to 4% over the full resonance mass range considered. Figure 4(b) shows the reconstruction and selection efficiency for spin-2 resonances produced via gluon-gluon and quark-antiquark initial states as a function of the resonance mass. For spin-2 resonances produced in gluon-gluon (quark-antiquark) initial states, the efficiencies range from about 22% (28%) to about 35% (54%) in the mass range from  Figure 4: Reconstruction and selection efficiency (including kinematic acceptance) for the X → Zγ final state as a function of the resonance mass m X (a) for a spin-0 resonance via gluon-gluon fusion, separately for the ee and the µµ categories, and (b) for a spin-2 resonance produced via either the gluon-gluon or the quark-antiquark initial states. The markers show the efficiencies obtained from simulation, while the curves represent the parameterisation used in the analysis. The efficiencies are given with respect to (a) X → Z(→ ee)γ and X → Z(→ µµ)γ, respectively, 250 GeV to 2.4 TeV. The efficiency differences between the spin-0 resonance produced via gluon-gluon fusion, the spin-2 resonance produced via gluon-gluon initial states and the spin-2 resonance produced via quark-antiquark initial states are primarily related to the different photon transverse momentum distributions between the different production mechanisms.

Background modelling
The background is mainly composed of non-resonant production of a Z boson in association with a photon (irreducible background), and of inclusive Z+jets events where a jet is misidentified as a photon (reducible background), and the relative contributions are determined using data as described below. Contributions from other sources, such as tt production, W/Z events, and, for the H → Zγ search, from other Higgs boson decays are expected to be negligible based on studies of simulated events. The background exhibits a smoothly falling distribution as a function of the invariant mass of the candidate Z boson and photon, m Zγ .
The estimated background composition is used to construct simulated background samples with the same composition as the data background. These samples are used in the optimisation of the selection criteria, the choice of analytical model of the background shape, and the estimation of the related systematic uncertainties. The searches rely only indirectly on the measured background composition since the background shape parameters are determined from the data.
The composition of the background is estimated using a combined binned fit to the calorimeter isolation distribution of the photon candidate in the signal region and in a control region enriched in Z+jets background. In the control region, photon candidates are required to fail the tight identification, but to pass a modified loose identification. It differs from the tight identification by removing the requirements on four out of nine shower shape variables which are the least correlated with the calorimeter isolation [93]. The calorimeter isolation distribution for photons and the contribution of true photons to the control region are determined from simulation, while the calorimeter isolation distribution for jets is determined in the fit and assumed to be the same in the signal and control regions. This assumption is supported by extensive studies performed in the context of earlier analyses [93]. The composition is determined in the inclusive selection for the H → Zγ search and the fraction of Z+γ events is found to be 0.838 ± 0.005 (stat.) ± 0.031 (syst.). For the high-mass resonance search, the fraction of Z+γ events is found to be 0.916 ± 0.009 (stat.) +0.013 −0.019 (syst.). The systematic uncertainties are estimated by varying the set of shower shape variables that are removed from the tight identification to define the modified loose identification [93]. The results of the composition estimate are cross-checked with a two-dimensional sideband technique [93] based on the calorimeter isolation of the photon candidate and whether or not the photon candidate satisfies the tight identification criteria (when the photon fails the tight identification it is still required to pass the modified loose identification), which gives consistent results.
The analytical model of the background and the m Zγ range used for the final fit are chosen to limit the bias in the extracted signal yield, while at the same time limiting the number of free parameters in the fit to avoid degradation of the sensitivity [1]. For each category used in either analysis, the bias (also referred to as spurious signal) is estimated as a function of the signal invariant mass by performing a signal+background fit to a m Zγ background-only distribution with small statistical fluctuations. The background-only distribution is constructed from the fast simulation of Z+γ events, and the contribution from Z+jets events is included by reweighting the Z+γ simulated distribution as follows: for each category, the shape of the Z+jets contribution is determined in a data control region defined by requiring that the photon candidate fails to satisfy the identification and isolation criteria. To smooth statistical fluctuations in the Z+jets shape, a first-order polynomial is fitted to the ratio of the Z+jets and Z+γ shapes, and the smoothed Z+jets shape is constructed from the fit result and the Z+γ shape. The reweighting of the Z+γ distribution to take into account the Z+jets contribution is determined from the smoothed Z+jets shape. The normalisation of the Z+jets contribution is determined from the number of events obtained when applying the selection and categorisation to the Z+γ and Z+jets simulation samples and the purity of the inclusive sample, measured as described above. The spurious signal is required to be less than 40% (20%) of the expected statistical uncertainty in the signal yield, which is dominated by the expected statistical uncertainty of the background, in the search for Higgs boson (high-mass resonance) decays to Zγ. The looser requirement for the H → Zγ search is chosen to improve the robustness of the procedure against statistical fluctuations in the simulated Z+γ event sample. If two or more considered functions satisfy this requirement, the function with the fewest number of parameters is chosen.
For the H → Zγ search, the fit range is also optimised on the basis of the spurious signal estimates, taking into account the spurious signal and the number of parameters of the chosen functions in all categories. A fit range from 115 GeV to 150 GeV is selected. A second-order Bernstein polynomial is chosen as the parameterisation for the VBF and high relative p T categories, and a fourth-order Bernstein polynomial is chosen for the other categories. For the chosen parameterisation, the largest spurious signal obtained in a window of 121-129 GeV is assigned as a systematic uncertainty in each category associated with the choice of background function and ranges from 1.7 events in the VBF category to 25 events in the µµ low p Tt category. The choice of background functions is validated by using second-and third-order polynomials for the smoothing of the Z+jets background shape and by varying the Z+γ purity by ±15%. The large variation of the Z+γ purity is chosen to cover the purity differences between the different categories and intended to also account for the additional uncertainty in the estimation of the Z+jets invariant mass distribution.

The high-mass resonance search considers as a background model a class of functions [22] given by
where x = m Zγ / √ s, N is a normalisation factor, k determines the number of terms considered in the exponent, and b and a k are determined by the fit. When testing on the background-only distribution constructed using the simulated Z+γ sample taking into account Z+jets contributions as discussed before, the spurious signal criterion is found to be satisfied for the full mass range for k = 0. The spurious signal used as an estimate of the systematic uncertainty is parameterised as a smooth function of the invariant mass. It ranges from 3.6 (6.1) events at 250 GeV to 0.01 (0.005) events at 2.4 TeV for the Z → µµ (Z → ee) channel. The choice of analytic function for the background shape is validated by using second-and third-order polynomials for the smoothing of the Z+jets background shape, by varying the Z+γ purity by ±5% (motivated by the range of purities estimated in the two categories), and by varying the PDFs in the Z+γ simulation using the uncertainties associated with the different eigenvectors of the PDF set.
The possibility of needing higher-order functions when fitting the selected analytical function to the data m Zγ distribution is further investigated by an F-test. The test statistic F defined as compares the fit qualities between less and more complex functions. The χ 2 0 (χ 2 1 ) is the χ 2 value of a binned fit with the less (more) complex parameterisation, p k is the number of free parameters of each fit, and n is the number of bins of the invariant mass distribution. Should the probability to find values of F more extreme than the one measured on data be less than 5%, the less complex parameterisation would be rejected in favour of the more complex parameterisation. The binning for the F-test is chosen to guarantee a sufficient number of events in each bin. For the H → Zγ search, the test is carried out to determine if there is any indication that a higher-order Bernstein polynomial is required. It does not lead to a change in the chosen parameterisation. For the high-mass search, the test is performed to determine whether or not the quality of the fit to data is improved significantly if using k = 1. The test confirms that the choice of k = 0 is adequate. The χ 2 per degree of freedom is 1.2 for 30 degrees of freedom (1.1 for 17 degrees of freedom) or better for the chosen parameterisations for the H → Zγ (high-mass resonance) search.

Systematic uncertainties
The experimental and theoretical uncertainties that are considered in the searches can be grouped into three classes: uncertainties associated with the parameterisation of the signal and background distributions (see Section 6.1), experimental uncertainties in the efficiency and acceptance affecting the expected event yields (see Section 6.2), and theoretical uncertainties in the modelling of the signal in the simulation (see Section 6.3). The nuisance parameters in the likelihood function (see Section 7) represent the uncertainties which are studied in each category using the simulated signal samples generated at m H = 125 GeV for the H → Zγ search, and m X = [300 − 2500] GeV at high-mass. The main experimental sources of uncertainty are summarised in Table 4.

Uncertainties from signal and background modelling
The uncertainties in the lepton and photon momentum and energy scale and resolution impact the modelling of the signal. Their impact is assessed by comparing the nominal m Zγ shape parameters with the m Zγ shape parameters after varying the lepton and photon momentum and energy scale and resolution by their uncertainties. Uncertainties in both the position (µ CB ) and width (σ CB ) of the signal m Zγ distribution are considered.
The systematic uncertainties in the muon momentum scale and resolution were determined from Z → µµ and J/ψ → µµ events using the techniques described in Ref. [77]. At m H = 125 GeV, the uncertainty in the muon momentum scale (resolution) varies σ CB by up to 0.5% (4.0%). In the high-mass search, the effect of the muon momentum scale (resolution) uncertainty is to change σ CB by up to 3.0% (4.0%). The typical effect of the muon momentum scale uncertainty is to change µ CB by < 0.1% of its nominal value.
The systematic uncertainties in the electron and photon energy scale and resolution follow those in Refs. [94,95]. The overall energy scale factors and their uncertainties were determined using Z → ee events. Compared to Ref.
[95], several systematic uncertainties were re-evaluated with the 13 TeV data, including uncertainties related to the observed LAr cell non-linearity, the material simulation, the intercalibration of the first and second layer of the calorimeter, and the pedestal corrections. At m H = 125 GeV, the uncertainty in the electron and photon energy scale (resolution) results in variation in σ CB between 0.6% and 3.5% (1.1% and 4.0%) depending on the category. For a high-mass resonance, σ CB varies between 1.0% and 4.0% (4.0% and 30%) due to uncertainties in the electron/photon momentum scale (momentum resolution). The variation in µ CB is less than 0.2% (0.6%) at m H = 125 GeV (at high masses).
For the H → Zγ search, an additional uncertainty in the assumed Higgs mass m H = 125.09 GeV is added in the fit, reflecting the 0.24 GeV [3] uncertainty in the measured Higgs boson mass.
The uncertainty due to the choice of background function is taken to be the signal yield (spurious signal) obtained when fitting the m Zγ spectra reconstructed from background-only distributions as discussed in Section 5.

Experimental uncertainties affecting the signal efficiency and acceptance
Experimental uncertainties affecting the signal efficiency and acceptance can be either correlated between all event categories (yield uncertainties) or anticorrelated between some of the categories (migration uncertainties) when they are related to how the signal populates the event categories.
The uncertainty in the combined 2015+2016 integrated luminosity is 3.2%, correlated between all categories. It is derived, following a methodology similar to that detailed in Ref.
[96], from a preliminary calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016.
A variation in the pile-up reweighting of the simulation is included to cover the uncertainty in the ratio of the predicted and measured inelastic cross sections in the fiducial volume defined by m > 13 GeV where m is the mass of the hadronic system [97]. The uncertainty in the signal efficiency is no more than 0.03% (0.2%) for the H → Zγ (high-mass resonance) search.
The uncertainties in the reconstruction, identification, isolation, and trigger efficiency measurements for muons, electrons and photons (see Section 4) are treated as fully correlated between all categories. They are determined from control samples of J/ψ → µµ and Z → µµ for muons, J/ψ → ee and Z → ee for electrons, and Z → γ, Z → ee, and inclusive photons for photons, using methods described in Refs. [77,78,80].
For the H → Zγ search, the uncertainties in the signal efficiency from the photon identification and isolation are found to be no more than 1.7% and 0.4%, respectively. The uncertainties in the signal efficiency from the electron reconstruction, identification, isolation, and trigger are found to be no more than 0.4%, 1.6%, 0.2%, and 0.1%, respectively. The uncertainties in the signal efficiency from the muon selection and trigger are determined to be no more than 1.6% and 3.5%, respectively. For the high-mass search, the uncertainties in the signal efficiency from photon identification and isolation are found to be no more than 2.6% and 0.6%, respectively. The uncertainties in the signal efficiency from the electron reconstruction, identification, isolation, and trigger are found to be no more than 1.0%, 2.6%, 3.5%, and 0.2%, respectively. The uncertainties in the efficiency from the muon selection and trigger are determined to be as large as 0.7% and 4.2%, respectively. The uncertainty due to the limited size of the simulated event samples ranges from 1.2% to 2.0% for the search for high-mass resonances.
In the H → Zγ search, the expected signal yield in the VBF category is affected by the jet energy scale and resolution and the jet vertex tagging efficiency. The corresponding uncertainties are anticorrelated with the other categories. Uncertainties in the jet energy scale and resolution are estimated from the transverse momentum balance in dijet, γ+jet and Z+jet events [84]. Uncertainties in the efficiency of the jet vertex tagging are estimated by shifting the associated corrections applied to the simulation by an amount allowed by the data. The uncertainties in the category acceptances are as large as 4.6%, 6.9%, and 4.8% from the data-driven jet calibration, the impact of the jet flavour composition on the calibration, and the jet vertex tagging.

Theoretical and modelling uncertainties
For the H → Zγ search, theoretical and modelling uncertainties in the SM predictions for Higgs boson production and the decay to the Zγ final state are taken into account and are summarised in Table 5. They fall into two classes: uncertainties in the total predicted cross sections, the predicted decay branching ratio and the total efficiency, correlated between all categories; and uncertainties in the event fractions per category, anticorrelated between certain categories.
Uncertainties related to the total acceptance and efficiency for H → Zγ events affect the extraction of the signal strength, the branching ratio of H → Zγ assuming SM Higgs boson production, as well as the product of the Higgs boson production cross section and the branching ratio of H → Zγ (see Section 7). The uncertainty in the total efficiency due to the modelling of multiple-parton interactions is estimated from the difference in efficiency with and without multiple-particle interactions for the gluon-gluon fusion simulation sample, and found to be 5.3%.
The uncertainties related to the predicted Higgs boson production cross section affect the extraction of the signal strength as well as the branching ratio of H → Zγ assuming SM Higgs boson production. The uncertainties in the predicted total cross sections of the different Higgs boson production processes due to the perturbative order of the calculation and the combined uncertainties in the PDFs and α s are 3.9% and 3.2% for gluon-gluon fusion production, respectively, and range from 0.4% to 3.8% for the other production processes for a Higgs boson mass of 125.09 GeV and a centre-of-mass energy of 13 TeV [29]. An additional 5.0% [98] uncertainty accounts for the effect, in the selected phase space of the γ final state, of the interfering H → γ decay amplitudes that are neglected in the calculation of Ref. [29]. They originate from internal photon conversion in Higgs boson decays to diphotons (H → γ * γ → γ) or from Higgs boson decays to dileptons with an off-shell lepton (H → * → γ) [99, 100].
The uncertainty in the predicted Higgs boson branching ratio to Zγ affects the extraction of the signal strength. The relative theoretical uncertainty in the predicted Higgs boson branching ratio is 5.9% [29].
Uncertainties in the modelling of kinematic distributions in the simulation of Higgs boson production processes affect the predicted event fractions in the different categories. The uncertainty in the modelling of the production of jets in gluon-gluon fusion production due to the perturbative order in QCD is estimated by scale variations in MCFM [101]. It accounts for the uncertainty in the overall normalisation of H + 2jets events as well as the uncertainty due to the use of ∆φ Zγ, j j , which serves to apply an implicit third-jet veto, in the VBF BDT. The estimation of this uncertainty uses an extension of the Stewart-Tackmann method [102,103]. It corresponds to 45% of the ggF contribution to the VBF category. Additional uncertainties are assigned to account for potential mismodelling of the variables that serve as input to the VBF BDT (see Section 4). They are estimated by reweighting the simulated ggF events to match the distributions in m j j , ∆η j j , p Tt , ∆R min Z/γ, j , and |η Zγ − (η j1 + η j2 )/2| obtained from MadGraph5_aMC@NLO. The resulting uncertainty in the ggF contribution to the VBF category is 15%. Uncertainties in the modelling of the Higgs boson p T spectrum are taken to be the envelope of the variations of the renormalisation, factorisation, and resummation scales obtained using HRes 2.3 [104] to simulate the p T spectrum. The resulting uncertainties are evaluated using the Stewart-Tackmann method [102,103] for the high relative p T and the p Tt categorisation and found to range from 8.4% to 22%. Uncertainties from the choice of PDF set and α s are evaluated using the combined error PDF set, which takes into account 30 variations of NNLO PDFs and two variations of α s , following the PDF4LHC recommendations [33] and are found to range from 0.2% to 2.0%. The uncertainty in the acceptance due to the modelling of multiple-parton interactions is estimated from the difference in acceptance with and without multiple-particle interactions for the gluon-gluon fusion simulation sample and ranges from 2.9% to 25%.

Statistical procedure
A profile-likelihood-ratio test statistic [105] is used to search for a localised excess over a smoothly falling background in the m Zγ distribution of the data, as well as to quantify its significance and estimate either its production cross section or signal strength.
The extended unbinned likelihood function L(α, θ) is given by the product of a Poisson term, constructed from the number of observed events, n, the expected event yield, N, and the probability density function of the invariant mass distribution for each candidate event i, f tot (m i Zγ , α, θ) [22]: where L tot is the total integrated luminosity, ε is the signal efficiency, and σ(pp → H) SM · B(H → Zγ) SM is predicted by the SM. The theoretical uncertainties are taken into account as described in Section 6.
The probability density function for the invariant mass ( f tot (m i Zγ , α, θ)) is built from the probablity density functions f sig and f bkg describing the signal and background invariant mass distributions, respectively: The index c indicates the category. The θ bkg are nuisance parameters that determine the shape of the background. The nuisance parameters associated with the uncertainties in the signal parameterisation, efficiency and acceptance are denoted by θ sig . Nuisance parameters associated with uncertainties in the event yield or the m Zγ resolution are assigned log-normal probability density functions, while nuisance parameters associated with the m Zγ signal peak position are assigned Gaussian probability density functions. The nuisance parameters associated with the spurious signal, θ spur , are assigned Gaussian probability density functions.
The probability that the background can produce a fluctuation greater than or equal to an excess observed in data is quantified by the p-value of the α = 0 hypothesis, p 0 , which can also be expressed in terms of number of Gaussian standard deviations and provides an estimate of the local significance of a possible deviation from the expected background. The global significance, corrected for the effect that a deviation can occur anywhere in the search region, is estimated taking into account the trial factors [106]. The compatibility between the data and increasing non-zero values of α is used to set upper limits at the 95% CL on σ(pp → H/X) · B(H/X → Zγ), B(H → Zγ), and the signal strength for pp → H → Zγ, respectively, using a modified frequentist (CL s ) method [107], by identifying the value α for which the value of CL s is equal to 0.05.
The results are derived using closed-form asymptotic formulae [105] for masses up to 1.6 TeV. Due to the small number of events at large m Zγ in the high-mass resonance search, the results for m X > 1.6 TeV are derived using ensemble tests. The expected cross-section limits obtained from the asymptotic formulae agree to better than 10% with those obtained from the ensemble tests up to m X of 1.7 TeV. At high m X , the observed (expected) central values are underestimated by 35% (23%) in the asymptotic approach.

Results
No evidence of a localised excess is visible near the anticipated Higgs mass m H = 125.09 GeV, as shown in Figure 5 where The observed 95% CL limit on σ(pp → H) · B(H → Zγ) is found to be 6.6 times the SM prediction, corresponding to the limit of 547 fb. Assuming SM Higgs boson production, the upper limit on the branching ratio of the Higgs boson decay to Zγ is found to be 1.0%. The expected 95% CL limit on σ(pp → H) · B(H → Zγ) assuming no (a SM) Higgs boson decay to Zγ is 4.4 (5.2) times the SM prediction.
The invariant mass distributions of events satisfying the high-mass selection are displayed for both categories (ee and µµ) in Figure 6 and compared to the background-only fit performed in the fit range 200 GeV < m Zγ < 2500 GeV. The highest masses measured in the eeγ and µµγ final states are 1.47 TeV and 1.57 TeV, respectively. In the fit range, no significant excess is observed with respect to the background-only hypothesis. The largest deviation is observed around m X = 960 GeV corresponding to a local significance of 2.7 σ. The global significance, evaluated using the search region of [250-2400] GeV in mass, is found to be 0.8 σ. The observed and expected upper limits on σ(pp → X) · B(X → Zγ) as a function of m X are shown in Figure 7. The observed limits vary between 88 fb and 2.8 fb for the mass range from 250 GeV to 2.4 TeV at the 95% CL for a spin-0 resonance produced via gluon-gluon fusion, while the expected limits range from 61 fb to 2.7 fb in the same mass range. For a spin-0 resonance produced via vector-boson fusion, the limits are up to 4% lower due to the slightly larger efficiency for this production process. The results are also interpreted in terms of spin-2 resonances. Observed and expected upper limits at the 95% CL on σ(pp → X) · B(X → Zγ) are derived and shown in Figure 8 for both the gg and qq processes. The observed limits for the gg (qq) process vary between 117 fb (94 fb) and 3.7 fb (2.3 fb) for the mass range from 250 GeV to 2.4 TeV, while the expected limits range between 82 fb (66 fb) and 3.6 fb (2.2 fb) in the same mass range.
The limits on σ(pp → X) · B(X → Zγ) for high-mass resonances are valid for resonances with a natural width that is small compared to the detector resolution. [GeV] TeV as a function of the high-mass spin-0 resonance's mass, assuming production via gluon-gluon fusion and using the narrow width assumption (NWA). For m X > 1.6 TeV results are derived from ensemble tests in addition to the results obtained using closed-form asymptotic formulae. The shaded regions correspond to the ±1 and ±2 standard deviation bands for the expected exclusion limit derived using asymptotic formulae.  Figure 8: The observed (solid line) and expected (dashed line) upper limit derived at the 95% CL on σ(pp → X) · B(X → Zγ) at √ s = 13 TeV as a function of the spin-2 resonance mass produced via (a) gluon-gluon initial states and (b) qq initial states modelled using the Higgs Characterisation Model (HCM), using the narrow width assumption (NWA). For m X > 1.6 TeV results are derived from ensemble tests in addition to the results obtained using closed-form asymptotic formulae. The shaded regions correspond to the ±1 and ±2 standard deviation bands for the expected exclusion limit derived using asymptotic formulae.

Conclusion
Searches for Zγ decays of the SM Higgs boson (H → Zγ) and of a narrow high-mass resonance (X → Zγ) in 36.1 fb −1 of pp collisions at √ s = 13 TeV at the LHC have been performed with the ATLAS experiment. The observed data are consistent with the expected background. No evidence for the H → Zγ and X → Zγ decays is observed and upper limits are set on σ(pp → H) · B(H → Zγ) at m H = 125.09 GeV and also on σ(pp → X) · B(X → Zγ) as a function of m X . For the Higgs mass of 125.09 GeV, the observed 95% CL upper limit on the σ(pp → H) · B(H → Zγ) is 6.6 times the SM prediction. The search for high-mass Zγ resonances was studied using both the spin-0 and spin-2 interpretations. The observed limit varies between 88 fb and 2.8 fb for the mass range from 250 GeV to 2.4 TeV for a spin-0 resonance, where a resonance produced in gluon-gluon fusion is used as a benchmark model. For spin-2 resonances, LO predictions from the Higgs Characterisation Model are used as benchmarks. The limits for spin-2 resonances range between 117 fb (94 fb) and 3.7 fb (2.3 fb) for the the mass range from 250 GeV to 2.4 TeV for a resonance produced via gluon-gluon (quark-antiquark) initial states. The corresponding expected limits for this mass range vary between 61 fb and 2.7 fb for a spin-0 resonance, and between 82 fb (66 fb) and 3.6 fb (2.2 fb) for a spin-2 resonance produced via gluon-gluon (quark-antiquark) initial states.
(Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [108].     [89] OPAL collaboration, K. Ackerstaff et al., Search for anomalous production of dilepton events with missing transverse momentum in e + e − collisions at √ s = 161 GeV and 172 GeV, Eur. Phys. J. C 4 (1998) 47, arXiv: hep-ex/9710010.    The ATLAS Collaboration 1 Department of Physics, University of Adelaide, Adelaide, Australia