Search for dark photons in decays of Higgs bosons produced in association with Z bosons in proton-proton collisions at √ s = 13 TeV The CMS collaboration

A search is presented for a Higgs boson that is produced in association with a Z boson and that decays to an undetected particle together with an isolated photon. The search is performed by the CMS Collaboration at the Large Hadron Collider using a data set corresponding to an integrated luminosity of 137 fb−1 recorded at a center-of-mass energy of 13 TeV. No significant excess of events above the expectation from the standard model background is found. The results are interpreted in the context of a theoretical model in which the undetected particle is a massless dark photon. An upper limit is set on the product of the cross section for associated Higgs and Z boson production and the branching fraction for such a Higgs boson decay, as a function of the Higgs boson mass. For a mass of 125 GeV, assuming the standard model production cross section, this corresponds to an observed (expected) upper limit on this branching fraction of 4.6 (3.6)% at 95% confidence level. These are the first limits on Higgs boson decays to final states that include an undetected massless dark photon.


Introduction
Following the discovery of a Higgs boson by the ATLAS and CMS Collaborations [1-3], a primary focus of the LHC physics program has been the study of the properties of the new particle. The observation of a sizable branching fraction of the Higgs boson to invisible or almost invisible final states [4][5][6][7] would be a strong sign of physics beyond the standard model (BSM). Studies of the new boson at a mass of about 125 GeV [8,9] show no significant deviation from the standard model (SM) Higgs boson hypothesis, and measurements of its couplings constrain its partial decay width to undetected decay modes [10,11]. Assuming that the couplings of the Higgs boson to W and Z bosons are smaller than the SM values, an upper limit of 38% has been obtained at 95% confidence level (CL) on the branching fraction of the 125 GeV Higgs boson to BSM particles [11]. This paper presents a search for a scalar boson H produced in association with a Z boson and decaying to an undetected particle together with a photon. Several BSM models predict Higgs boson decays to undetected particles and photons [7,12,13]. In this search, the target final state is Z(→ )H(→ γγ D ), where = e, µ, and γ D is a massless dark photon that couples to the Higgs boson through a charged dark sector [14][15][16][17], and is undetected in the CMS detector. The branching fraction to such an invisible particle and a photon, B(H → invisible + γ), can be as large as 5%, and be consistent with all model parameters and current LHC constraints [15]. A Feynman diagram for such a process is shown in figure 1. While the main focus is the case where the production cross section (σ ZH ) is assumed to be the same as that for the SM-like Higgs boson with a mass of 125 GeV, the same analysis is also used to search for heavy neutral Higgs bosons with masses between 125 and 300 GeV, since similar decays are also possible for potential non-SM scalar bosons.
In the SM, a similar signature to the signal process arises when the Higgs boson decays via H → Zγ → ννγ, which has a branching fraction of 3 × 10 −4 . Searches for the decay H → Zγ using Z → final states have yielded an upper limit at 95% CL on the product of the cross section and branching fraction of about four times the SM expectation [18,19]. With the available data set, the present search is not sensitive to this SM decay, but because of enhancements that may arise from BSM physics, the search may be sensitive to Higgs boson decays to invisible particles and photons. The analysis summarized in this paper uses proton-proton (pp) collision data collected at √ s = 13 TeV by the CMS detector in 2016-18 with a total integrated luminosity of 137 fb −1 . A similar search was performed by the CMS Collaboration using the data collected at √ s = 8 TeV [20], although that analysis investigated Higgs bosons produced both in gluon-gluon fusion and in association with a Z boson.
The main backgrounds in this analysis arise from WZ and ZZ production, where an electron is mis-identified as a photon, or where additional leptons are not identified because they fail either the lepton identification criteria or the kinematic selections. A second set of backgrounds are due to WW and top quark production, where the invariant mass of the lepton pair falls into the Z boson mass window. There are also small contributions from other multiboson production processes, such as Zγ. To enhance the discrimination between the potential signal and the remaining background processes, a binned maximum-likelihood fit to several signal and control regions is performed.

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 -2 -JHEP10(2019)139 pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter, 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 magnetic flux-return yoke outside the solenoid. 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. [21]. Events of interest are selected using a twotiered trigger system [22]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events, while the second level selects events by running a version of the full event reconstruction software optimized for fast processing on a farm of computer processors.

