A search for the $Z\gamma$ decay mode of the Higgs boson in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector

A search for the $Z\gamma$ decay of the Higgs boson, with $Z$ boson decays into pairs of electrons or muons is presented. The analysis uses proton$-$proton collision data at $\sqrt{s}$ = 13 TeV corresponding to an integrated luminosity of 139 fb$^{-1}$ recorded by the ATLAS detector at the Large Hadron Collider. The observed data are consistent with the expected background with a $p$-value of 1.3%. An upper limit at 95% confidence level on the production cross-section times the branching ratio for $pp\to H\to Z\gamma$ is set at 3.6 times the Standard Model prediction while 2.6 times is expected in the presence of the Standard Model Higgs boson. The best-fit value for the signal yield normalised to the Standard Model prediction is $2.0^{+1.0}_{-0.9}$ where the statistical component of the uncertainty is dominant.


Introduction
A new boson [1,2] was discovered in 2012 by the ATLAS and CMS Collaborations. The observed properties of the particle, such as its couplings to Standard Model (SM) elementary particles, its spin and its parity, are so far consistent with the predictions for a SM Higgs boson (H) [3][4][5][6][7]. The mass of this boson was determined by the ATLAS and CMS Collaborations to be m H = 125.09 ± 0.21(stat) ± 0.11(syst) GeV using the LHC Run 1 data set [8]. Subsequent measurements, by both collaborations during the LHC Run 2, are published [9,10] and are consistent with this value.
The SM Higgs boson can decay into Zγ through loop diagrams and the branching ratio is predicted to be B(H → Zγ) = (1.54 ± 0.09) × 10 −3 at m H = 125.09 GeV [11]. It can differ from the SM value for several scenarios beyond the SM, for example, if the Higgs boson were a neutral scalar of different origin [12,13], or a composite state [14]. Different branching ratios are also expected for models with additional colourless charged scalars, leptons or vector bosons that couple to the Higgs boson, due to their contributions via loop corrections [15][16][17].
Final states where the Z boson decays into electron or muon pairs can be efficiently triggered and clearly distinguished from background events produced in pp collisions. In addition, the Z(→ )γ ( = e or µ) final state can be reconstructed completely with good invariant mass resolution and relatively small backgrounds.
Previous searches for the H → Z(→ )γ decay by the ATLAS and CMS Collaborations use the full pp data sets collected at √ s = 7 and 8 TeV [18,19] and partial data sets collected at 13 TeV [20,21]. In all cases, no significant excess of events above the expected background is observed around the Higgs boson mass. Prior to the present study, the ATLAS Collaboration reported an observed (expected for a SM Higgs boson) upper limit on the production cross-section times branching ratio for pp → H → Zγ of 6.6 (5.2) times the SM prediction at 95% confidence level for a Higgs boson with m H = 125.09 GeV using a sample corresponding to an integrated luminosity of 36.1 fb −1 [20]. The CMS Collaboration reported an observed (expected for a SM Higgs boson) upper limit of 7.4 (6.0) times the SM prediction for a Higgs boson with m H = 125 GeV using a sample corresponding to an integrated luminosity of 35.9 fb −1 [21].
An updated search for decays of the Higgs boson into Zγ with the Z boson decaying into electron or muon pairs is detailed in this letter. The search uses pp collision data recorded at √ s = 13 TeV with the ATLAS detector at the LHC from 2015 to 2018, corresponding to a total integrated luminosity of 139 fb −1 . There are important improvements in this analysis compared with the previous one [20], including an increase in the size of the data set, an improved event categorisation, and optimised lepton and photon identification criteria. Events with at least one photon and two electrons or muons of opposite charge are classified into six mutually exclusive categories which are designed to maximise the sensitivity to the presence of the SM Higgs boson decaying into Zγ. The dominant background is the irreducible non-resonant production of Z bosons together with photons. A simultaneous fit to the reconstructed Zγ invariant mass distributions in all the categories is performed to extract the overall H → Zγ signal yield.

