Search for dark matter particles produced in association with a Higgs boson in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for dark matter (DM) particles is performed using events with a Higgs boson candidate and large missing transverse momentum. The analysis is based on proton-proton collision data at a center-of-mass energy of 13 TeV collected by the CMS experiment at the LHC in 2016, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The search is performed in five Higgs boson decay channels: h $\to \mathrm{b\bar{b}}$, $\gamma\gamma$, $\tau^{+}\tau^{-}$, W$^{+}$W$^{-}$, and ZZ. The results from the individual channels are combined to maximize the sensitivity of the analysis. No significant excess over the expected standard model background is observed in any of the five channels or in their combination. Limits are set on DM production in the context of two simplified models. The results are also interpreted in terms of a spin-independent DM-nucleon scattering cross section and compared to those from direct-detection DM experiments. This is the first search for DM particles produced in association with a Higgs boson decaying to a pair of W or Z bosons, and the first statistical combination based on five Higgs boson decay channels.


Introduction
A host of astrophysical and cosmological observations confirm [1][2][3][4] that dark matter (DM) exists and makes up 26.4% of the total energy density of the universe [5]. However, all of the existing evidence for DM is based only on its gravitational interaction. Whether DM interacts with standard model (SM) particles in any other way remains an open question. There are a number of beyond-the-SM theories suggesting a particle nature of DM [6]. Several types of particle candidates for DM are proposed in these models, all compatible with the observed relic density of DM in the universe [7]. A favored hypothesis is that the bulk of DM is in the form of stable, electrically neutral, weakly interacting massive particles (WIMPs) [8], with masses in a range between a few GeV and a few TeV, thus opening the possibility of DM production at high-energy colliders [9].
Traditionally, searches for DM at colliders involve a pair of WIMPs that recoil against a visible SM particle or a set of SM particles. Because of the lack of electric charge and the small interaction cross section, WIMPs do not leave a directly detectable signal, but in a hadron collider experiment their presence can be inferred via an imbalance in the total momentum in the plane transverse to the colliding beams ( p miss T ), as reconstructed in the detector. This scenario gives rise to a potential signature where a set of SM particles, X, are produced recoiling against the DM particles, represented by the p miss T (the "mono-X" signature). Recent searches at the CERN LHC considered X to be a hadronic jet [10,11], heavy-flavor quarks (bottom and top) [12, 13], a photon [14, 15], or a W or Z boson [11,[16][17][18].
The discovery of an SM-like Higgs boson [19][20][21] extended the possibility of probing DM at colliders, complementing other mono-X searches. In this paper we designate the state observed at 125 GeV by the symbol h, since in the context of the theoretical models considered below, it does not correspond to the SM Higgs boson. Here, we present a search for the pair production of DM particles in association with a Higgs boson resulting in the final state h + p miss T [22,23], referred to as the "mono-Higgs". While in a typical mono-X search, the X particle is emitted as initial-state radiation, this process is strongly suppressed for the case of the Higgs boson because of the smallness of both the Higgs boson Yukawa couplings to light quarks and its loop-suppressed coupling to gluons. Thus, the mono-Higgs production can be either a result of final-state radiation of DM particles, or of a beyond-the-SM interaction of DM particles with the Higgs boson, typically via a mediator particle. A number of searches have been carried out by the ATLAS and CMS Collaborations looking for the mono-Higgs signature in several Higgs boson decay channels, at center-of-mass energies of 8 and 13 TeV [24][25][26][27][28][29][30][31][32]. So far, none of these searches has observed a significant excess of events over the SM expectations.
In this paper, we describe the first search for mono-Higgs production in the W + W − and ZZ Higgs boson decay channels, as well as the combination of these searches with the previously published results in the bb [30,31], γγ [32], and τ + τ − [32] channels. (Hereafter, for simplicity we refer to bb, τ + τ − and W + W − as bb, ττ and WW, respectively.) All the analyses are based on a data sample of proton-proton (pp) collisions at √ s = 13 TeV collected in 2016 and corresponding to an integrated luminosity of 35.9 fb −1 .
Two simplified models of DM production recommended by the ATLAS-CMS Dark Matter Forum [33] are investigated. Figure 1 shows representative tree-level Feynman diagrams corresponding to these two models. The diagram on the left describes a type-II two Higgs doublet model (2HDM) [34,35] further extended by a U(1) Z group and referred to as the Z -2HDM [36]. In this model, the Z boson is produced via a quark-antiquark interaction and then decays into a Higgs boson and a pseudoscalar mediator A, which in turn can decay to a pair of Dirac fermion DM particles χ. The diagram on the right shows the production mechanism in the baryonic Z model [22], where Z is a vector boson corresponding to a new baryon number U(1) B symmetry. The Z boson acts as a DM mediator and can radiate a Higgs boson before decaying to a pair of DM particles. A baryonic Higgs boson h b is introduced to spontaneously break the new symmetry and to generate the Z boson mass via a coupling that is dependent on the h b vacuum expectation value. The Z boson couplings to quarks and the DM particles are proportional to the U(1) B gauge couplings. A mixing between the h b and h states allows the Z boson to radiate h, resulting in a mono-Higgs signature. In the Z -2HDM, the predicted DM production cross section depends on number of parameters. However, if the mediator A is produced on-shell, the kinematic distributions of the final-state particles depend only on the Z and A boson masses, m Z and m A . In this paper, a scan in m Z between 450 and 4000 GeV and in m A between 300 and 1000 GeV is performed. The values of m A below 300 GeV have been already excluded by the existing constraints on flavor changing neutral currents in the b → sγ transitions [34], and hence are not considered in the analysis. The masses of the 2HDM heavy Higgs boson and the charged Higgs boson are both fixed to the m A mass. The ratio of the vacuum expectation values of the two Higgs doublets, tan β, is varied from 0.4 to 10. The DM particle mass is fixed to 100 GeV, the A-DM coupling strength g χ is fixed to 1, and the Z coupling strength to quarks g Z is fixed to 0.8. The branching fraction of the decay of A to DM particles B(A → χχ) decreases as the mass of the DM candidate (m χ ) increases, for the range of m A considered in this analysis. However, since the relative decrease in B(A → χχ) is less than 7% as m χ increases from 1 to 100 GeV, the results shown in this paper for m χ = 100 GeV are also applicable to lighter DM particles.
The results are expressed in terms of the product of the signal production cross section and branching fraction B(A → χχ), where B(A → χχ) is ≈100% for m A = 300 GeV and decreases for m A greater than twice the mass of the top quark, where the competing decay A → tt becomes kinematically accessible. The contribution to the mono-Higgs signal from another process possible in the model, Z → Z(→ νν ) + h, is not considered in this analysis. Further details on the choice of the model parameters are given in Refs. [27,37]. We note that for the chosen set of parameters, the values of m Z within our sensitivity reach have been recently excluded by the ATLAS and CMS searches for dijet resonances at √ s = 13 TeV [38][39][40][41]. Nevertheless, we keep this benchmark, specifically developed for the LHC Run-2 searches [33], to allow a direct comparison with the results of other mono-Higgs searches. Given that the kinematic distributions of the final states depend only very weakly on the value of the g Z coupling, our results can be reinterpreted for lower g Z values, where the interplay between the mono-Higgs and the dijet analysis sensitivities changes.
For the baryonic Z model, m Z between 100 and 2500 GeV and m χ between 1 and 700 GeV are used for this study. The Z -DM coupling is fixed to g χ = 1 and the Z -quark coupling is fixed to g q = 0.25. The mixing angle between the baryonic Higgs boson and the SM-like Higgs boson is set to sin θ = 0.3, and the coupling between the Z boson and h is assumed to be proportional to m Z . The branching fractions of the Higgs boson decays are altered for m Z m h /2, because the decay h → Z Z ( * ) becomes kinematically accessible. Therefore the region m Z < 100 GeV, for which the modification of the h branching fractions is sizable, is not considered in the analysis. For both benchmark models, h is assumed to have a mass of 125 GeV. A considerable amount of p miss T is expected, as shown in Fig. 2. The reason that the p miss T spectrum is harder for the Z -2HDM is that the DM particles are produced via a resonant mechanism in this case, whereas for the baryonic Z model they are not. The difference in shape becomes more marked as m Z increases. In Fig. 2 (right) it can be seen that the shape of the p miss T distribution is almost independent of m χ in the baryonic Z model, and depends most strongly on m Z . Although the signal sensitivity in the h → bb channel is higher than in the other final states considered (γγ, ττ, WW, and ZZ) because of the channel's large branching fraction and manageable background in the large-p miss T region, the statistical combination of all five decay modes is performed to improve the overall sensitivity. The h → γγ and h → ZZ channels exhibit better resolution in the reconstructed Higgs boson invariant mass, while the h → ττ, h → WW, and h → ZZ channels benefit from lower SM backgrounds, which results in a higher sensitivity for signals with a soft p miss T spectrum.
In the h → bb channel analysis, the h is reconstructed from two overlapping b jets. Thus different approaches are used for the two models, because of the difference in the average Lorentz boost of the Higgs boson, which is higher in the Z -2HDM than in the baryonic Z model. The Higgs boson is reconstructed using a jet clustering algorithm with a distance parameter of 0.8 for the Z -2HDM and 1.5 for the baryonic Z model. For the baryonic Z model, a simultaneous fit of the distribution of the recoil variable in the signal region (SR) and the control regions (CRs) is performed to extract the signal. For the Z -2HDM, a parametric fit of the Z boson transverse mass is used to estimate the major backgrounds and to extract the signal.
The search in the h → γγ channel [32] uses a fit to the diphoton invariant mass distribution to extract the signal. This analysis is performed in two categories distinguished by the p miss T value, high ( >130 GeV) and low , in order to be sensitive to a large variety of possible signals.
The search in the h → ττ channel [32] is based on the combination of the events for the three τ lepton decay modes with the highest branching fractions: τ h τ h , µτ h , and eτ h , where τ h denotes a hadronically decaying τ lepton. After requiring a p miss T (>105 GeV) in order to suppress the background sufficiently, the signal is extracted by performing a simultaneous fit in the SR and in the CRs to the transverse mass of the Higgs boson reconstructed from the two τ leptons. In the h → WW channel search, the fully leptonic decays of the two W bosons are considered, requiring one lepton to be an electron and the other to be a muon, in order to reduce the contamination from the Z → e + e − and Z → µ + µ − backgrounds. The h → ZZ search is performed in the fully leptonic decay channel of the Z boson pair: h → ZZ → 4 . The analysis strategy follows closely the measurement of the Higgs boson properties in the same channel [42].
The paper is organized as follows. After a brief introduction of the CMS detector in Section 2, the data and simulated event samples are described in Section 3. The event reconstruction and the analysis strategy for each Higgs boson decay mode used in the statistical combination are detailed in Sections 4 and 5, respectively. The combination procedure and the main systematic uncertainties are described in Sections 6 and 7, respectively. The results are presented in Section 8, and the paper is summarized in Section 9.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters, made of steel and quartz fibres, extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system [43]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz in a time of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [44].
The pp collision data were collected at √ s = 13 TeV in 2016. The time spacing between adjacent bunches of 25 ns leads to an average number of pp interactions per bunch crossing of 23 assuming the pp inelastic cross section of 69.2 mb [45]. The integrated luminosity of the data sample used in all the analyses described in this paper corresponds to 35.9 fb −1 , after imposing data quality requirements.

Signal and background simulation
Signal samples for the five Higgs boson decay modes are generated at leading order (LO) in perturbative quantum chromodynamics (QCD) using the MADGRAPH5 aMC@NLO v2.3.0 generator [46,47], for both the Z -2HDM and baryonic Z model [33]. The Higgs boson is treated as a stable particle during the generation, and its decays are described subsequently using PYTHIA 8.212 [48].
A detailed description of the simulated samples used for the h → bb, h → γγ, and h → ττ analyses can be found in Refs. [30][31][32]. The production of a Higgs boson in association with a Z boson decaying to a pair of neutrinos is an irreducible background for all the final states considered. Other Higgs boson backgrounds originating from gluon-gluon fusion (ggF) and vector boson fusion (VBF) production modes are small. These backgrounds are simulated at next-to-LO (NLO) in QCD with POWHEG v2 [49][50][51].
The main nonresonant backgrounds in the h → WW analysis are from the continuum WW, single top quark, and top quark pair production. The continuum WW production is simulated in different ways: POWHEG [52] is used to generate qq → WW events at NLO precision, whereas gg → WW events are generated at LO using MCFM v7.0 [53][54][55]. The simulated qq → WW events are reweighted to reproduce the p WW T distribution from the p T -resummed calculation at next-to-NLO (NNLO) plus next-to-next-to-leading logarithmic precision [56,57]. The LO gg → WW cross section, obtained directly from MCFM, is further corrected to NNLO precision via a K factor of 1.4 [58]. Single top quark, tt, WZ, and Wγ * backgrounds are generated at NLO with POWHEG. Drell-Yan (DY) production of Z/γ * is generated at NLO using MADGRAPH5 aMC@NLO, and the p T spectrum of the dilepton pairs is reweighted to match the distribution observed in dimuon events in data. Other multiboson processes, such as Wγ, ZZ, and VVV (V = W or Z), are generated at NLO with MADGRAPH5 aMC@NLO. All samples are normalized to the latest available theoretical cross sections, NLO or higher [53,54,59].
In the h → ZZ analysis, the SM production mechanism constitutes a major background because this has the same experimental signature and satisfies the low p miss T threshold used in the analysis. It is simulated with POWHEG [49,50,60] in four main production modes: ggF, including quark mass effects [61]; VBF [62]; associated production with a top quark pair (tth) [63]; and associated production with a vector boson (Wh, Zh), using the MINLO HVJ [64] extension of POWHEG. In all cases, the Higgs boson is forced to decay via the h → ZZ → 4 ( = e, µ, or τ) channel. The description of the decay of the Higgs boson to four leptons is obtained using the JHUGEN 7.0.2 generator [65,66]. In the case of Zh and tth production, the Higgs boson is allowed to decay as h → ZZ → 2 + X, such that four-lepton events where two leptons originate from the decay of the associated Z boson or top quarks are also taken into account in the simulation. The cross sections for the processes involving SM Higgs boson production are taken from Ref. [67].
All processes are generated using the NNPDF3.0 [68] parton distribution functions (PDFs), with the precision matching the parton-level generator precision. The PYTHIA generator with the underlying event tune CUETP8M1 [69] is used to describe parton showering and fragmentation. The detector response is simulated using a detailed description of the CMS apparatus, based on the GEANT4 package [70]. Additional simulated pp minimum bias interactions in the same or adjacent bunch crossings (pileup) are added to the hard scattering event, with the multiplicity distribution adjusted to match that observed in data.

Event reconstruction
The particle-flow (PF) algorithm [71] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of photons is obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track [72]. The energy of muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
Electron candidates are required to have |η| < 2.5. Additional requirements are applied to reject electrons originating from photon conversions in the tracker material or jets misreconstructed as electrons. Electron identification criteria rely on observables sensitive to the bremsstrahlung along the electron trajectory and on the geometrical and momentum-energy matching between the electron track and the associated energy cluster in the ECAL, as well as on the ECAL shower shape observables and association with the primary vertex.
Muon candidates are reconstructed within |η| < 2.4 by combining information from the silicon tracker and the muon system. Identification criteria based on the number of measurements in the tracker and in the muon system, the fit quality of the muon track, and its consistency with its origin from the primary vertex are imposed on the muon candidates to reduce the misidentification rate.
For each event, hadronic jets are clustered from PF candidates using the infrared-and collinearsafe anti-k T algorithm [73,74], with a distance parameter of 0.4 (AK4 jets) or 0.8 (AK8 jets). Jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, on average, within 5 to 10% of the true momentum over the entire p T spectrum and detector acceptance. Pileup interactions can result in additional spurious contributions to the jet momentum measurement from tracks and calorimetric energy depositions. To mitigate this effect, tracks identified to be originating from pileup vertices are discarded and a correction based on the jet area [75] is applied to account for the neutral pileup particle contributions. Jet energy corrections are derived from simulation to bring the measured response of jets to that of particle-level jets on average. In situ measurements of the momentum balance in dijet, photon+jet, Z+jet, and multijet events are used to account for any residual differences in the jet energy scale (JES) between data and simulation [76]. The jet energy resolution (JER) amounts typically to 15% at p T = 10 GeV, 8% at 100 GeV, and 4% at 1 TeV. Additional selection criteria are applied to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures [77].
At large Lorentz boosts, the two b quarks from the Higgs boson decay may produce jets that overlap and make their individual reconstruction difficult. In this case, either the AK8 jets or larger-area jets clustered from PF candidates using the Cambridge-Aachen algorithm [78,79] with a distance parameter of 1.5 (CA15 jets) are used. To reduce the impact of particles arising from pileup interactions when reconstructing AK8 or CA15 jets, the four-vector of each PF candidate matched to the jet is scaled with a weight calculated with the pileup-per-particle identification algorithm [80] prior to the clustering. The CA15 jets are also required to be central (|η| < 2.4). The "soft-drop" jet grooming algorithm [81] is applied to remove soft, large-angle radiation from the jets. The mass of a groomed AK8 or CA15 jet are referred to as m SD .
To identify jets originating from b quark fragmentation (b jets), two b tagging algorithms are used. The combined secondary vertex (CSVv2) [82] and the combined multivariate analysis (cMVAv2) algorithms [82] are used to identify AK4 jets originating from b quarks by their characteristic displaced vertices. For the AK8 jets, subjets inside the jet are required to be tagged as b jets using the CSVv2 algorithm. A likelihood for the CA15 jet to contain two b quarks is derived by combining the information from the primary and secondary vertices and tracks in a multivariate discriminant optimized to distinguish CA15 jets originating from the h → bb decay from those produced by energetic light-flavor quarks or gluons [31].
Hadronically decaying τ leptons are reconstructed from jets using the hadrons-plus-strips algorithm [83]. This algorithm uses combinations of reconstructed charged hadrons and energy deposits in the ECAL to identify the three most common hadronic τ lepton decay modes: 1prong, 1-prong+π 0 (s), and 3-prong. The τ h candidates are further required to satisfy the isolation criteria with an efficiency of 65 (50)% and a misidentification probability of 0.
The p miss T is reconstructed as the negative vectorial sum of all PF particle candidate momenta projected on the plane transverse to the beams. Since the presence of pileup induces a degradation of the p miss T measurement (p miss T resolution varies almost linearly from 15 to 30% as the number of vertices increases from 5 to 30 [84]), affecting mostly backgrounds with no genuine p miss T , an alternative definition of p miss T that is constructed only using the charged PF candidates ("tracker p miss T ") is used in the h → WW analysis. In the rest of the paper, p miss T corresponds to the PF p miss T , unless specified otherwise.

Analysis strategy
In this section we briefly discuss the analysis strategies in the previously published [30-32] h → bb, h → γγ, and h → ττ, channels, and provide full descriptions of the new analyses in the h → WW and h → ZZ decay channels. The summary of all the decay channels contributing to the combination is presented in Table 1.

T channel
The events used in this final state are selected using a triggers that require large amount (> 90 or > 120 GeV) of p miss T , or H miss T defined as the magnitude of the vectorial sum of the transverse momenta of all jets with p T > 20 GeV in an event. The trigger selection is 96 (100)% efficient for events that subsequently have p miss T > 200 (350) GeV in the off-line reconstruction. As can be seen in Fig. 2, the Lorentz boosts of the Higgs boson are different for the Z -2HDM and baryonic Z model. The events with large boost in the Z -2HDM are reconstructed using a large-radius AK8 jet with p T > 200 GeV and |η| < 2.4. In addition, the h → bb topology is selected by requiring at least one subjet of the AK8 jet to be b tagged. The analysis considers separately two categories, distinguished by the number of b tagged subjets in the event, one or two, the latter being the high-purity category with higher sensitivity. For events with lower boost in the baryonic Z model, Higgs boson candidates are reconstructed using CA15 jets.
To select the h → bb candidates using the AK8 jet, one or both subjets are required to pass the loose b tagging criteria, which has an efficiency of 85%, and a misidentification rate of about 10% for jets originating from light-flavor quarks or gluons. In the case of the CA15 jets, a multivariate double b tagging algorithm [82] is used to discriminate the signal from the background of light-flavor jets [31], with an efficiency of 50% and a misidentification rate of Table 1: Summary of the individual channels entering the combination. Analyses are categorized based on the model, p miss T selection, and subsequent decay products listed here. The categorization is the same for both the Z -2HDM and the Baryonic Z model for all decay channels except, as indicated, h → bb. A dash ("-") in the last column implies that the analysis is presented in this paper.

Decay channel Final state or category Reference
10%. The AK8 (CA15) analysis requires the Higgs boson candidate mass to be in the 105-135 (100-150) GeV range to reduce nonresonant backgrounds. The difference in the two mass window requirements is primarily driven by the differences in the performance of the two algorithms and in the jet mass resolutions. For both analyses, the mass window was chosen to maximize the signal sensitivity. In order to further reduce the background contributions from W + jets and tt production, events with an electron, muon, photon (p T > 10 GeV), or τ h (p T > 18 GeV) candidates passing loose identification and isolation criteria are vetoed. Furthermore, in the AK8 analysis, the number of additional b tagged AK4 jets with p T > 20 GeV is required to be zero, while in the CA15 analysis, the number of AK4 jets with p T > 30 GeV, well-separated from the CA15 jet in the event, is required to be at most one. The sensitivity of the analyses is further enhanced by using jet substructure variables. The full details of the event selection for the AK8 and CA15 jet analyses can be found in Refs.

T channel
Signal candidate events in the h → γγ analysis are selected using a diphoton trigger with asymmetric p T thresholds of 30 and 18 GeV on the leading and subleading photons, respectively, and loose identification and isolation requirements imposed on both photon candidates. The diphoton invariant mass is further required to exceed 90 GeV.
Slightly higher thresholds of 30 (20) GeV on the leading (subleading) photon p T and of 95 GeV on the diphoton mass are used offline. The photon candidates are required to pass the isolation criteria if the spatial distance in η-φ plane (∆R = √ (∆η) 2 + (∆φ) 2 ) between the two photons exceeds 0.3. The isolation selection is not used for photons that are coming from the decay of a highly Lorentz-boosted Higgs boson, as the two photons are likely to be found in the isolation cone of one another. The analysis is performed in two categories distinguished by the value of p miss T : high-p miss T (>130 GeV) and low-p miss T (50-130 GeV).
The multijet background, with a large p miss T in an event originating from the mismeasurement of the energy of one or more jets, is reduced by allowing at most two jets with p T > 30 GeV. To suppress the contribution from the multijet background, the azimuthal separation between 5.3 The h(→ τ τ) + p miss T channel 9 the direction of any jet with p T > 50 GeV and p miss T is required to exceed 0.5 radians. Finally, to select signal-like events with the DM particles recoiling against the Higgs boson, the azimuthal separation between p miss T and the direction of the Higgs boson candidate reconstructed from the diphoton system is required to exceed 2.1 radians. More details of the event selection can be found in Ref. [32].

T channel
In the h → ττ analysis, the three final states with the highest branching fractions are analyzed: τ h τ h , µτ h , and eτ h . The events are selected online with a trigger requiring the presence of two isolated τ h candidates in the τ h τ h final state, and a single-muon (single-electron) trigger in the µτ h (eτ h ) final state. Electron, muon, and τ h candidates passing the identification and isolation criteria are combined to reconstruct a Higgs boson candidate in these three final states. The signal events are then selected with the requirements: p miss T > 105 GeV and visible p T of the ττ system > 65 GeV. To ensure that the ττ system originates from the Higgs boson, the visible mass of the ττ system is required to be less than 125 GeV. In order to reduce the contribution from multilepton and tt backgrounds, the events are vetoed if an additional electron, muon, or a b tagged jet is present. More details of the event selection can be found in Ref. [32].

T channel
The search in the h → WW decay channel is performed in the fully leptonic, opposite-sign, different-flavor (eµ) final state, which has relatively low backgrounds. The presence of the neutrinos and the DM particles escaping detection results in large p miss T in signal events. The selected eµ + p miss T events include a contribution from the h → WW → ττν τ ν τ process with both τ leptons decaying leptonically. Several background processes can lead to the same final state, dominated by tt and WW production.
Online, events are selected using a suite of single-and double-lepton triggers. In the offline selection, the leading (subleading) lepton is required to have p T > 25 (20) GeV. Electron and muon candidates are required to be well-identified and isolated to reject the background from leptons inside jets. Backgrounds from low-mass resonances are reduced by requiring the dilepton invariant mass (m ) to exceed 12 GeV, while backgrounds with three leptons in the final state are reduced by vetoing events with an additional well-identified lepton with p T > 10 GeV. The p miss T in the event is required to exceed 20 GeV in order to reduce the contribution from instrumental backgrounds and Z/γ * → τ + τ − decays. To suppress the latter background, the p T of the dilepton system is required to be greater than 30 GeV and the transverse mass of the dilepton and p miss T system, m h T , is required to be greater than 40 GeV. In order to reduce the Z/γ * → e + e − , µ + µ − or τ + τ − background with p miss T originating either from τ lepton decays or from mismeasurement of the energies of e, µ or additional jets, a variable p miss T,proj [85] is introduced. This is defined as the projection of p miss T in the plane transverse to the direction of the nearest lepton, unless this lepton is situated in the opposite hemisphere to p miss T , in which case p miss T,proj is taken to be p miss T itself. A selection using this variable efficiently rejects Z/γ * → background events, in which the p miss T is preferentially aligned with leptons. Since the p miss T resolution is degraded by pileup, a quantity p miss T,mp is defined as the smaller of the two p miss T,proj values: the one based on all the PF candidates in the event, and the one based only on the reconstructed tracks originating from the primary vertex. A requirement p miss T,mp > 20 GeV is effective in suppressing the targeted background. The above requirements define the event preselection.
The expected signal significance is enhanced by introducing two additional selections: m < 76 GeV and the distance in η-φ space between the two leptons ∆R < 2.5, as illustrated in Fig. 3. The first requirement exploits the fact that the invariant mass of the leptons coming from the h → WW decay tends to be low because of the presence of the two neutrinos in the decay chain and of the scalar nature of the Higgs boson. The second requirement utilizes the fact that the Higgs boson in signal events recoils against the DM particles and is highly boosted.

Background estimation
Since full kinematic reconstruction of the Higgs boson mass and p T is impossible in this decay channel because of the presence of undetected neutrinos and DM particles, to maximize the sensitivity of the search, a boosted decision tree (BDT) multivariate classifier has been trained for each of the two signal models. The BDT exploits the following input variables: Here, m defines the transverse mass of p miss T and the leading (subleading) lepton in the event, and ∆φ is the azimuthal angle between the directions of the two lepton momenta.
For both benchmark models, the BDT training considers processes with two prompt leptons and genuine p miss T (WW, tt, tW, and h → WW production) as the backgrounds. For the Z - The main background processes arise from top quark (tt and single top quark production, mainly tW), nonresonant WW events, and nonprompt leptons. The contribution of nonpromptlepton background in the SR is determined entirely from data, while the contributions of the top quark, WW, and Z/γ * → τ + τ − background are estimated using simulated samples. The normalizations of simulated backgrounds are obtained using dedicated CRs that are included in the maximum-likelihood fit used to extract the signal, together with the SR. Smaller backgrounds, WZ and Wγ * , are estimated using simulation after applying a normalization factor estimated in the respective CRs. The WZ CR is defined by requiring the presence of two opposite-sign, same-flavor leptons, compatible with the decay of a Z boson and one additional lepton of a different flavor, consistent with originating from a W boson decay. In the Wγ * CR, the two leptons produced by the decay of the virtual photon are required to have p T > 8 GeV and be isolated. Since the two leptons may be close to each other, the isolation is computed without taking into account the contribution of lepton tracks falling in the isolation cone. An additional lepton consistent with originating from the W decay is required. The WZ and Wγ * CRs are not used in the maximum-likelihood fit; instead, the normalization scale factors are extracted and directly applied to the corresponding simulated samples. The remaining backgrounds from diboson and triboson production are estimated directly from simulation.
The gg → W + W − and qq → W + W − backgrounds are estimated from simulation normalized as discussed in Section 3. The main feature of these processes is that, as the two W bosons do not originate in a decay of the Higgs boson, their invariant mass does not peak at the Higgs boson mass. For this reason, events in the corresponding CR are required to have a large dilepton invariant mass, achieved by inverting the SR m < 76 GeV requirement.
The estimation of the top quark background is performed in two steps. First, a top quark enriched CR is defined to measure a scale factor quantifying the difference in the b tagging efficiencies and mistag rates in data and simulation. This CR is obtained from the SR selection by inverting the b tagged jet veto. In second step, the scale factor is applied to the corresponding simulated samples with a weight per event that depends on the number, flavor, and kinematic distributions of jets.
The W + jets production contributes as a background in the h → WW analysis when a jet is misidentified as a lepton. A CR is defined to contain events with one isolated lepton and another lepton candidate that fails the nominal isolation criteria, but passes a looser selection. The probability for a jet satisfying this looser selection to pass the nominal one is estimated from data in an independent sample dominated by nonprompt leptons from multijet production. This probability is parameterized as a function of the p T and η of the lepton and applied to the events in the CR. In order to estimate the nonprompt lepton contamination in the SR, a validation region enriched in nonprompt leptons is defined with the same requirement as the SR, but requiring same-sign eµ pairs. The maximum discrepancy between data and prediction in the validation region, amounting to ≈30%, is taken as the uncertainty in the W + jets background prediction.
The Z/γ * → τ + τ − background is estimated from simulation, after reweighting the Z boson p T spectrum to match the distribution measured in data. The normalization of the simulated sample is estimated from data using events in the m h T < 40 GeV region. A normalization factor is then extracted from this region and applied to the SR.
The main difference between the present analysis and the measurement of the SM Higgs boson properties in the same channel [85] is in the signal extraction method. The latter analysis uses a multidimensional fit to the m h T , m , and p 2 T distributions, whereas a fit to the BDT discriminant distribution is used in the present analysis. The signal event topology is defined by the presence of four charged leptons (4e, 4µ, or 2e2µ) and significant p miss T produced by the undetected DM particles. The events are selected online with triggers requiring the presence of two isolated leptons (ee, µµ, or eµ), with asymmetric p T thresholds of 23 (17) GeV on the leading and 12 (8) GeV on the subleading electron (muon). Dilepton triggers account for most of the signal efficiency in all three final states. In order to maximize the signal acceptance, trilepton triggers with lower p T thresholds and no isolation requirements are added, as well as single-electron and single-muon triggers with isolated lepton p T thresholds of 27 and 22 GeV, respectively [42].

The h(→ ZZ
The reconstruction and selection of the Higgs boson candidates proceeds first by selecting two Z boson candidates, defined as pairs of opposite-sign, same-flavor leptons (e + e − , µ + µ − ) passing the selection criteria and satisfying 12 < m (γ) < 120 GeV, where the Z boson candidate mass m (γ) includes the contribution of photons identified as coming from final-state radiation [42]. The ZZ candidates are then defined as pairs of Z boson candidates not sharing any of the leptons. The Z candidate with the reconstructed mass closest to the nominal Z boson mass [86] is denoted as Z 1 , and the other one is denoted as Z 2 . All the leptons used to select the Z 1 and Z 2 candidates must be separated by ∆R( i , j ) > 0.02.
The leading (subleading) of the four leptons must have p T > 20 (10) GeV, and the Z 1 candidate must have a reconstructed mass m Z 1 above 40 GeV. In the 4e and 4µ channels, if an alternative Z i Z j candidate based on the same four leptons is found, the event is discarded if m Z i is closer to the nominal Z boson mass than m Z 1 . This requirement rejects events with an on-shell Z boson produced in association with a low-mass dilepton resonance. In order to suppress the contribution of QCD production of low-mass dilepton resonances, all four opposite-sign pairs that can be built with the four leptons (regardless of the lepton flavor) must satisfy m i j > 4 GeV and the four-lepton invariant mass must satisfy m 4 > 70 GeV. If more than one ZZ candidate passes the selection, the one with the highest value of the scalar p T sum of four leptons is chosen. The above requirements define the event preselection.
The m 4 distribution for selected ZZ candidates exhibits a peak around 125 GeV, as expected for both the SM Higgs boson production and signal. However, because of the much lower cross section, the potential signal is overwhelmed by the background after the SM Higgs boson selection, as shown in Fig. 4 (left). The distribution of p miss T for selected ZZ candidates is shown in Fig. 4 (right). After the preselection, the remaining background comes from the SM Higgs boson (mostly Vh), tt+V, and VV/VVV production. Another background dominated by the Z+jets production ("Z+X") [42] arises from secondary leptons misidentified as prompt because of the decay of heavy-flavor hadrons and light mesons within jets, and, in the case of electrons, from photon conversions or charged hadrons overlapping with photons from π 0 → γγ decays. The nonprompt-lepton background also contains smaller contributions from tt+jets, Zγ+jets, WZ+jets, and WW+jets events, with a jet misidentified as a prompt lepton. These backgrounds do not exhibit peak in the distribution of m 4 , and are reduced by applying a selection on the m 4 around the Higgs boson mass (115 < m 4 < 135 GeV), by rejecting events with more than four leptons, and by requiring the number of b tagged jets in the event to be less than two.

Background estimation
The dominant irreducible backgrounds from the SM Higgs boson and nonresonant ZZ production are determined from simulation, while the Z+X background is determined from data [42]. All other backgrounds are determined from simulation. Background contributions from the SM Higgs boson production in association with a Z boson or a tt pair, followed by the h → WW → 2 2ν decay, have been studied with simulated events and found to be negligible.
The Z+X background is estimated from data by first determining the lepton misidentification probability in a dedicated CR and then using it to derive the background contribution in the SR. The lepton misidentification probability is defined as the probability that a lepton passing a loose selection with relaxed identification or isolation criteria also passes the tight selection criteria. The misidentification probability is measured in a Z+lepton CR where the Z boson candidate (with the mass within 7 GeV of the nominal Z boson mass) is formed from the two selected leptons passing the tight identification criteria, and an additional lepton is required Table 2: Summary of the maximum number of additional objects allowed in an event for each analysis. A dash means that no restriction on the corresponding object is applied in the corresponding analysis.
to pass the loose selection. This sample is dominated by Z+nonprompt-lepton events. The electron and muon misidentification probabilities are measured as functions of the lepton candidate p T , its location in the barrel or endcap region of the ECAL or the muon system, and p miss T in the event, using Z(→ )+e and Z(→ )+µ events, respectively, in the Z+lepton CR. The misidentification probabilities are found to be independent of the charge of the lepton within the uncertainties.
The strategy for applying the lepton misidentification probabilities relies on two additional CRs. The first CR is defined by requiring that the two leptons that do not form the Z 1 candidate, pass only the loose, but not the tight identification criteria. This CR defines the "2 pass + 2 fail" (2P2F) sample and is expected to be populated by events that intrinsically have only two prompt leptons (mostly from DY production, with a small contribution from tt and Zγ events). The second CR is defined by requiring only one of the four leptons to fail the tight identification and isolation criteria and defines the "3 pass + 1 fail" (3P1F) sample, which is expected to be populated by the type of events that populate the 2P2F CR, but with different relative proportions, as well as by WZ+jets events with three prompt leptons.

Statistical combination of the search channels
The analyses in the five channels described above are almost completely statistically independent of each other, allowing these analyses to be combined without accounting for the possibility of events being selected in more than one final state. Whenever an explicit veto ensuring the strict mutual exclusivity of the channels is not placed in a particular analysis, it was checked that there are no overlapping events with the other channels. The summary of the vetoes on additional objects, namely electrons, muons, τ leptons, photons, jets, and b tagged jets, in each analysis is presented in Table 2. These selections not only reduce the major backgrounds, but also ensure the nearly complete mutual exclusivity of the analyses considered for the combination. The overlap in the SR is zero and for the CR it is less than 0.01%, i.e., it is much smaller than the systematic uncertainty in the analysis. For the Z -2HDM, the two parameters that we scan are m Z and m A . All five analyses contribute to the combination in the ranges 800 < m Z < 2500 GeV and 300 < m A < 800 GeV. For m Z < 800 GeV, it is not possible to perform the h → bb analysis efficiently, therefore only four other decay channels are used for the combination. For m Z > 2500 GeV and m A > 800 GeV the signal selection efficiency is significant only for the h → bb decay mode, hence only the h → bb channel contributes in this region.
For the baryonic Z model, the two parameters that we scan are m Z and m χ , and all five analyses are performed in the full phase space considered for the combination. Since the maximum sensitivity for all the analyses is achieved for m χ = 1 GeV, the comparison of individual analyses is shown only for this DM particle mass, to demonstrate the improvement in the sensitivity achieved in the combination of individual channels.

Systematic uncertainties
A number of systematic uncertainties are considered in the combination, broadly divided into two categories: theoretical and experimental. Theoretical uncertainties are considered fully correlated among all five channels. Only the systematic uncertainties attributed to the experimental sources that are correlated between different channels are described for the combined result in section 7.3. The details of all experimental systematic uncertainties in the h → bb analysis using AK8 jets are described in Ref.
[30] and those for the analysis using CA15 jets are described in Ref.
[31]; for the h → γγ and h → ττ channels they are given in Ref.
[32]; and for the h → WW and h → ZZ analyses they are discussed in this section.

T channel
The normalization and the kinematic shapes of the BDT discriminant distributions for the main backgrounds are derived from data CRs, and therefore systematic uncertainties in both the normalization and shapes are considered.
For the nonprompt-lepton background the uncertainty amounts to approximately 30%, and covers the uncertainty in the lepton misidentification rate, the dependence on the CR background composition, and the statistical component because of the finite event count in the CR.
The top quark background CR is included as an additional category in the signal extraction fit. The kinematic shapes of the top quark background are taken from simulation corrected for the b tagging scale factors, with the uncertainties covering the difference between the b tagging efficiency in data and simulation [82]. A similar procedure is applied for the DY background, by defining a CR in low-m T phase space, and to the nonresonant WW background, for which a high-m CR is defined. Experimental uncertainties are estimated by applying scale factors between data and simulation, and/or by smearing of certain kinematic variables in simulation, with the corresponding changes further propagated to all analysis variables. The signal acceptance uncertainty associated with the combination of single-lepton and dilepton triggers is measured to be 2%. The uncertainty in the ratio between the single top quark and top quark pair production cross sections, 8% at 13 TeV [87], has been also included, as it affects the top quark background yield from the maximum-likelihood fit used to extract the signal and dominant backgrounds. The uncertainty in the p T spectrum of the top quark has been applied to all the observables in order to cover the difference between the simulated and observed spectra [88], and is of the order of 1%.
The uncertainty in the Higgs boson branching fraction for the h → WW decay is about 1% [67]. The uncertainty in the NNLO K factor applied to the LO gg → WW cross section estimate is 15% [89]. The p WW T spectrum in the qq → WW sample has been reweighted to match the resummed calculation [56,57]. The associated shape uncertainties related to the missing higherorder corrections are modeled by varying the factorization, renormalization, and resummation scales up and down independently by a factor of 2 from their nominal values [56]. Finally, uncertainties arising from the limited size of the simulated samples are included for each bin of the BDT discriminant distributions, in each category. The main sources of the uncertainties affecting the analysis are listed in Table 3.

T channel
A source of systematic uncertainty in the nonprompt-lepton background estimate potentially arises from the difference in the composition of the SM background processes with nonprompt leptons (Zγ+jets, tt, Zγ+jets) contributing to the CRs where the lepton misidentification rate is measured and applied. This uncertainty can be estimated by measuring the misidentification rates in simulation for the 2P2F and 3P1F CRs. Half of the difference between the misidentification rates obtained from simulation in these two CRs is used as a measure of the systematic uncertainty in the lepton misidentification rate and is further propagated to the uncertainty in the nonprompt-lepton background, and amounts to 43% for the 4e, 36% for the 4µ, and 40% for the 2e2µ final states.
The uncertainty in the full signal selection efficiency is at the level of 1%. The uncertainty in the m 4 resolution from the uncertainty in the per-lepton energy resolution is about 20% [42] and affects the signal and all the backgrounds from Higgs boson production.
In addition, there are two types of systematic uncertainties related to the modeling of p miss T . The first uncertainty is related to the approximately Gaussian core of the resolution function for correctly measured jets and other physics objects and corresponds to the uncertainty in the genuine p miss T . The second uncertainty, attributed to significant mismeasurement of p miss T , is an uncertainty in the "mismeasured" p miss T . The uncertainties from the modeling of genuine p miss T are measured by varying the parameters associated with the corrections applied to p miss T and by propagating those variations to the p miss T calculation, after applying the full analysis selection. Each correction is varied up and down by one standard deviation of the input distribution. The corrections used in this calculation come from JES, JER, muon, electron, photon, and the unclustered energy scales.
The uncertainty in the mismeasured p miss T is obtained from a sample with significant contributions from misidentified leptons and mismeasured jets, obtained by requiring an oppositesign, same-flavor dilepton pair passing the Z 1 candidate selection, and an additional same-sign, same-flavor pair ("OS+SS" sample). This sample is enriched in misidentified leptons that form the same-sign pair and is expected to lead to significant mismeasurement of p miss T , not already covered by the uncertainties in the Gaussian core discussed above. We derive the mismeasured p miss T uncertainty from the comparison of the p miss T shapes in the "OS+SS" sample and in the SR, with a requirement that the m 4 be outside the Higgs boson invariant mass peak (|m 4 − 125 GeV| > 10 GeV). The uncertainty in mismeasured p miss T is applied to the Z+X sample only, since the effect is expected to be negligible when four genuine leptons are produced, as is the case for the signal and for most of the simulated background samples.
An uncertainty of 10% in the K factor used for the gg → ZZ prediction is applied [89]. A systematic uncertainty of 2% in the h → ZZ → 4 branching fraction [67] affects both signal and the SM Higgs boson background yields. Theoretical uncertainties in the tt+V background cross sections are taken from Ref. [90]. A summary of the experimental uncertainties is given in Table 4.

Systematic uncertainties in the combination
The uncertainties associated with the background normalization and fit parameters are assumed to be uncorrelated, whereas those associated with the standard object selection are considered fully correlated and are summarized in Table 5. In all five decay channels, a normalization uncertainty of 2.5% for simulated samples is used to account for the uncertainty in the measurement of the integrated luminosity [91]. Also fully correlated across all channels are the systematic uncertainties related to theoretical calculations of the Higgs boson production Table 5: Systematic uncertainties in the combination of channels, along with the type (rate/shape) of uncertainty affecting signal and background processes, correlated amongst at least two final states. For the rate uncertainties, the percentage of the prior value is quoted, while for shape uncertainties an estimate of the impact of systematic uncertainties on the yield is also listed. A dash ("-") implies that a given uncertainty does not affect the analysis. Whenever an uncertainty is present but kept uncorrelated in a particular channel, this is mentioned explicitly. The effect of the b jet mistag rate uncertainty is very small in the h → bb Z'-2HDM analysis and hence it is added to the effect of the b tagging efficiency uncertainty in quadrature. cross section, PDFs, and renormalization and factorization scale uncertainties estimated using the recommendations of the PDF4LHC [92] and LHC Higgs Cross Section [67] working groups, respectively. These uncertainties range from 0.3 to 9.0%.
Uncertainties from imprecise knowledge of the JES are evaluated by propagating the uncertainties in the JES for individual jets in an event, which depend on the jet p T and η, to all the analysis quantities. The uncertainties in the selection of b tagged AK4 jets are taken into account using the uncertainties in the b tagging efficiency and misidentification rate estimated from the difference between data and simulation [82]. The uncertainty due to the difference in the performance of electron, muon, and τ lepton identification between data and simulation is taken into account for individual decay channels and considered fully correlated in the statistical combination. An uncertainty of 1-3% in the electron energy scale and an uncertainty of 0.4-1.0% in the muon energy scale are considered to be correlated in the combination.

Results
The event selection described in Section 5 has been used to discriminate the mono-Higgs signal from backgrounds in each channel. The observed yields in data and the expected event yields for the signal and background processes in the h → bb, h → γγ, and h → ττ channels can be found in Refs. [30][31][32]. The corresponding yields for the h → WW and h → ZZ analyses are discussed in Section 8.1. Tables 6, 7 and figures 5, 6 show one signal mass hypothesis for each model, normalized to the respective cross section. For the Z -2HDM, the signal is normalized to the cross section calculated for mass values of Z and A bosons of 1200 and 300 GeV, respectively, and for g Z = 0.8, tan β = 1. For the baryonic Z model, the signal is normalized to the cross section corresponding to the Z and m χ masses of 500 and 1000 GeV, respectively, and for g χ = 1, g q = 0.25.   The expected background yields and the observed number of event in data, along with the expected yields for two signal benchmarks in the h → WW and h → ZZ channels, are summarized in Tables 6 and 7, respectively. Figure 5 shows the BDT discriminant distribution for the expected backgrounds and observed events in data for the h → WW analysis. Benchmark signal contributions in the Z -2HDM (left) and baryonic Z (right) model are also shown, scaled by the factors of 500 and 100, respectively, for better visibility. Figure 6 shows the p miss T distribution of the expected backgrounds and observed events in data for the h → ZZ analysis. Benchmark signal contributions are also shown. For both analyses, the total uncertainty, given by a quadratic sum of the statistical and systematic components, is shown. The bottom panels show the ratios of data to the total background prediction with their total uncertainties. The potential signal is extracted from the fit to the BDT discriminant (p miss T ) spectrum with a signal-plus-background hypothesis for the h → WW (h → ZZ) channel. The profile likelihood ratio is used as a test statistic, in an asymptotic approximation [93]. Data agree well with the expected background and no signal is observed in either channel. Limits on the model parameters at 95% confidence level (CL) are set using the modified frequentist CL s criterion [94][95][96] with all the nuisance parameters profiled.
The observed and expected upper limits on the DM candidate production cross section are shown in Fig. 7 for the h → WW (upper) and h → ZZ (lower) channels for the Z -2HDM with m A = 300 GeV (left) and for the baryonic Z model with the value of m χ fixed at 1 GeV (right). All other model parameters are fixed to the values described in Section 1. The upper limits for the h → ZZ analysis already include the statistical combination of all three final states used. The h → WW analysis excluded the region from 780 to 830 GeV for m A = 300 GeV in the Z -2HDM.

Results of the statistical combination
The observed and expected upper limits at 95% CL on the DM production cross section normalized to the predicted cross section, as a function of m Z , from the combination of all five channels are shown in Fig. 8 for the Z -2HDM with m A = 300 GeV (left) and for the baryonic For the Z -2HDM, the combination is dominated by the h → bb analysis for m Z > 800 GeV. However, the h → bb analysis has no sensitivity for m Z values below 800 GeV, and a combination of the h → γγ and h → ττ channels plays a significant role in this region of the model parameter space. The range of m Z excluded at 95% CL spans from 500 to 3200 GeV for m A = 300 GeV.
For the baryonic Z model, the combination results are also dominated by the h → bb channel, but the h → γγ and h → ττ channels also provide a nonnegligible contribution in constraining the model parameters. The range of m Z excluded at 95% CL spans from 100 to 1600 GeV for m χ = 1 GeV. Figure 9 shows the observed and expected 95% CL exclusion contours on σ/σ th in the m Z -m A and m Z -m χ planes for the Z -2HDM (left) and baryonic Z (right) model, respectively.
The results for the Z -2HDM are also interpreted in the m Z -tan β plane for three different m A values: 300, 400, and 600 GeV. Since the shape of the p miss T distribution does not change with tan β, and affects only the product of the Z production cross section and branching fraction to the mono-h channel, the limit shown in Fig. 9 (left) can be simply rescaled for different values of tan β, from 0.5 to 10. These limits, in the m Z -tan β plane, are shown in Fig. 10 which also couples to the SM quarks. A point in the parameter space of this model is determined by four variables: the DM particle mass m χ , the mediator mass m med , the mediator-DM coupling g χ , and the universal mediator-quark coupling g q . The couplings for the present analysis are fixed to g χ = 1.0 and g q = 0.25, following the recommendation of Ref.
[37]. The results are interpreted in terms of 90% CL limits on the spin-independent (SI) cross section σ SI for the DM-nucleon scattering. The value of σ SI for a given set of parameters in the s-channel simplified DM model is given by [37]: where µ nDM is the reduced mass of the DM-nucleon system and f (g q ) is the mediator-nucleon coupling, which depends on g q . The resulting σ SI limits, as a function of m χ are shown in    Figure 9: The upper limits at 95% CL on the observed and expected σ/σ th in the m Z -m A and m Z -m χ planes for the Z -2HDM (left) and baryonic Z model (right), respectively. The region enclosed by the contours is excluded using the combination of the five decay channels of the Higgs boson for the following benchmark scenarios: g Z = 0.8, g χ = 1, tan β = 1, m χ = 100 GeV, and m A = m H = m H ± for the Z -2HDM, and g χ = 1, g q = 0.25 for the baryonic Z model.    Figure 11: The upper limits at 90% CL on the DM-nucleon spin-independent scattering cross section σ SI , as a function of m χ . Results obtained in this analysis are compared with those from the CMS dijet analyses [39, 41] and from several direct-detection experiments: CRESST-II [97], CDMSLite [98], PandaX-II [99], LUX [100], XENON-1T [101], and CDEX-10 [102]. Fig. 11. Results obtained in this analysis are compared with those from the CMS dijet analyses 1 [39, 41] and from several direct-detection experiments. For the chosen set of parameters, the cross section limit from the present analysis is more stringent than the direct-detection limits for m χ between 1 and 5 GeV.

Summary
A search for dark matter particles produced in association with a Higgs boson has been presented, using a sample of proton-proton collision data at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 . Results from five decay channels of the Higgs boson, h → bb, h → γγ, h → τ + τ − , h → W + W − , and h → ZZ, are described, along with their statistical combination. No significant deviation from the standard model prediction is observed in any of the channels or in their combination. Upper limits at 95% confidence level on the production cross section of dark matter are set in a type-II two Higgs doublet model extended by a Z boson and in a baryonic Z model. The results in the baryonic Z model are also interpreted in terms of the spin-independent dark matter nucleon scattering cross section. This is the first search for DM particles produced in association with a Higgs boson decaying to a pair of W or Z bosons, and the first statistical combination based on five Higgs boson decay 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 centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.  [16] ATLAS Collaboration, "Search for an invisibly decaying Higgs boson or dark matter candidates produced in association with a Z boson in pp collisions at √ s = 13 TeV with the ATLAS detector", Phys. Lett. B 776 (2018) 318, doi:10.1016/j.physletb.2017.11.049, arXiv:1708.09624.