Search for associated production of dark matter with a Higgs boson decaying to bb or γγ at √ s = 13 TeV

: A search for dark matter is performed looking for events with large missing transverse momentum and a Higgs boson decaying either to a pair of bottom quarks or to a pair of photons. The data from proton-proton collisions at a center-of-mass energy of 13 TeV, collected in 2015 with the CMS detector at the LHC, correspond to an integrated luminosity of 2.3 fb 1 . Results are interpreted in the context of a Z฀-two-Higgs-doublet model, where the gauge symmetry of the standard model is extended by a U (1) Z group, with a new massive Z฀ gauge boson, and the Higgs sector is extended with four additional Higgs bosons. In this model, a high-mass resonance Z฀ decays into a pseudoscalar boson A and a light SM-like scalar Higgs boson, and the A decays to a pair of dark matter particles. No significant excesses are observed over the background prediction. Combining results from the two decay channels yields exclusion limits in the signal cross section in the m Z − m A phase space. For example, the observed data exclude the Z฀ mass range from 600 to 1860 GeV, for Z฀ coupling strength g Z = 0.8, the coupling of A with dark matter particles g χ = 1, the ratio of the vacuum expectation values tan β = 1, and m A = 300 GeV. The results of this analysis are valid for any dark matter particle mass below 100 GeV. Abstract: A search for dark matter is performed looking for events with large missing transverse momentum and a Higgs boson decaying either to a pair of bottom quarks or to a pair of photons. The data from proton-proton collisions at a center-of-mass energy of 13 TeV, collected in 2015 with the CMS detector at the LHC, correspond to an integrated luminosity of 2.3 fb − 1 . Results are interpreted in the context of a Z ′ -two-Higgs-doublet model, where the gauge symmetry of the standard model is extended by a U(1) Z ′ group, with a new massive Z ′ gauge boson, and the Higgs sector is extended with four additional Higgs bosons. In this model, a high-mass resonance Z ′ decays into a pseudoscalar boson A and a light SM-like scalar Higgs boson, and the A decays to a pair of dark matter particles. No signiﬁcant excesses are observed over the background prediction. Combining results from the two decay channels yields exclusion limits in the signal cross section in the m Z ′ m A phase space. For example, the observed data exclude the Z ′ mass range from 600 to 1860 GeV, for Z ′ coupling strength g Z ′ = 0 . 8, the coupling of A with dark matter particles g χ = 1, the ratio of the vacuum expectation values tan β = 1, and m A = 300 GeV. The results of this analysis are valid for any dark matter particle mass below 100 GeV.

The CMS collaboration 29