ATLAS detector and data sample
The ATLAS detector [22,23] 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 (ID) surrounded by a thin superconducting solenoid, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer (MS) incorporating three large air-core toroidal magnets with eight coils each.
A two-level trigger system [24] 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 high-level trigger which runs algorithms similar to those in the offline reconstruction software, reducing the event rate to approximately 1 kHz from the maximum L1 rate of 100 kHz.
The events were collected with triggers requiring either one or two electrons or muons in the event. Due to the increasing luminosity, the transverse momentum (p T ) thresholds were increased slightly during the data-taking periods. At the highest instantaneous luminosity recorded, the lowest p T threshold for the single-muon trigger was 26 GeV. For the dimuon trigger, asymmetric p T thresholds of 22 GeV and 8 GeV were used. The p T threshold was 26 GeV for the single-electron trigger and 17 GeV for both electrons in the dielectron trigger. These lowest-threshold triggers were complemented by triggers with higher thresholds but looser lepton identification criteria. For H → Zγ events that pass the full analysis selection, described in Section 4, the trigger selection is 95.6% efficient for the eeγ final state and 92.2% efficient for the µµγ final state. After applying trigger and data quality requirements, the integrated luminosity of the data used in this search corresponds to 139 ± 2.4 fb −1 . The average number of pp interactions per bunch crossing (pile-up) ranged from about 13 in 2015 to about 39 in 2018, with a peak instantaneous luminosity of 2 × 10 34 cm −2 s −1 . The production of the SM Higgs boson via gluon-gluon fusion (ggF), via vector-boson fusion (VBF), in association with a vector boson (VH, where V is a W or a Z boson) and with a top-quark pair (ttH) was modelled with the P -B v2 Monte Carlo event generator [61][62][63][64], as described in Ref. [65] and summarised in Table 1. P 8 [27] was used to simulate the H → Zγ decay as well as to provide parton showering, hadronisation and the underlying event. Other Higgs boson production processes are not considered as their contributions to the total Higgs boson production cross-section are of the order of 0.1% or less. All production modes of the Higgs boson contribute to the signal in this analysis and their relative yields were fixed to the SM predictions. Contributions from H → µµ (where the reconstructed photon originates from QED final-state radiation) were evaluated using samples produced in the same manner and are considered as a potential background in this analysis. The impact of interference between Higgs bosons decays with the same final-state signature (H → γ * γ and H → µµ) is expected to be negligible [66] and is neglected.

Simulation samples
Additional samples of Higgs bosons produced by gluon-gluon fusion are used for studies of theoretical uncertainties. A sample with multiple parton interactions disabled is used to study the uncertainties in the signal acceptance related to the modelling of non-perturbative quantum chromodynamics (QCD) effects. A sample generated with M G 5_aMC@NLO [67] using the NNPDF30 PDF [68] set, which includes up to two jets at next-to-leading-order (NLO) accuracy in QCD using the FxFx merging scheme [67,69], is used to study the ggF acceptance in the analysis categories.
The background in this analysis originates mainly from non-resonant production of a Z boson and a photon (Zγ), with a smaller contribution from the production of Z bosons in association with jets (Z+jets), with one jet misidentified as a photon.
A large sample of background Zγ events was simulated with the S v2.2.2 [70] generator using a fast simulation of the calorimeter response [71]. It was produced at NLO precision in QCD for up to one additional parton and leading-order (LO) accuracy in QCD for up to three additional partons using the NNPDF3.0nnlo PDF set [72]. The matrix elements were matched and merged with the S parton shower [73,74] using the MEPS@LO prescription [75][76][77][78].
The electroweak production of Zγ j j with jets originating from the fragmentation of partons arising from electroweak vertices is also considered. It was generated at LO accuracy in QCD using M -G 5_aMC@NLO 2.3.3 with no additional partons in the final state, which is orthogonal to the Zγ simulation with Sherpa. The NNPDF30 LO PDF set [72] was used for the generation of the events, and the hadronisation, parton shower and the underlying event of the events was modelled using P 8.212 with the A14 parameter tune [79].
The background originating from Z+jets is estimated using a data-driven technique which is validated using a sample simulated with the P -B v1 MC generator [62][63][64]80] at NLO accuracy in QCD. It was interfaced to P 8.186 for the modelling of the parton shower, hadronisation and the underlying event with parameters set according to the AZNLO tune [81]. The CT10 PDF set [82] was used for the hard-scattering processes, whereas the CTEQ6L1 PDF set [83] was used for the parton shower. Table 1: Higgs boson MC samples produced with P -B v2 with the techniques used and their precision in α s for the event generation (gen.). The precision of the total cross-section used in the sample normalisation is also reported. The cross-section of the qq →ZH sample is normalised to take into account contributions from both qq →ZH and gg →ZH.

