Search for Z$\gamma$ resonances using leptonic and hadronic final states in proton-proton collisions at $\sqrt{s}=$ 13 TeV

A search is presented for resonances decaying to a Z boson and a photon. The analysis is based on data from proton-proton collisions at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 35.9 fb$^{-1}$, and collected with the CMS detector at the LHC in 2016. Two decay modes of the Z boson are investigated. In the leptonic channels, the Z boson candidates are reconstructed using electron or muon pairs. In the hadronic channels, they are identified using a large-radius jet, containing either light-quark or b quark decay products of the Z boson, via jet substructure and advanced b quark tagging techniques. The results from these channels are combined and interpreted in terms of upper limits on the product of the production cross section and the branching fraction to Z$\gamma$ for narrow and broad spin-0 resonances with masses between 0.35 and 4.0 TeV, providing thereby the most stringent limits on such resonances.


Introduction
One of the key aspects of the CERN LHC physics program is the search for new resonances predicted in theories beyond the standard model (SM). Given the fairly stringent limits already set on masses of such resonances in fermionic decay channels (e.g., via dilepton or dijet searches), it is particularly interesting to explore bosonic decay channels, which can dominate if the couplings of a new resonance to fermions are suppressed. Examples of such signatures are decays into a pair of massive bosons: VV and VH, where V represents either a W or a Z boson, and H refers to the recently discovered Higgs boson [1-3]. The latest results from these searches at the LHC are described in Refs. [4][5][6][7][8][9][10][11][12][13][14][15] for the VV channels and in Refs. [10,[16][17][18][19][20][21] for the VH channels.
Diboson decays involving photons, i.e., Wγ, Zγ, and γγ channels, are also important, as the search in the γγ channel demonstrated by contributing significantly to the discovery of the Higgs boson by the ATLAS and CMS Collaborations in 2012 [1-3]. While a resonance decaying to diphotons cannot be a vector or an axial vector, due to the Landau-Yang theorem [22,23], having one of the two bosons massive alleviates this constraint. Thus, charged (neutral) bosons of spin 0, 1, or 2 can be sought in the Wγ (Zγ) channel, leading to a broad search program. Such resonances are predicted in many extensions of the SM, such as technicolor [24] and little Higgs [25] models, as well as models with an extended Higgs boson sector [26,27] or with extra spatial dimensions [28,29].
In this paper, we describe a search for Zγ resonances in the leptonic ( , where refers to e or µ) and hadronic decay channels of the Z boson, as well as the combination of these channels. While the results in this paper are interpreted in terms of spin-0 resonances, they are broadly applicable to spin-1 and spin-2 states, as the signal acceptance depends only weakly on the spin of the resonance [30]. Similar searches in a combination of leptonic and hadronic decay channels of the Z boson have been recently published by ATLAS [31] at √ s = 13 TeV and by CMS at √ s = 8 and 13 TeV [32], based on significantly smaller integrated luminosities. Other searches for Zγ resonances have been performed only in the dilepton channels. These include searches by the L3 Collaboration at the CERN LEP [33] and the D0 Collaboration at the Fermilab Tevatron [34,35]. At the LHC, they have been done by ATLAS at √ s = 7, 8, and 13 TeV [30,36,37], and CMS at √ s = 8 and 13 TeV [38], as well as by ATLAS and CMS using the combined 7 and 8 TeV data [39,40], and by ATLAS using the 13 TeV data [30] in the context of a search for the H → Zγ decay.
The present search is for a resonance with a relatively narrow width appearing as an excess over the smooth Zγ invariant mass (m Zγ ) spectrum constructed from an energetic photon and the Z → or Z → qq decay products. While a search in the leptonic channels has lower SM backgrounds, resulting in higher sensitivity to resonance masses <1 TeV, at larger mass values, where backgrounds are small, the hadronic channels, with their higher branching fraction, dominate the sensitivity. The backgrounds in both channels are determined directly from fits to data.
The Z boson decays in the leptonic channels are reconstructed using an electron or a muon pair. The dominant backgrounds in the γ channel are the irreducible contribution from continuum Zγ production and the reducible backgrounds from either final-state radiation in Z → events or from Z boson production in association with one or more jets (Z+jets), where a jet is misidentified as a photon.
The Z boson decay products in the hadronic channels can be reconstructed either as two wellseparated small-radius jets, or as a single large-radius jet (J) resulting from the merging of the two quark jets because of the large Lorentz boost of a Z boson produced in the decay of a heavy resonance. The fraction of events corresponding to the merged topology, which has low background from SM sources, increases with the mass of the resonance. To optimize signal relative to background, in this paper we consider just the merged jet topology (Jγ), and thus the search in the hadronic channels is focused on relatively high resonance masses where the trigger for these events is efficient. We use jet substructure techniques to infer the presence of two subjets, and dedicated b tagging algorithms to identify those subjets that originate from b quark fragmentation. This provides the means to distinguish a signal from the dominant background from prompt-photon and QCD multijet production, with one of the jets spuriously passing jet substructure requirements (and, in the latter case, with another jet mimicking a photon).

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
The silicon tracker measures charged particles within the pseudorapidity range |η| < 2.5. It consists of 1440 silicon pixel and 15 148 silicon strip detector modules. For nonisolated particles of 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90  µm in the transverse (longitudinal) impact parameter [41] The ECAL consists of 75 848 lead tungstate crystals, which provide coverage in pseudorapidity |η| < 1.48 in a barrel region (EB) and 1.48 < |η| < 3.0 in two endcap regions (EE). Preshower detectors consisting of two planes of silicon sensors interleaved with a total of 3X 0 of lead are located in front of each EE detector.
In the region |η| < 1.74, the HCAL cells have widths of 0.087 in pseudorapidity and 0.087 in azimuth (φ). In the η-φ plane, and for |η| < 1.48, the HCAL cells map on to 5 × 5 arrays of ECAL crystals to form calorimeter towers projecting radially outwards from close to the nominal interaction point. For |η| > 1. 74, the coverage of the towers increases progressively to a maximum of 0.174 in ∆η and ∆φ. Within each tower, the energy deposits in ECAL and HCAL cells are summed to define the calorimeter tower energies, subsequently used to provide the energies and directions of hadronic jets.
The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7% for nonshowering electrons in the barrel region to 4.5% for showering electrons in the endcaps [42].
Muons are measured in the pseudorapidity range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. Matching muons to tracks measured in the silicon tracker results in a relative transverse momentum resolution for muons with 20 < p T < 100 GeV of 1.3-2.0% in the barrel and better than 6% in the endcaps. The p T resolution in the barrel is better than 10% for muons with p T up to 1 TeV [43].
Events of interest are selected using a two-tiered trigger system [44]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [45].

Data sets and event selection
The data used in this search correspond to an integrated luminosity of 35.9 fb −1 recorded by the CMS experiment at √ s = 13 TeV in 2016. The high instantaneous luminosity delivered by the LHC results in additional interactions in the same or neighboring bunch crossings (pileup) as the hard scattering interaction. The average number of pileup interactions in the 2016 data set is around 23.
In the eeγ channel, the selected events are required to pass a double-photon trigger with the transverse momentum p T > 60 GeV and pseudorapidity |η| < 2.5 requirements on both photon candidates. Since the photon trigger requirements do not include any track veto, this trigger is equally efficient in selecting photon and electron candidates. A combination of single-muon triggers requiring p T > 50 GeV and |η| < 2.4 on a muon candidate are used in the µµγ channel. In the Jγ channel, we require a logical "OR" of several triggers with the separate requirements: the scalar sum of transverse energies of all reconstructed jets (H T ) is above 800 or 900 GeV; a jet is present with the transverse energy above 500 GeV; and a photon candidate is present with p T > 165 or 175 GeV and |η| < 2.5. We determine the selection efficiency for these trigger combinations using unbiased data samples collected with different triggers. The triggers are found to be 98-100% efficient with respect to the offline selection, for the entire mass range used, in all three channels. The small residual inefficiency is taken into account when calculating the signal acceptance.
Simulated signal events of spin-0 resonances decaying to Zγ are generated at leading order (LO) in perturbative QCD using PYTHIA 8.205 [46] with the CUETP8M1 [47,48] underlyingevent tune. Several samples are generated with masses ranging from 0.3 to 4.0 TeV. Two resonance width assumptions were used in the simulation: one, termed "narrow", has its width (Γ X ) set to 0.014% of the resonance mass (m X ), and the second, referred to as "broad", has Γ X /m X = 5.6%. The first choice corresponds to a resonance with a natural width much smaller than the detector resolution. The second choice facilitates a direct comparison with the previous CMS publications [32,38]. We assume no interference between signal and the SM nonresonant Zγ production.
Simulated background events do not enter the analyses directly, as the backgrounds are obtained from fits to data, but are used to assess the accuracy of the background model and to optimize event selection. Standard model nonresonant Z( )γ production, which is expected to be the dominant background process in the γ channel, is generated at next-to-leading order (NLO) accuracy using the MADGRAPH5 aMC@NLO 2.3.3 generator [49,50]. The Z( )+jets events with a jet misidentified as a photon, which constitute a subdominant source of background, are generated at LO using MADGRAPH5 aMC@NLO, as are the dominant γ+jets and QCD multijet events, as well as subdominant hadronically decaying W+jets and Z+jets backgrounds in the Jγ channel. All background events are processed with PYTHIA for the descrip-tion of fragmentation and hadronization.
All simulated samples were produced using NNPDF3.0 [51] parton distribution functions (PDFs), processed with the full CMS detector model based on GEANT4 [52], and reconstructed with the same suite of programs as used for collision data. Pileup effects are taken into account by superimposing minimum bias events on the hard scattering interaction. The simulated samples are reweighted to match the reconstructed vertex multiplicity distribution observed in data.
The particle-flow (PF) event algorithm [53] aims to reconstruct and identify each individual particle in an event, based on an optimized combination of information from the various elements of the CMS detector. The energy of photons is directly obtained from the ECAL measurement, corrected for zero-suppression effects. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons is obtained from the curvature of the reconstructed muon track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energy deposits.
The events must contain at least one reconstructed primary vertex with at least four associated tracks, with transverse (longitudinal) coordinates required to be within 2 (24) cm of the nominal collision point. The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [54,55] with the tracks assigned to the vertex as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the transverse momenta of those jets.
Electron candidates must pass loose identification criteria based on the shower shape variables, on the ratio of energy deposits in the associated HCAL and ECAL cells, on the geometrical matching between the energy deposits and the associated track, and on the consistency between the energy reconstructed in the calorimeter and the momentum measured in the tracker [42].
Muon candidates are reconstructed from tracks found in the muon system that are associated with the tracks in the inner tracking systems. One muon candidate is required to pass a loose identification [56]. Another muon candidate is required to pass a tighter identification based on the numbers of associated hits found in the pixel and strip trackers, on the numbers of hits and track segments in the muon detector, and on criteria for the matching between the silicon detector track and the muon track segments [56].
Leptons are required to be isolated from other energy deposits in the event. This is expected for signal leptons from Z boson decays, but is not the case for backgrounds from nonprompt leptons originating, e.g., from b hadron decays. The relative isolation of a lepton is defined as the scalar sum of the transverse momenta of all relevant PF candidates within a cone around the lepton, divided by the p T of the lepton candidate. For an electron, the cone size ∆R = (∆η) 2 + (∆φ) 2 depends on its p T : for p T ≤ 50 GeV, 10 GeV p T , for 50 < p T ≤ 200 GeV, and 0.05, for p T > 200 GeV. (1) The electron isolation is based on the photons, and charged and neutral hadrons found in the isolation cone. Charged hadrons originating from pileup vertices are excluded from the sum. The contribution to the isolation sum from neutral pileup particles is accounted for by using the average energy density method [57]. The varying isolation cone radius in Eq.
(1) takes into account the aperture of b hadron decays as a function of their p T , and reduces the inefficiency from accidental overlap of electrons from Z boson decays and jets. For muons, a fixed cone of a size ∆R = 0.3 is used, and the isolation is based on all charged-particle tracks within the isolation cone, excluding the candidate muon track. In the case of two spatially close muons in the event, with overlapping isolation cones, both muons are excluded from each isolation sum. This procedure, together with the use of a variable cone size for electron isolation, ensures high lepton identification efficiency even in the topologies where a Z boson has a high Lorentz boost, as expected for Z bosons produced in a decay of a heavy resonance. The relative isolation of electron and muon candidates is required to be less than 0.1.
Photon identification is based on a multivariate analysis, employing a boosted decision tree algorithm [58,59]. The inputs to the algorithm include shower shape variables, isolation sums computed from PF candidates in a cone of radius ∆R = 0.3 around the photon candidate, and variables that account for the dependences of the shower shape and isolation variables on the pileup [60]. In addition, a conversion-safe electron veto [60] is applied. Photon candidates are required to pass a working point that corresponds to a typical photon reconstruction and identification efficiency of 90%, in the photon p T range used in the analysis.
In the Jγ channel, large-cone jets are used to reconstruct hadronically decaying highly Lorentz boosted Z boson candidates. Jets are reconstructed from PF candidates clustered using the antik T algorithm [54] with a distance parameter of 0.8. Charged hadrons not originating from the primary vertex are not considered in the jet clustering. Corrections based on the jet area [57] are applied to remove the energy contribution of neutral hadrons from pileup interactions. The energies of the jets are further corrected for the response function of the calorimeter. These corrections are extracted from simulation and confirmed with in situ measurements using the energy balance in dijet, multijet, γ+jet, and leptonically decaying Z+jet events [61,62]. Additional quality criteria are applied to jets to remove rare spurious noise patterns in the calorimeters, and also to suppress leptons misidentified as jets. The jet energy resolution amounts typically to 15% at 10 GeV, 8% at 100 GeV, and 4% at 1 TeV. Jets must have p T > 200 GeV and |η| < 2.0. The requirement on the jet η suppresses background from γ+jets and QCD multijet events, and ensures that the core of the jet is within the tracker volume of the detector (|η| < 2.5). The latter requirement is important for subsequent b quark tagging.
Events in the γ channel are required to have two same-flavor leptons (electrons or muons) and a photon. Additionally, leptons in the µµγ channel are required to have opposite electric charge. This requirement is not used in the eeγ channel due to a nonnegligible probability to misreconstruct the charge of an electron candidate because of an energetic bremsstrahlung. The leading electron (muon) is required to have p T > 65 (52) GeV and |η| < 2.5 (2.4). The subleading lepton is required to have p T > 10 GeV and to be found in the same pseudorapidity range as the leading lepton. The photon in the eeγ (µµγ) channel is required to satisfy p T > 65 (40) GeV and |η| < 2.5. Electrons and photons in the ECAL barrel-endcap transition region (1.44 < |η| < 1.57) are excluded from the analysis. In the eeγ channel, the p T thresholds on the electrons and photons in the ECAL endcap region are increased to 70 GeV, in order to ensure a fully efficient trigger. Photons are required to be separated from lepton candidates by ∆R > 0.4, to reduce the background from final-state radiation in Z → events. The invariant mass of the dilepton system is required to be 50 < m < 130 GeV. The minimum dilepton mass requirement suppresses contributions from pp → γγ * events, where an internal conversion of a photon produces a lepton pair. Finally, we require the ratio of the photon p T to m Zγ to be greater than 0.27. This requirement suppresses backgrounds due to jets misidentified as photons, without significant loss in the signal efficiency and without introducing a bias in the m Zγ spectrum. We search for resonances in the m Zγ spectrum above 300 (250) GeV in the electron (muon) channel.
For the Jγ channel, the photon candidates are required to have p T > 200 GeV and to fall within the barrel fiducial region of the ECAL (|η| < 1.44). Events with a photon reconstructed in the endcap region suffer from high γ+jets background and do not add to the sensitivity of the analysis; therefore they are not considered. Photon candidates in the event are required to be separated from large-radius jets by a distance of ∆R > 1.1, which guarantees that the photon isolation cone is not contaminated with the jet constituents.
To identify Z boson candidates in the Jγ channel, the reconstructed large-radius jet mass, evaluated after applying a jet pruning algorithm [63,64], is used. The jet pruning reclusters the jet constituents and eliminates soft, large-angle QCD radiation, which otherwise contributes significantly to the jet mass. The pruning algorithm reclusters each jet starting from its original constituents with the Cambridge-Aachen (CA) algorithm [65] and discards soft and wideangle recombinations in each step of the iterative CA procedure. The pruned jet mass (m pruned J ) is computed from the sum of the four-momenta of the remaining constituents, which are corrected with the same factor as has already been used in the generic jet reconstruction described above. A detailed description of the pruning algorithm can be found in Ref. [66]. For the signal selection, we require a Z candidate to have 75 < m pruned J < 105 GeV.
Finally, in the Jγ channel a requirement is imposed on the ratio of photon p T to the reconstructed Zγ mass of p T /m Zγ > 0.34, with the cutoff chosen, based on a study of simulated signals and backgrounds, to maximize the discovery potential for a narrow resonance. We search for resonances with masses m Zγ > 650 GeV in this channel.
To further discriminate against the QCD multijet and γ+jets backgrounds in the Jγ channel, we categorize the events according to the likelihood of a large-radius jet to contain subjets originating from b quark fragmentation and to contain exactly two subjets. In order to do so, we employ subjet b tagging and N-subjettiness [67] variables (τ N ). The N-subjettiness observable measures the spatial distribution of jet constituents relative to candidate subjet axes in order to quantify how well the jet can be divided into N subjets. Subjet axes are determined by a one-pass optimization procedure, which minimizes N-subjettiness [68]. In particular, the ratio of 2-to 1-subjettiness, τ 21 = τ 2 /τ 1 , offers an excellent separation between the QCD jets and jets from vector boson decays [69], which tend to have lower τ 21 values than the former.
To infer the presence of b quark subjets within a large-radius jet, pruned jets are split into two subjets by reversing the final iteration of the jet clustering algorithm. These subjets are classified according to the probability of their originating from b quarks, based on results from the combined secondary vertex (CSVv2) b tagging algorithm [70,71]. A jet is identified as being consistent with the Z → bb decay when at least one of its subjets satisfies the medium operating point of the CSVv2 algorithm, and the other subjet satisfies the loose operating point. The medium and loose operating points correspond to 70 and 85% (20 and 50%) in the b jet tagging efficiency for p T < 300 GeV (p T = 1 TeV), and 1-2% and 10-15% misidentification probability of a light-flavor jet, respectively. If an event contains a Z → bb candidate, it is classified as "b tagged". For the rest of the events, if the large-radius jet has τ 21 < 0.45, we classify the event as "τ 21 tagged". Otherwise, the event is assigned to the "untagged" category. These three categories are mutually exclusive and are combined for the final result. The additional classification according to the τ 21 value enhances signal sensitivity by 10-15% at low to intermediate signal masses (up to ∼2 TeV), relative to a previous analysis in the hadronic channels [32].

Background and signal modeling 4.1 Background modeling
Simulations in the γ channels indicate that 80-90% of the background remaining after the full event selection is from SM Z boson production accompanied by initial-state photon radiation, with the remainder mostly from Z+jets events, with a jet misreconstructed as a photon. The background m Zγ distributions fall steeply and smoothly with increasing mass, in these channels. Likewise, studies for the Jγ channel based on simulated background samples and on the lower sideband of the jet mass distribution (50 < m pruned J < 70 GeV) in data show that the invariant mass distribution m Zγ of the SM background also falls smoothly in this channel, and that the distributions of kinematic observables derived from the lower jet mass sideband match those for the signal selection.
The background is measured directly in data, through unbinned maximum-likelihood fits to the observed m Zγ distributions, performed separately in each channel. The background in each channel is parametrized with an empirical function. Different families of functions inspired by the ones used in searches for beyond-the-SM phenomena in the dijet, multijet, diphoton, and VV channels at hadron colliders are evaluated in the signal region using simulation ( γ channels) or the lower jet mass sideband (Jγ channel). Examples of these functions are: [72], g(x) = P 0 (1 − x) P 1 /x P 2 +P 3 (ln x) [73], h(x) = P 0 x P 1 exp(P 2 x + P 3 x 2 ) [74], where x = m Zγ / √ s, √ s is the center-of-mass energy (13 TeV), and the number of the fit parameters P i shown is the maximum order considered. The choice of the order within a family of fitting functions for the background distribution is made independently in each channel using the Fisher F-test [75], which balances the quality of the fit against the number of parameters. The choice among the families of functions is optimized based on the results of the bias test described below. The same function [76] was chosen in both the γ and Jγ channels: where P 0 is a normalization parameter, and P 1 , P 2 describe the shape of the invariant mass distribution.
In the γ channels, the absence of significant bias in the fit to background is verified by generating a large number of pseudo-experiments using the simulated background shapes, fitting them with different background models, and measuring the difference between the input and fitted background yields in various m Zγ windows within the entire search range. A pull variable is defined in each window by the difference between the input and fitted yields, divided by the combined statistical uncertainties in the data and the fit. If the absolute value of the median in the pull distribution is found to be >0.5, an additional uncertainty is assigned to the background parametrization. A modified pull distribution is then constructed, increasing the statistical uncertainty by an extra term, denoted as the bias term. The bias term is parametrized as a smooth function of m Zγ and tuned to make the absolute value of the median of the modified pull distribution to be <0.5 in all mass windows. This additional uncertainty is included in the likelihood function by adding to the background model a component with a distribution that is the same as the signal, but a normalization coefficient distributed as a Gaussian of mean zero, and a width equal to the integral of the bias term over the signal mass window, defined as the full width at half maximum. The inclusion of this additional component takes into account the possible mismodeling of the background shape. The bias term in the γ analysis corresponds to ≈0.3 events/GeV at m Zγ = 400 GeV and smoothly falls to ≈2 × 10 −4 events/GeV around m Zγ = 2 TeV. The observed m Zγ invariant mass spectra in data are shown in Fig. 1 for the eeγ (left) and µµγ (right) channels. The results of the fits and their uncertainties at 68% confidence level (CL) are shown by the lines and the shaded bands, respectively.  Figure 1: Observed m Zγ invariant mass spectra in the eeγ (left) and µµγ (right) channels. The best fits to the background-only hypotheses are represented by the red lines, with their 68% CL uncertainty bands given by the gray shadings. Several narrow and broad signal benchmarks with arbitrary normalization are shown on top of the background prediction with the dashed lines. The lower panels show the difference between the data and the fits, divided by the uncertainty, which includes the statistical uncertainties in the data and the fit. For bins with a small number of entries, the error bars correspond to the Garwood confidence intervals [77].
Similarly, in the Jγ channel, the mass spectra in the three analysis categories, derived either from the low jet mass sideband in data or from simulated background samples, are fitted with a variety of alternative functions to generate pseudo-data sets. Additionally, in a set of pseudoexperiments, signals with different mass values and cross sections close to the expected 95% CL limits are injected. The full spectra are fitted with the chosen function of Eq. (2) together with a signal model (discussed in Section 4.2), and the signal cross section is extracted. Distributions of the difference between the data and the fits, divided by the overall uncertainty for the obtained signal cross section, are constructed, and their shapes are found to be consistent with a normal distribution with a mean less than 0.5 and a width consistent with unity.
Thus, any possible systematic bias from the choice of the functional form in the region of low background is proven to be small compared to the statistical uncertainty from the accuracy of the measurements. In the region of large background, the uncertainty in the signal efficiency (discussed in Section 5) completely dominates the effect of the background uncertainty on the limits. Thus we assign the statistical uncertainty in the fits as the only uncertainty in the background predictions. The observed m Zγ invariant mass spectra in the signal region (75 < m pruned J < 105 GeV) for the b-tagged, τ 21 -tagged, and untagged categories, along with the corresponding fits, are shown in Fig. 2.

Signal modeling
The signal distribution in m Zγ is obtained from the generated events that pass the full selection. The signal shape is parametrized with a Gaussian core and two power-law tails, namely an extended form of the Crystal Ball (CB) function [78]. We find this functional form to provide an adequate description for both narrow and broad signals in the entire mass range used in the analysis. To derive the signal shapes for the intermediate mass values where simulation points are not available, a linear morphing [79] of the shapes obtained from the simulation is used. The typical mass resolution for narrow-width signal events is 1% for the eeγ channel, 1-2% for the µµγ channel, depending on the mass of the resonance, and 3% in the Jγ channel.
The product of signal acceptance and efficiency for a narrow resonance in the eeγ (µµγ) channel rises from about 27 (42)% at m Zγ = 0.35 TeV to about 46 (55)% at m Zγ = 2 TeV, and remains steady until 4 TeV. In the Jγ channel, the product of signal acceptance and efficiency for narrow resonances increases from 7 (3)% at 0.65 TeV to 9 (9)% at 4 TeV in the untagged (τ 21 -tagged) category, and is between 2 and 3% for the b-tagged category for the resonance masses between 0.65 and 4 TeV.
For a broad resonance the product of signal acceptance and efficiency is similar to the narrowresonance case up to 2 TeV. At large resonance masses (>2 TeV), however, the effect of rapidly falling PDFs introduces a lower tail in the signal mass distribution. The exact characteristics of this tail are quite sensitive to the resonant line shape. We therefore truncate the mass distribution of the resonance to correspond to the core of the line shape, defined as a window centered on the maximum of the CB function with a width given by ±5 times the CB function parameter σ, describing the standard deviation of its Gaussian core. The tails outside this window are conservatively discarded in the signal acceptance calculations and when fitting the data. As a result, the product of signal acceptance and efficiency falls to 2% at 4 TeV for the eeγ and µµγ events, to 0.2% in the untagged and τ 21 -tagged categories, and to <0.1% in the b-tagged category.

Systematic uncertainties
The statistical uncertainty in the fits of the background function to data is taken for the background uncertainty in all channels.
The following systematic uncertainties in the signal are defined below, and summarized in Table 1: • Integrated luminosity: the uncertainty in the CMS integrated luminosity is based on cluster counting in the silicon pixel detector and amounts to 2.5% [80].
• PDFs: a 1.0-3.5% uncertainty in the signal efficiency that takes into account the variation in the kinematic acceptance of the analysis is estimated using replicas of the NNPDF3.0 set, following the PDF4LHC prescription [81]. The uncertainty in the signal cross section due to the PDF choice is not considered.
• Pileup: the uncertainty due to the pileup description in the signal simulation, evaluated by changing the total inelastic cross section governing the average multiplicity of pileup interactions by ±5.0% [82], translates to a 1.0% uncertainty in the signal acceptance in all channels.
• Trigger: the uncertainty due to the trigger efficiency differences in data and simulation in the γ analysis is estimated with dedicated studies with leptons from Z boson decays and amounts to 1.0 (3.5)% for the eeγ (µµγ) channel. In the Jγ channel, a 2.0% uncertainty covers the variation of the trigger efficiency across the mass range probed in the analysis.
• Photon efficiency: the systematic uncertainty due to the differences in the photon identification efficiency between data and simulation is evaluated with Z → ee events in which the electrons are used as proxies for photons, and amounts to 1.5% [60].
• Lepton efficiency: the systematic uncertainty due to the differences in the lepton identification efficiency in data and simulation is evaluated with Z → ee (µµ) events and amounts to 2.5 (2.0)% in the eeγ (µµγ) channel.
• b tagging efficiency: the uncertainty due to the difference in the b tagging efficiency in data and simulation is estimated from control samples in data and simulation enriched in b quarks [71], and translates into a 15-32% uncertainty in the signal yield in the Jγ channel. It is anticorrelated between the b-tagged and the other two categories, as it induces signal migration between the categories.
• τ 21 tagging efficiency: to account for the difference between the τ 21 distributions in data and simulation, a scale factor of 0.97 ± 0.06 [69] is introduced for simulated signal samples. This translates into an uncertainty of 10-12% in the signal yield in the τ 21 -tagged category and is anticorrelated with that in the untagged category.
• Electron and photon energy scale and resolution: the electron and photon energy scale is known with 0.1-5.0% precision, depending on the energy. This uncertainty is based on the accuracy of the energy scale at the Z boson peak and its extrapolation to higher masses, and translates into a 0.2-4.6 (0.1-2.3)% correlated uncertainty in the m Zγ scale in the eeγ channel (µµγ and Jγ channels). The uncertainty in the electron and photon energy resolution based on the Gaussian smearing evaluated at the Z boson peak translates to a 10 (5)% uncertainty in the m Zγ resolution in the eeγ channel (µµγ and Jγ channels).
• Muon momentum scale and resolution: the muon momentum scale is measured with 0.1-5.0% precision up to p T = 200 GeV, with an additional 0.1-6.0% uncertainty at higher values, resulting in a 0.1-4.6% uncertainty in the m Zγ mass scale in the µµγ channel. A 10% uncertainty in the m Zγ resolution in the µµγ channel is conservatively assigned to account for the uncertainty in the muon momentum resolution.
• Jet energy scale (JES), jet mass scale (JMS), jet energy resolution (JER), and jet mass resolution (JMR): the uncertainties [61,62,66] are propagated to all the relevant quantities, and affect both the signal yield and its shape. The overall effect of these uncertainties added in quadrature corresponds to approximately 5.0% uncertainty in the signal yield, as determined by changing the four-momenta of the jets accordingly and carrying out the full analysis with the modified quantities.

Results
The data are consistent with the background-only expectations in all channels. We set upper limits on the production cross section of heavy spin-0 resonances using the asymptotic approximation [83] of the modified frequentist CL s method [84][85][86], with a likelihood ratio used as a test statistic, and uncertainties incorporated as nuisance parameters with log-normal (normalization) or Gaussian (shape) priors. The limits are set in the mass range between 0.35 (0.30) and 4.0 TeV in the eeγ (µµγ) channel and 0.65-4.0 TeV in the Jγ channel. We note that the asymptotic approximation tends to produce lower cross section limits than the exact CL s calculations in the regions with low background. We tested that the difference is at most 4 (30%) for resonance masses below 1 TeV (around 3 TeV).  approximately 3.0 (2.1) standard deviations for a narrow resonance. The limits on σ(X → Zγ), obtained by combining the two leptonic search channels and taking into account the leptonic branching fraction of the Z boson decays [87], are shown in Fig. 4. The rapid increase in the limit for a broad resonance with a mass above approximately 3 TeV is due to a significant lowmass tail in the resonance line-shape extending outside the truncation window, as discussed in Section 4.2.

The Jγ channel
The observed and expected 95% CL upper limits on the product of signal cross section and branching fraction in the Zγ channel, σ(X → Zγ) for narrow and broad resonances in the b-tagged, τ 21 -tagged, and untagged categories are presented in Fig. 5. The results based on the combination of the three categories for both narrow and broad resonances are shown in Fig. 6. The combination includes the (anti)correlation of systematic uncertainties between the three categories. We observe a small deviation at a mass ≈2 TeV with local significance of 2.7 (3.6) standard deviations for the narrow (broad) resonance width hypothesis. The global significance of this excess is 1.8 (2.8) standard deviations.

Combination of the γ and Jγ channels
The results based on the combination of the γ and Jγ channels are shown in Fig. 7, assuming uncorrelated uncertainties between the leptonic and hadronic channels, except for the uncertainties in the integrated luminosity, PDFs, and photon energy scale, which are taken as fully correlated among all the channels. These are the most stringent limits on resonances decaying in the Zγ channel to date in the mass range probed.

Summary
A search is presented for resonances decaying to a Z boson and a photon. The analysis is based on data from proton-proton collisions at a center-of-mass energy of 13 TeV, corresponding to  decay modes of the Z boson are investigated. In the leptonic channels, the Z boson candidates are reconstructed using electron or muon pairs. In the hadronic channels, they are identified using a large-radius jet, containing either light-quark or b quark decay products of the Z boson, via jet substructure and advanced b tagging techniques. The results from these channels are combined and interpreted in terms of upper limits on the product of the production cross section and the branching fraction to Zγ for narrow (broad) spin-0 resonances with masses between 0.35 and 4.0 TeV, ranging from 50 (100) to 0.3 (1.5) fb. These are the most stringent limits on such resonances to date.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus        [41] CMS Collaboration, "Description and performance of track and primary-vertex reconstruction with the CMS tracker", JINST 9 (2014) P10009, doi:10.1088/1748-0221/9/10/P10009, arXiv:1405.6569.