Introduction
Astrophysical observations have provided strong evidence for the existence of dark matter (DM) in the universe [1]. However, its underlying nature remains unknown and cannot be accommodated within the standard model (SM). The recent discovery of a Higgs boson with mass of about 125 GeV by the ATLAS and CMS experiments [2][3][4] provides an additional handle to probe the dark sector beyond the SM. As explained below, in the analyses presented here, it is assumed that there are five physical Higgs bosons, and that the new state corresponds to the light neutral CP-even state h. If DM has origin in particle physics, and if other than gravitational interactions exist between DM and SM particles, DM particles (χ) could be produced at the CERN LHC. One way to observe DM particles would be through their recoil against a SM particle X (X = g, q, γ, Z, W, or h) that is produced in association with the DM. This associated production of DM and SM particles is often referred to as mono-X production. The SM particle X can be emitted directly from a quark or gluon as initial-state radiation, or through a new interaction between DM and SM particles, or as final-state radiation. The Higgs boson radiation from an initial-state quark or gluon is suppressed through Yukawa or loop processes, respectively. A scenario in -1 - which the Higgs boson is part of the interaction producing the DM particles gives monoh searches a uniquely enhanced sensitivity to the structure of couplings between the SM particles and the dark matter [5][6][7]. At the LHC, searches for DM in the mono-h channel have been performed by the ATLAS Collaboration using data corresponding to integrated luminosities of 20 fb −1 at √ s = 8 TeV and 3.2 fb −1 at √ s = 13 TeV, through the decay channels h → bb [8,9] and h → γγ [10].
In this paper, a search for DM is presented in the mono-h channel in which the Higgs boson decays to either a pair of bottom quarks (bb) or photons (γγ). The results have been interpreted using a benchmark "simplified model" recommended in the ATLAS-CMS Dark Matter Forum, which is described in ref. [11]: a Z ′ -two-Higgs-doublet-model (Z ′ -2HDM) [7], where a heavy Z ′ vector boson is produced resonantly and decays into a SM-like Higgs boson h and an intermediate heavy pseudoscalar particle A, which in turn decays into a pair of DM particles, as shown in figure 1.
In the Z ′ -2HDM model, the gauge symmetry of the SM is extended by a U(1) Z ′ group, with a new massive Z ′ gauge boson. A Type-2 2HDM [12,13] is used to formulate the extended Higgs sector. A doublet Φ u couples only to up-type quarks, and a doublet Φ d couples to down-type quarks and leptons. Only Φ u and right-handed up-type quarks u R have an associated charge under the U(1) Z ′ group, while Φ d and all other SM fermions are neutral. After electroweak symmetry breaking, the Higgs doublets attain vacuum expectation values v u and v d , resulting in five physical Higgs bosons: a light neutral CPeven scalar h, assumed to be the observed 125 GeV Higgs boson, a heavy neutral CP-even scalar H, a neutral CP-odd scalar A, and two charged scalars H ± . The analysis in this paper is performed in the context of the so-called alignment limit where the h has SM-like couplings to fermions and gauge bosons, and the ratio of the vacuum expectation values tan β = v u /v d > 0.3, as implied from the perturbativity limit of the Yukawa coupling [7,14] of the top quark, the h-H mixing angle α is related to β by α = β − π/2.
The benchmark model is parametrized through six quantities: (i) the pseudoscalar mass m A , (ii) the DM mass m χ , (iii) the Z ′ mass m Z ′ , (iv) tan β, (v) the Z ′ coupling strength g Z ′ , and (vi) the coupling constant between the A and DM particles g χ .
Only the masses m A and m Z ′ affect the kinematic distributions of the objects in the final states studied in this analysis. In fact, when A is on-shell, i.e. m A > 2m χ , the distri--2 -butions have little dependence on m χ . The remaining parameters modify the production cross section of Z ′ , branching fraction, and decay widths of the Z ′ and the A, resulting in only small changes to the final-state kinematic distributions. This paper considers a Z ′ resonance with mass between 600 and 2500 GeV and an A with mass between 300 and 800 GeV, while the mass of DM particles m χ is less than or equal to 100 GeV. The parameters tan β and g χ are fixed at unity and two different assumptions on g Z ′ are evaluated as described in more detail later. Values of m A below 300 GeV are excluded by constraints on flavor changing neutral currents from measurements of b → sγ [13], and are not considered here.
The branching fraction for decays of A to DM particles, B(A → χχ), decreases as m χ increases; for the range of m A considered in this paper, the relative decrease of B(A → χχ) is less than 7% as m χ increases from 0 to 100 GeV. Therefore, although signals with m χ = 100 GeV are considered in this search, the results are valid for any value of dark matter particle mass below 100 GeV.
The results presented here consider only A decays to DM particles and the final signal cross section σ(Z ′ → Ah → χχh) includes the value of B(A → χχ). With the assumed dark matter particle mass, the value of B(A → χχ) is ≈ 100% for m A = 300 GeV. The branching fraction starts to decrease for m A greater than twice the mass of the top quark as the decay A → tt becomes kinematically accessible. For example, if m A = 400 (800) GeV, B(A → χχ) reduces to 54 (42)%.
The quantity p miss T , calculated as the negative vectorial sum of the transverse momentum (p T ) of all objects identified in an event, represents the total momentum carried by the DM particles. The magnitude of this vector is referred to as p miss T . For a given value of m Z ′ , the p T of the A decreases as m A increases. Therefore, the p miss T spectrum softens with increasing m A . A comparison of the p miss T distributions for three values of m A is shown in figure 2.
The signal cross section is calculated for two assumptions on g Z ′ : (i) a fixed value of g Z ′ = 0.8, as considered in ref. [9] and recommended in ref. [11], and (ii) using the maximum value from electroweak global fits and constraints from dijet searches [7]: yielding g Z ′ = 0.485 for m Z ′ = 1 TeV, and g Z ′ = 0.974 for m Z ′ = 2 TeV. It can be seen from eq. (1.1) that g Z ′ = 0.8 is the maximum allowed value of g Z ′ for tan β = 1 and m Z ′ = 1.7 TeV (the best reach of LHC as estimated by ref. [7]). Note that this analysis does not consider the contribution of another decay that gives a similar mono-h signature: , is a function of tan β and m Z ′ and does not depend on g Z ′ since the value of g Z ′ cancels in the ratio. The h → bb decay mode has the largest branching fraction (≈58%) of all, but suffers from relatively poor mass resolution of about 10%, and while the h → γγ branching fraction is small (≈0.2%), the channel benefits from the high precision in reconstructed diphoton mass, with a resolution of about 1-2%.  In the h → bb channel, the fact that the p T of the h should increase with m Z ′ and decrease with m A is exploited. The minimum separation in the pseudorapidity and azimuth (η, φ) plane between the decay products of h scales as m h /p h T , where p h T is the transverse momentum of the h boson. The allowed mass ranges of m Z ′ and m A imply a very wide range of values for p h T and consequently a wide range in the separation of the decay products. Analysis in this channel is therefore divided into two regimes: (i) a resolved regime where the h decays to two distinct reconstructed b jets, and (ii) a Lorentzboosted regime where the h is reconstructed as a single fat jet. For each mass point, the analysis with best sensitivity for the expected limit is used as the final result. The signal extraction is performed through a simultaneous fit to the signal-and background-enriched control regions.
The search in the h → γγ channel is performed by seeking an excess of events over the SM prediction in the diphoton mass spectrum, after requiring a large p miss T . Control samples in data are used to estimate the reducible background, which mainly consists of diphoton SM production. A counting approach is used to estimate the potential signal.
The paper is organized as follows. After a brief introduction to the CMS detector in section 2, the data and simulated events used for the analysis are described in section 3. The event reconstruction is detailed in section 4. Section 5 describes the analysis strategy for both Higgs boson decay channels. The description of the most relevant systematic uncertainties affecting the analysis is found in section 6. Finally, the results of the search are reported in section 7, and the summary is presented in section 8.

The CMS detector
The central feature of the CMS detector is a superconducting solenoid, of 6 m internal diameter, providing an axial magnetic field of 3.8 T along the beam direction. Within the -4 -solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL). Charged particle trajectories are measured by the silicon pixel and strip tracker system, covering 0 ≤ φ ≤ 2π in azimuth and |η| < 2.5. The electromagnetic calorimeter, which surrounds the tracker volume, consists of 75,848 lead tungstate crystals that provide coverage in pseudorapidity |η| < 1.48 in the barrel region (EB) and 1.48 < |η| < 3.0 in two endcap regions (EE). The EB modules are arranged in projective towers. A preshower detector consisting of two planes of silicon sensors interleaved with a total of three radiation lengths of lead is located in front of the EE. In the region |η| < 1.74, the HCAL cells have widths of 0.087 in pseudorapidity and azimuth. In the (η, φ) plane and for |η| < 1.48, the HCAL cells map on to 5×5 ECAL crystal arrays to form calorimeter towers projecting radially outwards from the nominal interaction point. For |η| > 1.74, the coverage of the towers increases progressively to a maximum of 0.174 in ∆η and ∆φ. Extensive forward calorimetry complements the coverage provided by the barrel and endcap calorimeters. Muons are measured in gas-ionization detectors embedded in the steel return yoke. A more detailed description of the CMS detector can be found in ref. [15].