Event selection, reconstruction and categorisation
Events are required to have two same-flavour opposite-charge leptons, ( = e, µ), and at least one photon candidate arising from a primary vertex candidate, reconstructed from charged particles (tracks) in the ID with transverse momentum p T > 500 MeV. This primary vertex candidate is required to be the one with the largest sum of the squared transverse momenta of the associated tracks.
Muon candidates are required to have a high-quality track in the ID or the MS, satisfy the medium identification criteria, be within |η| < 2.7 and have p T > 10 GeV [31]. Electron and photon candidates are reconstructed from topological clusters of energy deposits in the EM calorimeter, and in the case of an electron, a track in the ID matched to the cluster [30,91]. Electron and photon candidates in the transition region between the barrel and endcap EM calorimeters, 1.37 < |η| < 1.52, are excluded.
Electrons are required to satisfy loose likelihood-based identification criteria [30], have p T > 10 GeV and be within |η| < 2.47, while photon candidates are required to satisfy p T > 10 GeV, |η| < 2.37 and the tight identification criteria. The lepton and photon candidates are also required to be isolated from additional activity in the tracking detector and in the calorimeters. Contributions to the energy deposited in the calorimeters originating from the underlying events and pile-up are corrected for on an event-by-event basis using the method described in Refs. [92][93][94].
In order to ensure that muon and electron candidates originate from the primary vertex, it is required that the longitudinal impact parameter, ∆z 0 , computed relative to the primary vertex position, satisfies |∆z 0 · sin θ| < 0.5 mm, where θ is the polar angle of the track. Additionally, to suppress leptons from heavy-flavour decays the significance of the transverse impact parameter d 0 calculated relative to the measured beam-line position must satisfy Jets are reconstructed from topological clusters [95] using the anti-k t algorithm [96] with a radius parameter of 0.4. They are required to have p T > 25 GeV and |η| < 4.4. Jets produced in pile-up interactions are suppressed by requiring that those with p T < 60 GeV and |η| < 2.4 pass a selection based on a jet vertex tagging algorithm [97], which is 92% efficient for jets originating from the hard interaction.
An overlap removal procedure is applied to the selected lepton, photon and jet candidates. If two electrons share the same track, or the separation between two electron energy clusters satisfies |∆η| < 0.075 and |∆φ| < 0.125, then only the highest-p T electron is retained. Electron candidates that fall within ∆R = 0.02 of a selected muon candidate are also discarded. In order to suppress the events arising from QED final-state radiation (FSR), photon candidates within a ∆R = 0.3 cone around the leptons of the Z boson candidate are rejected. Jet-lepton and jet-photon overlap removal are also performed by removing the jet if its axis is within a cone of size ∆R = 0.2 around one of the leptons or the photon.
The Z boson candidates are reconstructed from two same-flavour opposite-charge leptons satisfying the selection criteria. The leptons of the Z boson candidate are additionally required to be consistent with being accepted by at least one of the triggers that the event passed. It is required that the p T of the leptons that are associated with the single-lepton or dilepton triggers is 1 GeV above the trigger threshold. The invariant mass of the Z boson candidates is required to be between 50 GeV and 101 GeV before applying the mass resolution improvements described in the following.
The resolution of the invariant mass of the Z → µµ candidates is improved by 3% by correcting the muon momenta for collinear FSR (∆R < 0.15), using identified photons in the EM calorimeter [65]. A constrained kinematic fit is applied to the dilepton invariant mass [98] for all Z boson candidates. This fit makes use of the measured values of the Z boson mass and width, and leads to a 14% improvement in the mass resolution of the Higgs boson candidates.
After the FSR correction and applying the kinematic fit, Z boson candidates are required to have an invariant mass within 10 GeV of the Z boson mass, m Z = 91.2 GeV [99]. If an event has multiple Z boson candidates which pass all requirements, the candidate with the mass closest to the Z boson mass is chosen. Fewer than 1% of the simulated signal events that pass the final H → Zγ selection have more than two leptons (dominated by events produced via VH and ttH) and less than 0.1% of the signal events have a Z boson candidate that does not match a true Z boson.
The Higgs boson candidate is reconstructed from the Z boson candidate and the highest-p T photon candidate in the event. The invariant mass of the final-state particles, m Zγ after correction for FSR and by the kinematic fit, is required to satisfy 105 GeV < m Zγ < 160 GeV. Finally, to reduce background contamination and simplify the background modelling, an additional requirement is placed on the transverse momentum of the photon such that p γ T /m Zγ > 0.12. For SM H → Z(→ )γ events, the reconstruction and selection efficiency (including kinematic acceptance) is 20.4% varying by a maximum of 2% depending on the production mode.
In order to improve the sensitivity to a H → Zγ signal, the selected events are classified into six mutually exclusive categories with different expected signal-to-background ratios and mass resolutions. Categories are defined according to the lepton flavour and event kinematics. Additionally, a boosted decision tree (BDT) trained to separate VBF signal events from other Higgs boson production modes and backgrounds is used to define a category of events with at least two jets. If there are more than two jets in an event, the two highest-p T jets are used. The kinematic variables used in the categorisation and as input to the BDT are: • The p T of the highest-p T jet, p j1 T . • The pseudorapidity difference between the two jets, ∆η j j .
• The minimum ∆R between the Z boson or photon candidate and either of the two jets, ∆R min γ or Z, j .
• The invariant mass of the two jets, m j j .
• The azimuthal separation between the Zγ system and the system formed by the two jets, ∆Φ Zγ, j j .
• The azimuthal angle between the dilepton system and the photon, ∆Φ Z,γ .
• The component, p Tt , 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 . This quantity is strongly correlated with the transverse momentum of the Zγ system, but has better experimental resolution [101,102].
Events with two or more jets with ∆η j j > 2 that have a BDT output score larger than 0.87 are classified into a VBF-enriched category. Of the remaining events, those that satisfy the requirement on p γ T /m Zγ > 0.4 are classified into a High relative p T category while the others are separated into four categories depending on the lepton flavour and a cut at p Tt = 40 GeV. The boundaries of the categories are selected to maximise the expected signal significance.
VBF events are estimated to constitute 72% of the signal in the VBF-enriched category. The High relative p T and High p Tt categories are expected to be enriched in VBF, VH and ttH events as these production modes have on average higher Higgs boson p T than ggF production. As the continuum Zγ background has on average lower Higgs boson candidate p T than the SM Higgs boson, the signal-to-background ratio is expected to be higher in these categories than in the other categories, as shown in Table 2. The table also summarises the observed number of events in data in the Zγ mass range of 105-160 GeV and the expected number of signal (S 68 ) and background (B 68 ) events in a m Zγ window containing 68% of the expected signal. In addition the width of the window containing 68% of the SM signal (w 68 ), which quantifies the mass resolution and is defined as the difference between the 84th and the 16th percentile of the signal mass distribution, is reported. B 68 is estimated from the simulated background samples normalised to the data yields in each category excluding the 120-130 GeV mass region. The categorisation improves the expected sensitivity, based on simple number counting, by approximately 50%.