Data samples and event reconstruction
The data used in this search were collected in three separate LHC operating periods in 2016, 2017, and 2018. The three data sets are analyzed independently, with appropriate calibrations and corrections to take into account the different LHC running conditions and CMS detector performance.
Monte Carlo (MC) simulated events are used to model the expected signal and background yields. Three sets of simulated events for each process are needed to match the different data taking conditions in the three different years. The next-to-leading-order (NLO) powheg v2 [23][24][25][26][27] generator is used to simulate the ZH signal process at NLO in quantum chromodynamics (QCD), as well as the tt, tW, and diboson processes. The BSM Higgs boson production cross sections as a function of the Higgs boson mass (m H ) for the ZH process are taken from refs. [28,29]. The signal samples are generated for masses of 125, 200, and 300 GeV. Production of ttW, ttZ, ttγ, and triple vector boson (VVV) events is generated at NLO in QCD using the MadGraph5 amc@nlo 2.2.2 (2.4.2) generator for 2016 (2017 and 2018) [30][31][32] samples. The NNPDF 3.0 NLO [33] (NNPDF 3.1 next-tonext-to-leading-order [34]) parton distribution functions (PDFs) are used for simulating all 2016 (2017 and 2018) samples. For all processes, the parton showering and hadronization are simulated using pythia 8. 226 (8.230) in 2016 (2017 and 2018) [35]. The modeling of the underlying event is generated using the CUETP8M1 [36,37] and CP5 tunes [38] for simulated samples corresponding to the 2016 and 2017-18 data sets, respectively.
All MC simulation events are processed through a simulation of the CMS detector based on Geant4 [39] and are reconstructed with the same algorithms as used for data. Additional pp interactions in the same and nearby bunch crossings, referred to as pileup, are also simulated. The distribution of the number of such interactions in the simulation is adjusted to match the one observed in the data. The average number of pileup interactions was 23 (32) in 2016 (2017 and 2018).
For this search, collision events were collected using single-electron and single-muon triggers that require the presence of an isolated lepton with transverse momentum (p T ) larger than 24 and 27 GeV, respectively. In addition, a set of dilepton triggers with lower p T thresholds were used, ensuring a trigger efficiency above 99% for events that would satisfy the subsequent offline selection.

JHEP10(2019)139
Information from all subdetectors is combined and used by the CMS particle-flow (PF) algorithm [40] for particle reconstruction and identification. Jets are reconstructed by clustering PF candidates using the anti-k T algorithm [41] with a distance parameter 0.4. Jets are calibrated in the simulation, and separately in data, accounting for energy deposits of neutral particles from pileup and any nonlinear detector response [42]. Jets with p T > 30 GeV and |η| < 4.7 are considered in the analysis. The effect of pileup is mitigated through a charged-hadron subtraction technique, which removes the energy of charged hadrons not originating from the primary interaction vertex (PV) [43]. The PV is defined as the vertex with the largest value of summed physics-object p 2 T . Here, the physics objects are the jets clustered using the jet finding algorithm [41,44] with the tracks assigned to the vertex as inputs, and the associated missing transverse momentum, taken as the negative vector p T sum of those jets.
Events are discarded if they contain a jet with p T > 20 GeV and |η| < 2.4 that is consistent with the fragmentation of a b quark. The combined secondary vertex (CSVv2) b tagging algorithm [45] is used for the 2016 data set, while the DeepCSV algorithm [45] is used for the 2017 and 2018 data sets. For the chosen working points, the efficiency to select b quark jets is about 62 (72)% for CSVv2 (DeepCSV) and the rate for incorrectly b tagging jets originating from the hadronization of gluons or u, d, s quarks is about 1%.
The vector p miss T is defined as the negative vector p T sum of all PF particle candidates. The magnitude of p miss T is the missing transverse momentum p miss T . Corrections to jet energies due to detector response are propagated to p miss T [46]. Events with possible contributions from beam halo processes or anomalous signals in the calorimeters are rejected using dedicated filters [46].
Electrons and muons are reconstructed by associating a track reconstructed in the tracking detectors with either a cluster of energy in the ECAL [47] or a track in the muon system [48]. Electron and muon candidates must pass certain identification criteria to be further selected in the analysis. For the "loose" identification, they must satisfy p T > 10 GeV and |η| < 2.5 (2.4) for electrons (muons). At the final stage of the lepton selection, the medium working points, following the definitions provided in refs. [47,48], are chosen for the identification criteria, including requirements on the impact parameter of the candidates with respect to the PV and their isolation with respect to other particles in the event [49].
Finally, photon candidates are reconstructed from energy deposits in the ECAL [50] within |η| < 2.5. The identification of the candidates is based on shower shape and isolation variables, and the medium working point, as described in ref. [50], is chosen to select those candidates. In addition, a standard "conversion-safe electron veto" [50] is applied to reject electrons misidentified as photons.