Data and simulated samples
The analysis is performed with pp collision data at √ s = 13 TeV collected by the CMS experiment at the LHC during 2015, corresponding to an integrated luminosity of 2.3 fb −1 . The MadGraph5 amc@nlo v2.3.0 [16] generator is used to generate the mono-h signal at leading order (LO) as predicted by the Z ′ -2HDM model described in section 1. In the MadGraph5 amc@nlo generation, a vector particle Z ′ that decays to a SM-like Higgs boson h with mass 125 GeV is produced resonantly together with a heavy pseudoscalar particle A that decays into a pair of DM particles. The decay of the SM-like Higgs boson is handled by pythia 8.205 [17].
The associated production of a SM Higgs boson and a Z boson (Zh) is a small but irreducible background for both decay channels. The Vh (Zh and Wh) processes are simulated using powheg v2.0 [18,19] and MadGraph5 amc@nlo for qq and gluon-gluon fusion, respectively. In the h → γγ channel, additional resonant but reducible backgrounds are considered. These backgrounds include the SM Higgs boson, produced through gluon fusion (ggh), through vector boson fusion (VBF), and in association with top quarks (tth). All of these resonant backgrounds are modeled at next-to-leading order (NLO) in simulation. The VBF Higgs boson samples are generated using powheg [20], while the ggh and tth samples are generated with MadGraph5 amc@nlo.
The dominant background processes for the h → bb decay channel are events with top quarks and W/Z bosons produced in association with jets. The tt events, produced via the strong interaction, and electroweak production of single top quarks in the t-and tW-channels are generated at NLO with powheg [21][22][23][24][25]. The s-channel process of single top quark production is generated with MadGraph5 amc@nlo. Differential measurements of top quark pair production show that the measured p T spectrum of top quarks is softer than the one produced in simulation. Scale factors to correct for this effect are -5 -derived from previous CMS measurements [26,27]. The sum of top quark pair events and single top quark events is referred to as "Top quark background" in the rest of the paper. The W and Z boson production in association with jets is simulated at LO with MadGraph5 amc@nlo. Up to four additional partons in the matrix element calculations are included. The MLM matching scheme [28] is used as an interface to the parton shower generated with pythia. The cross sections for W+jets and Z+jets processes are normalized to the next-to-next-to-leading order cross section, computed using fewz v3.1 [29]. Moreover, to improve the description of the distribution of high p T W+jets and Z+jets processes, events are reweighted using the generated p T of the vector boson to account for NLO quantum chromodynamics (QCD) and electroweak (EW) contributions [30][31][32]. The small background from diboson (WW, WZ, and ZZ) processes, labeled as VV in the rest of the paper, is simulated with pythia.
For the h → γγ decay channel, several nonresonant background sources can mimic the signal when an event has mismeasured p miss T and two photons with an invariant mass close to the mass of the SM-like Higgs boson. These sources include contributions from dijet and multijet events, EW processes such as t, tt, Z, ZZ, or W bosons produced in association with one or two photons, γγ, γ+jet, and Drell-Yan (DY) production in association with jets, where the Z boson decays to pairs of electrons and neutrinos. These backgrounds are generated with MadGraph5 amc@nlo, with the exception of the ZZ sample, which is generated with powheg [33]. These nonresonant background samples are not used for the background estimation, but are used to optimize the selection.
All simulated samples use the NNPDF 3.0 PDF sets [34]. The parton showering and hadronization are performed with pythia using the CUETP8M1 tune [35,36]. For the h → bb decay channel, to perform systematic studies in the boosted regime, an additional signal sample is generated with MadGraph5 amc@nlo, parton-showered and hadronized by herwig++ v2.7.1 [37] using the UE-EE-5C tune [38,39]. The samples are processed through a Geant4-based [40] simulation of the CMS detector. All samples include the simulation of "pileup" arising from additional inelastic proton-proton interactions in the same or neighboring bunch crossings. An average of approximately ten pileup interactions per bunch crossing is included in the simulation with a separation between bunches of 25 ns. The simulated pileup distribution is reweighted to match the corresponding observed distribution in the analyzed data.