Signal and background modelling
The signal and background yields are extracted from a fit to the m Zγ distribution observed in data by assuming parametric models for both the signal and backgrounds. For the signal, the expected acceptance and parameters that describe the shape are obtained from simulated signal samples. For the background, the models are chosen using simulated background samples and the values of their parameters are determined by a fit to data mass spectra. ), a class of functions given by

and a sum of power functions (Σ
Zγ ) where f i and p i are the free parameters of the models. The choice of analytical model of the background and the m Zγ range used for the final fit is optimised in each category using the background template mentioned earlier. The optimisation procedure includes a limit on the amount of bias in the extracted signal yield (also referred to as spurious signal) and a requirement on the fit quality, and it prefers models with fewer free parameters to maximise the statistical sensitivity to the signal [1]. For each category used in the analysis, the bias due to the spurious signal is estimated by performing a signal-plus-background fit to the m Zγ background-only distribution with m H varied between 123 GeV and 127 GeV. The maximum number of signal events derived from these fits in each category constitutes the spurious-signal systematic uncertainty. A requirement that the spurious signal be less than 50% of the expected statistical uncertainty in the signal yield is applied when selecting the background modelling function. Stricter requirements on the spurious signal are not possible due to the statistical uncertainty in the m Zγ background template. In addition, the χ 2 probability of the background-only fit is required to be larger than 1%. The fit range is optimised in each analysis category by varying the lower and upper bounds in 5 GeV steps within the ranges 105-115 GeV and 140-160 GeV, respectively. The optimal fit range and function are selected to achieve the highest signal significance while fitting the expected mass distribution of background plus SM signal. The significance evaluation also includes the spurious-signal systematic uncertainty. The selected background functional form and fit range in each category are detailed in Table 3.