Event selection
The signal topology consists of two oppositely charged same-flavor high p T isolated leptons, electrons or muons, compatible with a Z boson decay, large p miss T , an isolated high p T photon, and little jet activity. The signal cross section is several orders of magnitude lower -4 -JHEP10(2019)139 than that of the major reducible background processes, and therefore a stringent selection is required to obtain a sample of sufficient purity. To be consistent with the expected topology, the selection requires a leading (subleading) lepton with p T > 25 (20) GeV and at least one photon with transverse momentum p γ T > 25 GeV. To reduce background processes where the lepton pair is not from the decay of a Z boson, the dilepton mass must be compatible with that of a Z boson within 15 GeV of the pole mass m Z [51]. For the purpose of rejecting the bulk of the Zγ background as well as processes with little or moderate boost, a p miss T greater than 110 GeV and a transverse momentum of the dilepton system p T larger than 60 GeV are required.
To reduce the background from WZ events with a third lepton from the W boson decay, events are removed if, in addition to the two leptons satisfying the full selection criteria, there are any loosely identified leptons. To suppress the top quark background, events are rejected if any jet passes the b tagging selection (b jet veto) described above, or if there are more than two identified jets in the event.
The signal topology is characterized by a dilepton system ( − → ) with large p T balanced in the transverse plane by the p miss T + p γ T system from the Higgs boson decay. Therefore, to reject most of the background from Zγ events with misreconstructed p miss T , the azimuthal angle between the − → and p miss ) is required to be greater than 2.5 rad, the quantity |p is required to be smaller than 0.4, and the azimuthal angle between the leading jet and p miss T (∆φ jet, p miss T ) should be greater than 0.5 rad. In addition, the mass of the dilepton and photon system (m γ ) must be greater than 100 GeV to reject resonant Zγ events, where the photon originates from final-state radiation. Finally, the transverse mass of the p miss T and photon system, defined as m T , must be smaller than 350 GeV, which rejects events where the dilepton and photon objects are weakly correlated, or where the photon momentum is mismeasured. The quantity ∆φ p miss T , p γ T is the azimuthal angle between p miss T and the photon. A summary of the selection for the analysis is shown in table 1.

Background estimation
A combination of methods based on control samples in data and simulation is used to estimate background contributions. Uncertainties related to the theoretical and experimental predictions are taken into account, as described in section 7. Background contributions are categorized depending on whether they produce at least one lepton pair from the decay of a Z boson (resonant contributions) or no such lepton pair (nonresonant contributions). The expected yield for the irreducible background from pp → Z(→ )H(→ Zγ) → ννγ is less than 0.1 events and is consequently ignored in the analysis.

Nonresonant dilepton backgrounds
The contribution from the nonresonant dilepton backgrounds, mostly WW and top quark processes, is estimated by exploiting the lepton flavor symmetry in the final states of these processes [52]. A control region based on the e ± µ ∓ final state, whose branching -5 -JHEP10 (2019)  fraction is twice that of either same flavor lepton pair final state, is used to estimate these backgrounds in the e + e − and µ + µ − channels. This region is completely dominated by this nonresonant dilepton background. The method considers the differences between the electron and muon identification efficiencies when extrapolating from the different-flavor to the same-flavor final states. The resulting predictions agree with the number of background events estimated by applying the same method to the simulation. The chosen eµ control region contains 3 events that satisfy the full analysis selection, to be compared with an expectation of 2.8 ± 0.5 (stat) from the simulation.

Resonant background with genuine missing transverse momentum
The resonant background with genuine missing transverse momentum in which an electron is mis-identified as a photon is dominated by the WZ → eν process. In this case, the background comes from events where the electron from the W boson decay is wrongly identified as a photon. The electron to photon misidentification rate is measured in Z → ee events by comparing the ratio of eγ to ee pairs consistent with the Z boson mass, as reconstructed in data and simulation. The average misidentification rate is 1-5%, with the larger values corresponding to higher photon pseudorapidity |η γ |.
Background processes with two leptons and a genuine hard photon are estimated with the simulation. These events arise primarily from the WZ → ν process (where the lepton from the W boson decay is not identified) and ZZ → 2 2ν. In both cases an additional hard photon must be radiated.
To assess the normalization of the WZ → ν background, a control region is selected in data by applying the full selection on events where the selected lepton from the W boson decay plays the role of the photon. In this region, 231 events are observed in data, while the simulation predicts 241 ± 4 (stat) events.