Event reconstruction
A global event reconstruction is performed using the particle-flow (PF) [41][42][43] algorithm, which optimally combines the information from all the subdetectors and produces a list of stable particles, namely muons, electrons, photons, charged and neutral hadrons.
The reconstructed interaction vertex with the largest value of i p 2 Ti , where p Ti is the transverse momentum of the i th track associated with the vertex, is selected as the primary event vertex. This vertex is used as the reference vertex for all objects reconstructed using the PF algorithm. The offline selection requires all events to have at least one primary -6 -vertex reconstructed within a 24 cm window along the z-axis around the mean interaction point, and a transverse distance from the mean interaction region less than 2 cm.
Jets are reconstructed from the PF candidates, after removing charged hadrons originating from pileup vertices, using the anti-k T clustering algorithm [44] with distance parameters of 0.4 (AK4 jet) and 0.8 (AK8 jet), as implemented in the FastJet package [45]. In order to improve the discrimination of signal against multijet background, the pruning algorithm described in refs. [46,47], which is designed to remove contributions from soft radiation and pileup, is applied to AK8 jets. The pruned jet mass (m pruned corrected ) is defined as the invariant mass associated with the four-momentum of the pruned jet, after the application of the jet energy corrections [48]. Corrections to jet momenta are further propagated to the p miss T calculation [49]. In addition, tracks with p T > 1 GeV, |η| < 2.5, and with longitudinal impact parameter |d Z | < 0.1 cm from the primary vertex are used to reconstruct the track-based missing transverse momentum vector, p miss T,trk . The jets originating from the decay of b quarks are identified using the combined secondary vertex (CSV) algorithm [50,51], which uses PF jets as inputs. The algorithm combines the information from the primary vertex, track impact parameters, and secondary vertices within the jet using a neural network discriminator. The loose (medium) working point (WP) used in this analysis has a b jet selection efficiency of 83% (69%), a charm jet selection efficiency of 28% (20%), and a mistag rate for light-flavor jets of ≈10% (1%) [50]. The AK8 jets are split into two subjets using the soft-drop algorithm [52,53]. The CSV algorithm is tested and validated for AK4 and AK8 jets [50]. The working points for the analyses of the resolved and boosted regimes were chosen by maximizing the expected significance. The loose WP of the subjet b tagging algorithm is used for the boosted regime, whereas the medium WP of the AK4 jet b tagging algorithm is used for the resolved regime, since the background is higher in this case.
Photons are reconstructed in the CMS detector from their energy deposits in the ECAL, which come from an electromagnetic shower involving several crystals. The energy is clustered at the ECAL level by building a cluster of clusters, supercluster (SC), which is extended in the φ direction because of the strong magnetic field inside the detector, which deflects the electron and positron produced if the photon converts in the tracker [54]. In order to achieve the best photon energy resolution, corrections are applied to remove channel-to-channel response variations and to recover energy losses due to incomplete containment of the shower or conversions, as detailed in ref. [55]. Additional residual corrections are made to the measured energy scale of the photons in data (≤1%) and to the energy resolution in simulation (≤2%) based on a detailed study of the mass distribution of Z → e + e − events. The uncertainties in the measurements of the photon energy scale and resolution are taken as systematic uncertainties as described in section 6. This process is outlined for the 8 TeV data set in ref. [55]. Values are adjusted for the 13 TeV data set.
Electron reconstruction requires the matching of a supercluster in the ECAL with a track in the silicon tracker. Identification criteria [56] are based on the ECAL shower shape. Muons are reconstructed by combining two complementary algorithms [57]: one in which tracks in the silicon tracker are matched to a muon track segment, and another -7 -in which a global track fit is performed, seeded by the muon track segment. Further identification criteria are imposed on muon candidates to reduce the number of misidentified hadrons. Hadronically decaying τ leptons (τ h ) are reconstructed using the hadron-plusstrips algorithm [58], which uses the charged-hadron and neutral-electromagnetic objects to reconstruct intermediate resonances into which the τ lepton decays.

Event selection and background estimation
This analysis searches for excesses over the background-only prediction in events with large p miss T and a system of two b-tagged jets or two photons that has a reconstructed invariant mass close to the mass of the SM-like Higgs boson h. In the h → bb decay channel, the analysis relies on fitting the p miss T distribution simultaneously in the signal region (SR), defined after selecting a mass window around the Higgs boson mass, and in background-enriched control regions (CRs). For the h → γγ decay channel, a simple analysis is performed where the signal and resonant background contributions are estimated by counting the number of simulated events in the SR, while the nonresonant background is extrapolated from the data in a low-p miss T region. In the following sections, the event selection and analysis strategy are described in detail for the two channels separately.

The channel h → bb
A search for DM produced in association with h → bb is performed in a resolved regime, where events are required to have at least two AK4 jets, and in the Lorentz-boosted regime where one AK8 jet is required. In addition, p miss T is required to be large because it is a key signature of the signal events and it provides strong rejection against the large reducible backgrounds described in section 3.

Event selection
The trigger used in the selection of signal-like events requires p miss T > 90 GeV and H miss T > 90 GeV, where H miss T is defined as the magnitude of the vectorial sum of the p T of all jets in the event with p T > 20 GeV. An additional trigger with a p miss T > 170 GeV requirement is used to achieve higher efficiency. In this way, events with either high p miss T or high H miss T will pass the trigger. For events passing the selection criteria that have p miss T > 170(200) GeV for the resolved (boosted) analysis, the trigger efficiency is found to be greater than 98%. The p miss T threshold for the analysis of the resolved regime is set slightly lower to enhance the signal efficiency in this region of phase space, where the p miss T distribution is softer. Event filters are used to remove spurious high p miss T events caused by instrumental noise in the calorimeters, or beam halo muons. It has been verified that the efficiency of these filters for accepting signal events is very close to 100%. The main part of the event selection consists of Higgs boson tagging. This selection is different for the resolved and boosted analyses. In the resolved regime, events are required to have two AK4 jets with p T > 30 GeV and |η| < 2.4. These two jets are used to reconstruct the Higgs boson candidate, which is required to have p T > 150 GeV. Each of the two AK4 jets in the resolved regime is required to pass the b tagging selection, whereas in the boosted regime, -8 -the two subjets inside an AK8 jet must both pass the b tagging selection. In the boosted regime, the decay products from the Higgs boson are merged. Therefore, an AK8 jet with p T greater than 200 GeV is used to reconstruct the Higgs boson. If more than one Higgs boson candidate is reconstructed, the ambiguity is resolved by selecting the candidate with the highest p T . Backgrounds due to hadronic jets are further reduced by constraining the reconstructed Higgs boson candidate mass, m bb , to be between 100 and 150 GeV. For the resolved regime, the Higgs boson candidate mass is reconstructed using two b-tagged AK4 jets. For the boosted regime, the corrected pruned mass of the AK8 jet with two b-tagged subjets is used as the Higgs boson candidate mass.
Multijet events can act as a source of background when the energy of one of the jets is mismeasured. Therefore, the absolute difference between the azimuthal angles of the vector p miss T and any other AK4 jet with p T > 30 GeV is required to be greater than 0.4 radians. Multijet background is further reduced in the resolved analysis by requiring the azimuthal angle difference between the p miss T and p miss T,trk to be less than 0.7 radians. Events are rejected if they have any isolated electron (muon) with p T > 10 GeV and |η| < 2.5 (2.4) or any τ h candidates with p T > 20 GeV and |η| < 2.3 [56,58,59]. In addition, the events must not have any additional loose AK4 b-tagged jet or more than one additional AK4 jet with p T > 30 GeV and |η| < 4.5. These vetoes considerably reduce the background from semileptonic top decay modes and leptonic decays of W+jets.
The product of the detector acceptance and selection efficiency varies from 1 to 29%, depending on the values of m Z ′ and m A . The average p miss T increases with m Z ′ and decreases with m A . The overall selection efficiency, shown in table 1, follows the same trend.