Systematic uncertainties
The dominant experimental uncertainty in the signal yield is the spurious signal from the choice of background model. The spurious signal corresponds to as much as 50% of the statistical error in the expected signal yield per category, due to a limited number of simulated background events. It introduces a 28% systematic uncertainty in the signal strength, defined as the ratio of the observed to expected signal yield; however, it only increases the total uncertainty in the expected signal strength by 5.6% because of the large statistical uncertainty. The contribution from H → µµ is estimated using simulated events and amounts to about 1.7% inclusively, and up to 3.3% in individual categories. An uncertainty in this contribution taken from the latest ATLAS measurement [106] results in an uncertainty of 2.1% in the expected signal yield. The combined uncertainty in the signal yield in any category due to the reconstruction, identification, isolation, and trigger efficiency measurements [30,31] is no more than 2.6%, 2.4% and 1.6% for photons, electrons and muons, respectively. Pile-up also affects the lepton and photon identification efficiency but contributes a negligible amount to the uncertainty (0.2%). The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [107, 108].
The theoretical uncertainties in the predicted signal yield originate from uncertainties in the predicted branching ratio (5.7%) [11] as well as from uncertainties in the modelling of the production cross-section and kinematics of the Higgs boson due to missing higher-order QCD calculations (5.3%), that are dominated by uncertainties in the QCD renormalisation and factorisation scales (5.2%). Smaller effects originate from the parton shower modelling uncertainty (1.3%), PDFs (2.5%) and α s (1.9%). The uncertainties in the Higgs boson event kinematics due to missing higher-order QCD calculations impact the distribution of signal events amongst the analysis categories. They are evaluated using an extension of the Stewart-Tackmann method [11,109], based on inputs from Refs. [110][111][112]. Details of how the uncertainty in the acceptance of ggF events in the VBF-enriched category and all other categories is evaluated can be found in Refs. [113] and [20], respectively. Additionally, to account for the uncertainties in the modelling of jet kinematics in ggF events the category acceptance is compared with the acceptance derived from the M G 5_aMC@NLO sample. The Higgs boson and jet kinematics particularly affect the ggF signal acceptance in the VBF-enriched category (37%) and High relative p T category (21%). The effect of parton shower modelling, PDFs and α s on the distribution of signal events amongst the analysis categories is less than 11%, 1% and 2%, respectively. The expected signal yield in the VBF category is also affected (14%) by the jet energy scale, jet energy resolution and vertex tagging efficiency [114].
The uncertainty in the modelling of the signal shape varies between analysis channels. The uncertainty in the mass resolution (σ CB ) is dominated by the uncertainty in the electron and photon energy resolution (< 3.4%) and in the muon momentum resolution (< 3.6%). The uncertainty in the signal position (µ CB ) from the uncertainty in electron, photon and muon calibration (< 0.15%) is less than the uncertainty in the assumed Higgs boson mass of 0.19% [8]. The impact of the signal model uncertainty on the signal strength is less than 2%.

Results
A profile-likelihood-ratio test statistic [115] is used to search for a localised excess of events above the expected background by performing a fit to the m Zγ spectra in the various categories. In the same manner as was done in previous searches for H → Zγ [20], the likelihood is built from the product of Poisson probability terms across all categories with two contributions: non-resonant background, and Higgs boson signal. The likelihood includes terms for the systematic uncertainties discussed in Section 6 implemented as nuisance parameters. The nuisance parameters describe the systematic uncertainties, which are parameterised as Gaussian or log-normal priors. Upper limits are set on the Higgs boson production cross-section at 95% confidence level (CL) using the modified frequentist formalism [116].
The invariant mass distributions of the Zγ events for the various categories are shown in Figure 1 with the background-only fit superimposed. The expected Higgs boson signal normalised to 20 times the SM prediction for m H =125 GeV is also shown. At m H =125.09 GeV, the observed (expected with a SM Higgs boson) p-value under the background hypothesis is 0.013 (0.123), which corresponds to a significance of 2.2 σ (1.2 σ). A weighted sum of all categories with the fitted signal-plus-background model superimposed is shown in Figure 2. The events are weighted by ln(1 + S 68 /B 68 ), where S 68 and B 68 are defined in Section 4.
The best-fit value for the H → Zγ signal strength, defined as the ratio of the observed to the predicted SM signal yield, is found to be 2.0 ± 0.9(stat.) +0.4 −0.3 (syst.) = 2.0 +1.0 −0.9 (tot.) with an expected value of 1.0 ± 0.8(stat.) ± 0.3(syst.) assuming the presence of the SM Higgs boson. The total uncertainty is dominated by the statistical component from the data. The systematic component of the total uncertainty is dominated by the spurious-signal uncertainties.
The observed 95% CL limit on the H → Zγ signal strength is found to be 3.6 times the SM prediction compared with an expected value of 1.7 (2.6) assuming no (SM) Higgs boson decays into Zγ. The observed upper limit on σ(pp → H) · B(H → Zγ) is 305 fb at 95% CL. Assuming the SM Higgs boson production cross-section, the upper limit at 95% CL on B(H → Zγ) is found to be 0.55%.      [11] D.