Resonant background with no genuine missing transverse momentum
The background from Zγ or Z + jets events is predicted by the simulation to be less than 5% of the total background, because of the stringent selection used. One of the data control regions used to verify that the background is correctly estimated selects events with lower p miss T than the default selection. Within the uncertainties, good agreement between data and simulation is found. To estimate the overall ZZ normalization, and also to emulate the Zγ process, a four-lepton sample is selected in data, and the full analysis selection is performed, with one of the Z boson dilepton pairs as a photon. In this control region, 5.1 ± 0.2 (stat) events are expected from simulation, while 7 events are observed in data.

Signal extraction
After applying the event selection, the 2016, 2017, and 2018 data sets are treated individually in order to maximize the sensitivity of the combination, since the signal-to-background ratio is different in each case. On the other hand, the electron and muon channels are merged because they show a similar signal-to-background ratio.
To discriminate between the potential signal and the remaining background processes, a binned maximum-likelihood fit to the m T spectrum is performed. The signal spectrum shows a Jacobian peak with an end-point at m T ∼ m H , while the background processes have either a flat distribution or display an increase towards lower values of m T . Since the contamination from electrons misidentified as photons is larger at large |η γ | values, improved sensitivity is achieved by considering separately events with the selected photon at low-or high-|η γ |.
In the maximum-likelihood fit, each bin of the m T distribution is separated into a low-|η γ | (|η γ | < 1) and a high-|η γ | (|η γ | > 1) bin, for the signal region and the eµ, WZ, and ZZ control regions. For each individual bin, a Poissonian likelihood term is used to describe the fluctuation of the yields around the expected central value, which is given by the sum of the contributions from signal and background processes. The uncertainties affect the overall normalizations of the signal and background yields, as well as the shapes of the predictions across the distributions of the observables. Correlations between systematic uncertainties in different categories are taken into account. Uncertainties that purely affect the normalization within a category are incorporated as nuisance parameters with log-normal priors, while those associated with changes in shapes are assigned probability density functions described by a polynomial interpolation with a Gaussian constraint. The total likelihood is defined as the product of the likelihoods of the individual bins and the probability density functions for the nuisance parameters, including the product of the likelihood for the individual years. In summary, the maximum-likelihood fit function, L, can be written as: where i runs over the three data-taking periods, j runs over the m T bins, k runs over the two |η γ | values, P(N | λ) is the Poisson probability, θ are nuisance parameters for the systematic uncertainties, andθ are their default values. The values N SR Obs,(i,j,k) , N eµ Obs,(i,k) , N 3 Obs, (i,k) , and N 4 Obs,(i,k) are the observed data events in the signal region, and the eµ, WZ, and ZZ control regions, respectively. The parameters µ, µ Nonres , µ WZ , and µ ZZ are the fit normalization factors for the signal, nonresonant, WZ, and ZZ processes, respectively. The values N ZH , N Nonres , N WZ , N ZZ , and N Other are the expected number of events for the signal, nonresonant, W, ZZ, and remaining processes, respectively, for the different regions. This approach follows that of ref.
[53], where more details can be found.
The m T distributions for the eµ, WZ, and ZZ control regions are shown in figure 2. This analysis fits the m T distributions for two regions of |η γ |, a procedure that improves the expected limits by about 30 to 50% compared with the results from simply counting the contents of a single m T bin for each |η γ | region, as was done in ref. [20]. The improvement from splitting the data into two regions of |η γ | is about 4%.