Analysis strategy and background estimation
Several CRs are used to correct the background normalizations with dedicated scale factors. For both resolved and boosted regimes, the selection criteria of these CRs are kept as close as possible to those of the SR, except for the inversion of the additional object vetoes (leptons, jets) and the Higgs boson mass window. This makes the CRs orthogonal to the SR.
For the resolved regime, three CRs are specified: Z(→ νν)+jets, top quark, and W+jets. The b tagging selection in all the CRs is the same as in the SR in order to minimize the b tagging systematic uncertainties when extrapolating the background scale factors measured in the CRs to the SR. The Z(→ νν)+jets CR is defined with the same selection as the SR, except for the inversion of the reconstructed Higgs boson mass requirement. The W+jets and top quark CRs are defined by removing the mass selection and requiring exactly one isolated electron (muon) with p T > 10 GeV and |η| < 2.5 (2.4). Events with one additional AK4 jet are placed in the top quark CR, whereas events with no additional AK4 jets enter the W+jets CR.
For the boosted regime, the Z(→ νν)+jets CR is defined by inverting the mass requirement for the AK8 jet. Owing to the low event count and very similar topology between the W+jets and top quark backgrounds it is difficult to construct two separate CRs for W+jets and top quark backgrounds. Hence, the single-lepton CR, a combination of mainly W+jets and top quark events, is defined using the same selection as that for the signal, -9 - Zj  but requiring exactly one isolated electron (muon) with p T > 10 GeV and |η| < 2.5 (2.4) and removing the mass requirement. Figure 3 shows the Higgs boson candidate mass for the resolved and boosted regimes. They correspond to the simultaneous fit of the p miss T distributions in the SR and background enriched CRs to extract the signal. Data-to-simulation ratios for pre-fit and post-fit background predictions are shown in the lower panels of all figures 3-6. Figure 4 shows the comparison of data and simulation for the main observable, p miss T , in the W+jets, top quark, and Z(→ νν)+jets CRs for the resolved regime. The comparison between data and simulated samples for the boosted regime is shown in figure 5 for the single-lepton CR and the Z(→ νν) mass sideband region. Figure 6 shows the p miss T distributions in three bins in the SR that are used for the final signal extraction. These three bins were chosen to optimize the expected limits. The selected signal and background events are compared to data and fit simultaneously in the SR and CRs in three p miss T bins, separately for the resolved and the boosted regimes. The simultaneous fit of SR and background-enhanced CRs is performed correlating the scale factors and systematic uncertainties as described in section 6. The measured datato-simulation post-fit scale factors are compatible with unity within the total combined statistical and systematic uncertainty. In particular, for the resolved regime, the scale factors for the backgrounds are 1.23 ± 0.17 for Z(→ νν)+jets, 1.    Pre-fit Post-fit Figure 6. Post-fit distribution of p miss T expected from SM backgrounds and observed in data for the resolved (left) and the boosted (right) regimes in the signal region with three different m Z ′ signal points overlaid. Other parameters for this model are fixed to m χ = 100 GeV and tan β = g χ = 1. The cross sections for the signal models are computed assuming g Z ′ = 0.8. The bottom panels show the data-to-simulation ratios for pre-fit (red markers) and post-fit (black markers) background predictions with a hatched band corresponding to the uncertainty due to the finite size of simulated samples and a gray band that represents the systematic uncertainty in the post-fit background prediction (see section 6). The last bin includes all events with p miss T > 350 (500) GeV for the resolved (boosted) regime.

The channel h → γγ
The h → γγ search is performed using a diphoton selection. A set of requirements is applied to ensure good-quality photon candidates. Additional kinematic requirements on the objects in the final state are applied to reduce the background. The diphoton invariant mass and p miss T are used as the discriminating variables to estimate the signal.

Event selection
Diphoton triggers with asymmetric transverse energy thresholds (30/18 GeV) are used to select events with the diphoton invariant mass above 95 GeV. The trigger selection uses a very loose photon identification based on the cluster shower shape and loose isolation requirements (both defined in detail in ref. [55]), and a requirement that the ratio of hadronic-to-electromagnetic energy of the photon candidates is less than 0.1.
The main source of background for photons, which arises from jets with high electromagnetic energy content, is rejected by considering the ratio of energies deposited by the photon candidate in the hadron and electromagnetic calorimeters and the spread of the energy deposition in the η direction, as described in [55]. In addition, misidentified photons are rejected using the isolation variables Iso Ch , Iso γ , and Iso Neu calculated by summing the p T of the charged hadrons, photons and neutral hadrons, respectively, in a cone of radius ∆R = 0.3. In the photon identification, Iso Neu and Iso γ are corrected for the median transverse energy density (ρ) of the event to mitigate the effects of pileup [60].
The photons in the EB (i.e. the photons with |η| ≤ 1.44) and photons in the EE (1.566 ≤ |η| ≤ 2.5) have different selection criteria, equivalent to those used in refs. [61,62]. The working point chosen for this analysis corresponds to 90.4% (90.0%) photon ID efficiency in the EB (EE), while the misidentification rate in the EB (EE) is 16.2% (18.7%) for objects with p T > 20 GeV.
A high-quality interaction vertex, defined as the reconstructed vertex with the largest number of charged tracks, is associated to the two photons in the event. The efficiency of selecting the correct vertex for all generated mass points, defined as the fraction of signal events with well reconstructed vertices that have a z position within 1 cm of the generator-level vertex, is approximately 78%.
The optimal signal selection is chosen by studying the discriminating power of variables such as the p T /m γγ of each photon, p miss T , and the p T of the diphoton system (p Tγγ ). A selection on p T that scales with m γγ is chosen such that it does not distort the m γγ spectrum shape. The p Tγγ variable, included because it has a better resolution than p miss T , has a distribution of values that are on average larger for signal than for background events, given that the Higgs boson is expected to be back-to-back in the transverse plane with the p miss T .
In addition, two geometrical requirements are applied to enhance the signal over background discrimination and to veto background events with mismeasured p miss T : • the azimuthal separation between the p miss T and the Higgs boson direction (reconstructed from the two photons) |∆φ(γγ, p miss T )| must be greater than 2.1 radians.
-13 - • the minimum azimuthal angle difference between the p miss T and the jet direction in the event min(|∆φ(jet, p miss T )|) must be greater than 0.5 radians. The jet direction is derived by considering all the jets reconstructed from the clustering of PF candidates by means of the anti-k t algorithm [44] with a distance parameter of 0.4. Jets are considered if they have a p T above 50 GeV in the |η| range below 4.7 and satisfy a loose set of identification criteria designed to reject spurious detector and reconstruction effects.
The set of selection criteria that maximizes the expected significance for each Z ′ mass point is studied. The optimized selection for the m Z ′ = 600 GeV and m A = 300 GeV sample maintains a large efficiency for the other signal mass points, while the backgrounds remain small. Therefore a common set of criteria is used for all signal masses with m Z ′ between 600 and 2500 GeV and m A between 300 and 800 GeV. The chosen kinematic selections include p T1 /m γγ > 0.5, p T2 /m γγ > 0.25, p Tγγ > 90 GeV, p miss T > 105 GeV. Events are vetoed if they have any muons or more than one electron present. This allows the analysis to be sensitive to events where an electron originating from conversion of the photon before reaching the ECAL is identified outside the photon supercluster. Standard lepton identification requirements are used [56,59]. This requirement is 100% efficient for the signal and reduces significantly the EW background contributions.
The SR of this analysis is defined as the region with 120 < m γγ < 130 GeV and p miss T above 105 GeV. The distribution of m γγ for the selected events before the p miss T requirement is shown in figure 7 for the full mass range considered in this analysis: 105 < m γγ < 180 GeV. Also shown is the p miss T distribution of the selected events after the m γγ SR selection. It can be seen that after applying the requirement that m γγ has to be close to the Higgs boson mass, the SM background contribution in the high-p miss T region is close to zero and the DM signal is well separated from the background distribution.

Background estimation
The final state with a γγ pair and large p miss T has two classes of background: resonant and nonresonant. The contributions from each class are treated differently.
Resonant backgrounds arise from decays of the SM Higgs boson to two photons. They appear as an additional peak under the expected signal peak and are evaluated with the MC simulation by counting the number of expected events from all SM Higgs production modes in the SR.
The contribution of the nonresonant backgrounds (N bkg SB ) in the sideband (SB) region, mostly multijets and EW processes with mismeasured large p miss T and misidentified photons, is evaluated from the data by counting the number of events in the m γγ sidebands 105 < m γγ < 120 GeV and 130 < m γγ < 180 GeV, with p miss T > 105 GeV in both cases. Then N bkg SB is scaled by a transfer factor α to take into account the relative fraction between the number of events in the m γγ SR and SB region. The expected number of nonresonant background events in the SR is given by:  The derivation of α relies on the knowledge of the background shape f bkg (m γγ ) as follows: and is evaluated by performing a fit to the m γγ distribution in a CR of the data. In this analysis, the low-p miss T CR, with p miss T < 105 GeV, is used. The fit to data in the low-p miss T region used to calculate α is shown in figure 8. In this case the negligible contribution of the resonant SM Higgs boson processes is not considered. The data are fit with a background-only model using an analytic power law function: where the parameter a, the normalization, and b are free parameters, defined as positive. The fit is performed with an unbinned maximum likelihood technique. The function defined in eq. (5.3) was chosen after examining several models and performing a bias study using nonresonant background MC to evaluate any possible background mismodeling, following the procedure described in ref. [63]. It has been verified that the fitted parameters of the power law function are compatible within the uncertainties with both data and simulation.
To derive a robust estimate of α, several fits to both data and simulated background events are performed using different analytic functions and looking at different CRs of p miss T . Within the uncertainties, α is independent of the p miss T CR used and is consistent between data and simulation. The fitted shape of the low-p miss T CR in data is taken as -15 -