Efficiencies and systematic uncertainties
Several sources of systematic uncertainty are taken into account in the maximum-likelihood fit. For each source of uncertainty, the effects on the final distributions are considered correlated.
The assigned uncertainties in the integrated luminosity measurements for the data used in this analysis are 2.5, 2.3, and 2.5% for the 2016, 2017, and 2018 data samples [54-56], respectively. They are treated as uncorrelated across the three data sets.
The simulation of pileup events assumes a total inelastic pp cross section of 69.2 mb, with an associated uncertainty of 5% [57,58], which has an impact on the expected signal and background yields of about 1%.
Discrepancies in the lepton and photon reconstruction and identification efficiencies between data and simulation are corrected by applying scale factors to all MC simulation samples. These scale factors are determined using Z → events in the Z boson peak region that were recorded with unbiased triggers [47,48]. The scale factors depend on the p T and η of the lepton, and are within 2% of unity for both electrons and muons. The uncertainty in the determination of the trigger efficiency leads to an uncertainty smaller than 1% in the expected signal yield. The lepton momentum scale uncertainty is computed by varying the momenta of the leptons in simulation by their uncertainties, and repeating the analysis selection. The resulting yield uncertainties are ≈1% for both electrons and muons. The above procedure is applied also to determine the scale factors for photons, and the yield uncertainty for photon candidates is ≈1.5%. These uncertainties are treated as correlated across the three data sets.
The uncertainty in the calibration of the jet energy scale directly affects the acceptance of the jet multiplicity requirement and the p miss by shifting the jet energy in the simulation up and down by one standard deviation. The uncertainty in the jet energy scale is 2-5%, depending on p T and η [42], and the impact on the expected signal and background yields is about 3%. In this analysis, b tagging is used to reject events with genuine b quark jets in the final state, since signal events have no b quarks to first order in the decay channel of interest. The b tagging efficiency in the simulation is corrected using scale factors determined from data [45]. These values are estimated separately for correctly and incorrectly identified jets. Each set of values results in the b tagging efficiency of about 1-4%, and the impact on the expected signal and background yields is about 1%. The uncertainties in the jet energy scale and b tagging are treated as uncorrelated across the three data sets.
The theoretical uncertainties due to the choice of the QCD renormalization and factorization scales are estimated by varying these scales independently up and down by a factor of two [59,60]. The variations of the PDF set and the strong coupling constant are used to estimate the corresponding uncertainties in the yields of the signal and background processes, following refs. [33,61]. The combined impact on the expected yields from the -9 - ZH 300 (product of acceptance and efficiency) 3.9 ± 0.2 (10.20 ± 0.51%) above sources is about 4%. The statistical uncertainty associated with the limited number of simulated events is also considered as part of the systematic uncertainty, leading to an impact on the expected yields of about 5%. These systematic uncertainties are much smaller than the statistical uncertainty because of the limited size of the data sample, and the effect of all systematic uncertainties reduces the sensitivity by less than 4%.

Results
The numbers of observed and expected events after applying the full selection requirements are shown in table 2. The signal size is chosen for illustration purposes to be 0.1σ ZH , to have a rounded number relatively close to the minimum where this analysis is expected to have sensitivity, and to avoid quoting large expected yields. The product of acceptance and selection efficiency increases at larger m H values because of the larger p T values for all objects in the events. The m T distributions for events with |η γ | < 1 and |η γ | > 1 after the event selection are shown in figure 3. Agreement between the data and the background-only prediction is observed. By using the fit strategy described in section 6, upper limits as a function of m H are derived for the product of σ ZH and B(H → invisible + γ). For m H = 125 GeV, this result can be interpreted as an upper limit on B(H → invisible + γ) assuming the production rate for an SM Higgs boson [29]. The upper limits at 95% CL are calculated using a modified frequentist approach with the CL s criterion [62,63] and asymptotic method for the test statistic [53,64]. The observed (expected) 95% CL upper limit at m H = 125 GeV on B(H → invisible+γ) is 4.6 (3.6 +2.0   range from ∼40 to ∼4 fb as m H increases from 125 to 300 GeV. These limits also apply to other models where a scalar particle decays to a photon and light invisible particles.

Summary
A search is presented for a Higgs boson produced in association with a Z boson and decaying to an undetected particle together with an isolated photon. The analysis is based on a data set recorded by the CMS experiment in 2016-18 at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 137 fb −1 . No significant excess of events above the expectation from standard model backgrounds is found. The results are used to place limits on the product of the cross section for associated ZH production and the branching -11 -

JHEP10(2019)139
fraction for such decays of the Higgs boson, in the context of a theoretical model where the undetected particle is a massless dark photon. The observed and expected upper limits at 95% confidence level at m H = 125 GeV on B(H → invisible + γ), assuming standard model ZH associated production, are 4.6 and 3.6%, respectively. Allowing for deviations from standard model ZH production, the product of σ ZH and B(H → invisible + γ) is excluded above ∼40 to ∼4 fb, for m H ranging from 125 to 300 GeV. These are the first limits on Higgs boson decays to final states that include an undetected massless dark photon.

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: BMBWF and FWF ( Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. [51] Particle Data Group collaboration, Review of particle physics, Phys. Rev . D 98 (2018) 030001 [INSPIRE].
[57] ATLAS collaboration, Measurement of the inelastic proton-proton cross section at √ s = 13 TeV with the ATLAS detector at the LHC, Phys. Rev. Lett. 117 (2016)