Systematic uncertainties
The systematic uncertainties common to the two Higgs boson decay channels are as follows.
An uncertainty of 2.7% is used for the normalization of simulated samples in order to reflect the uncertainty in the integrated luminosity measurement in 2015 [64]. In the h → bb analysis an uncertainty of 2% is estimated in the signal yield for p miss T above 170 GeV by varying the parameters describing the trigger turn-on. For the h → γγ analysis the trigger uncertainty (approximately 1%) is extracted from Z → e + e − events using the tag-andprobe technique [65]. The following uncertainties in clustered and unclustered calorimetric energy affect the p miss T shapes and the normalization of the signal and background yield predictions: the JES for each jet is varied within one standard deviation as a function of p T and η, and the efficiency of the event selection is recomputed to assess the variation on the normalization and p miss T shape for signal and backgrounds; the signal acceptance and efficiency are recomputed after smearing the energy of each jet to correct for the difference in jet energy resolution between the data and simulation (≈5%); the systematic uncertainties in the calibration of unclustered energy in the calorimeter are propagated as normalization and shape uncertainties in the p miss T calculation. The total effect of the systematic uncertainty in the signal yield, considering all of these variations on p miss T is approximately 3% for the h → bb analysis and less than 1% for the h → γγ analysis. Among the three sources, the JES is the one that most affects the signal yield.
-16 - The following systematic uncertainties only affect the h → bb decay channel: the b tagging scale factors are applied consistently to jets in signal and background events. An average systematic uncertainty of 6% per b jet, 12% per c jet, and 15% per light quark or gluon jet is used to account for the normalization uncertainty [50]. The pruned mass distribution of the AK8 jet is not perfectly reproduced by simulation. Therefore, a control region, with a large number of events enriched in boosted hadronically decaying W bosons reconstructed as AK8 jets, is used to measure the systematic uncertainty due to this effect, giving an estimated value of 5%. Moreover, different hadronization algorithms (pythia and herwig++) give slightly different shapes for the pruned mass distribution. Therefore, an additional uncertainty of 10% is assigned to account for the difference between simulations. For the boosted regime, the same background normalization scale factor is used for W+jets and top quark backgrounds. The uncertainty in the relative normalization of these two processes is 30%. An uncertainty of 2% is measured by varying the lepton efficiency scale factors within one standard deviation and recomputing the signal selection efficiency. For W+jets, Z(→ νν)+jets and top quark backgrounds, variations in the renormalization and factorization scales directly affect the normalization and shape of the p miss T distribution. A variation of approximately 5% is found for the yields of these backgrounds in the signal region. The uncertainty in the signal acceptance and p miss T shape due to the choice of PDFs is measured following the method described by the PDF4LHC group [66]. A variation of approximately 3% is found in the signal yields. The effect of electroweak corrections as described in section 3 is studied by recomputing the normalization and shapes for the W+jets and Z(→ νν)+jets backgrounds, by alternately removing the corrections or doubling them. An uncertainty of 20% is assumed for the single top quark , SM Higgs boson, and diboson production rates. Uncertainties due to the finite size of the signal and background simulated samples are included in the normalization and shape, such that each bin of the final fitted distributions is affected independently.
In summary, for h → bb, the overall uncertainties related to background determination methods, simulation, and theory inputs are estimated to be 10% in the background contributions in the SR. The impact of the uncertainty in the major background contributions (W+jets, Z(→ νν)+jets and top quarks) in the SR is reduced by constraining the normalizations of these processes in data with the simultaneous fit of p miss T shapes in the SR and CRs. The major sources of systematic uncertainties that affect the fit are JES uncertainties, b tagging uncertainties, and the statistical uncertainty in the simulated Z(→ νν)+jets and W+jets background samples. The effect of the remaining uncertainties on the final fit is ≈1%.
The following systematic uncertainties affect only the h → γγ analysis: as shown in eq. (5.1), the predicted number of nonresonant background events in the SR is evaluated from the number of observed events in the m γγ sidebands in the high-p miss T region (N bkg SB ) multiplied by a transfer factor α obtained by fitting the m γγ distribution in the low-p miss T control region. Therefore two different systematic uncertainties are assigned to this procedure, one for N bkg SB and one for α. The first systematic uncertainty takes into account the fact that N bkg SB is statistically limited. Secondly, a 20% systematic uncertainty is assigned to reflect the imperfect knowledge of the background m γγ shape in the low-p miss T region, hence -17 -on the knowledge of the α factor. This uncertainty is obtained by performing the fit to the m γγ distribution using several analytic functions, using data rather than using simulated events, and using other p miss T CRs. An observed peak above the diphoton continuum in the m γγ distribution around the SM Higgs boson mass would have a SM h → γγ contribution. In order to extract the DM signal, the resonant background contribution has to be evaluated and subtracted. The SM Higgs boson contribution is affected by both theoretical and experimental systematic uncertainties. For each SM Higgs boson production mechanism (ggh, VBF, tth, Vh), the uncertainties on the PDFs and α s , provided in ref. [67], are addressed using the procedure from the PDF4LHC group [66]. The size of the systematic uncertainty is computed for each process and category separately by checking the effect of each weight on the final event yield. An additional uncertainty on the h → γγ branching fraction of 5% is included following ref. [67]. A 1% photon energy scale uncertainty is assigned. This number takes into account the knowledge of the energy scale at the Z boson peak and of its extrapolation to higher masses. The uncertainty on the photon resolution correction factors is evaluated by raising and lowering the estimated additional Gaussian smearing measured at the Z boson peak by 0.5% in quadrature. The photon identification uncertainty is taken as an uncertainty in the data-to-simulation scale factors, which can be as large as 2%, depending on the p T and the η of the photon.
The h → γγ decay channel results are only marginally affected by systematic uncertainties as statistical uncertainties dominate the analysis.

Results
For the event selection described in section 5, the predicted signal acceptances multiplied by the efficiencies (Aǫ) are listed in table 1 for the two decay channels. Table 2 shows, for the h → bb channel, the SR post-fit yields for each background and signal mass point along with the sum of the statistical and systematic uncertainties for the resolved and boosted regimes. The total background uncertainty is approximately 10% and mainly driven by the systematic uncertainty.
For the h → γγ channel, when applying the event selection to the data, two events are observed in the m γγ sidebands and are used to evaluate the magnitude of the nonresonant background as described in section 5.2.2. This yields an expected number of 0.38 ± 0.27 (stat) nonresonant background events in the SR. Expected resonant background contributions are taken from the simulation as detailed in section 5.2.2 and are 0.057 ± 0.006 (stat) events considering both the Vh production (dominant) and the gluon fusion mode. Zero events are observed in the SR in the data.
Since no excess of events has been observed over the SM background expectation in the signal region, the results of this search are interpreted in terms of an upper limit on the production of DM candidates in association with a Higgs boson in the process Z ′ → Ah → χχh. The upper limits are computed at 95% confidence level (CL) using a modified frequentist method (CL s ) [67][68][69] computed with an asymptotic approximation [70]. A profile likelihood ratio is used as the test statistic in which systematic uncertainties are modeled as nuisance parameters. These limits are obtained as a function of m Z ′ and m A -18 -  Table 1. The product of acceptance and efficiency (with statistical uncertainty) for signal in the SR, after full event selection for the h → bb (upper) and the h → γγ (lower) decay channels. The systematic uncertainty for h → bb (h → γγ) is approximately 10% (5%). For h → bb, the value shown here is either for the resolved regime or for the boosted regime, depending on which is used for the calculation of the limit on σ (Z ′ → Ah → χχh), as shown in figure 10 left.
for both Higgs boson decay channels and for the combination of the two. The two decay channels are combined using the branching ratios predicted by the SM. In the combination of the two analyses, all signal and p miss T -related systematic uncertainties as well as the systematic uncertainty in the integrated luminosity are assumed to be fully correlated. Figure 9 (left) shows the 95% CL expected and observed limits on the dark matter production cross section σ(Z ′ → Ah → χχh), for h → bb and h → γγ for m A = 300 GeV. These results, obtained with m χ = 100 GeV, can be considered valid for any dark matter particle mass below 100 GeV since the branching fraction for decays of A to DM particles, B(A → χχ), decreases as m χ increases. As shown in figure 9, for the phase space parameters considered for this model (g χ and tan β equal to unity), results of the combined analysis are mainly driven by the h → bb channel. The combination with the h → γγ channel provides a 2-4% improvement in terms of constraints on the model for the low Z ′ mass values. Future iterations of this search will explore additional phase space regions of the Z ′ -2HDM model, i.e. larger values of tan β, where the h → γγ channel becomes more sensitive than h → bb [7].  Table 2. Post-fit background event yields and observed numbers of events in data for 2.3 fb −1 in both the resolved and the boosted regimes for the h → bb analysis. The expected numbers of signal events for m A = 300 GeV, scaled to the nominal cross section with g Z ′ = 0.8, are also reported. The statistical and systematic uncertainties are shown separately in that order.
Figures 9 (right) and 10 show the 95% CL expected and observed upper limits on the signal strength σ 95%CL (Z ′ → Ah → χχh)/σ theory (Z ′ → Ah → χχh). For m A = 300 GeV, the Z ′ mass range from 600 to 1780 GeV is expected to be excluded with a 95% CL when the signal model cross section is calculated using g Z ′ = 0.8, while the observed data, for m A = 300 GeV, exclude the Z ′ mass range from 600 to 1860 GeV. When the signal model cross section is calculated using the constrained g Z ′ , the expected exclusion range is 830 to 1890 GeV, and the observed exclusion range is 770 to 2040 GeV. Figure 10 shows the expected and observed upper limits on the signal strength for the h → bb and h → γγ decay channels. Figure 11 shows the upper limits on the signal strength combining the results from both the h → bb and h → γγ decay channels.    Figure 9. Left: the expected and observed 95% CL limits on dark matter production cross sections for h → bb and h → γγ for m A = 300 GeV. The exclusion region is shown for two g Z ′ values. The dark green and light yellow bands show the 68% and 95% uncertainties on the expected limit. Right: the expected and observed 95% CL limits on the signal strength for m A = 300-800 GeV are shown.

Summary
Other parameters for this model are fixed to m χ = 100 GeV and tan β = g χ = 1. The theoretical cross section (σ th ) used for the right-hand plot is calculated using g Z ′ = 0.8.  Figure 10. The observed (expected) 95% CL limits on the signal strength (as in figure 9 right), separately for the h → bb (left) and h → γγ (right) decay channels, and for m A = 300-800 GeV and m Z ′ = 600-2500 GeV. Other parameters for this model are fixed to m χ = 100 GeV and tan β = g χ = 1. The theoretical cross sections are calculated using g Z ′ = 0. 8 Figure 11. The observed (expected) 95% CL limits on the signal strength (as in figure 9 right) for the combination of h → γγ and h → bb decay channels, and for m A = 300-800 GeV and m Z ′ = 600-2500 GeV. Other parameters for this model are fixed to m χ = 100 GeV and tan β = g χ = 1. The theoretical cross sections times branching fractions are calculated using g Z ′ = 0.8. decays to two dark matter candidates. Two distinct channels are studied, where the Higgs boson decays to two b quarks or two photons.
No significant deviation is observed from the standard model background. With optimized selections, limits on the signal cross section σ(Z ′ → Ah → χχh) are calculated for various values of m Z ′ and m A assuming g χ and tan β equal to one. The limits are valid for any dark matter particle mass below 100 GeV. For m A = 300 GeV, the observed data exclude the Z ′ mass range of 600 to 1860 GeV for g Z ′ = 0.8, and the range 770 to 2040 GeV for the constrained value of g Z ′ . This is the first result on a search for dark matter produced in association with a Higgs boson at √ s = 13 TeV that combines results from the h → bb and h → γγ channels.

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 centres 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: the Austrian