Search for dark matter produced in association with a Higgs boson decaying to $\gamma\gamma$ or $\tau^+\tau^-$ at $\sqrt{s} =$ 13 TeV

A search for dark matter particles is performed by looking for events with large transverse momentum imbalance and a recoiling Higgs boson decaying to either a pair of photons or a pair of $\tau$ leptons. The search is based on proton-proton collision data at a center-of-mass energy of 13 TeV collected at the CERN LHC in 2016 and corresponding to an integrated luminosity of 35.9 fb$^{-1}$. No significant excess over the expected standard model background is observed. Upper limits at 95% confidence level are presented for the product of the production cross section and branching fraction in the context of two benchmark simplified models. For the Z'-two-Higgs-doublet model (where Z' is a new massive boson mediator) with an intermediate heavy pseudoscalar particle of mass $m_\mathrm{A} =$ 300 GeV and $m_\mathrm{DM} =$ 100 GeV, Z' masses from 550 GeV up to 1265 GeV are excluded. For a baryonic Z' model, with $m_\mathrm{DM} =$ 1 GeV, Z' masses up to 615 GeV are excluded. Results are also presented for the spin-independent cross section for the dark matter-nucleon interaction as a function of the mass of the dark matter particle. This is the first search for dark matter particles produced in association with a Higgs boson decaying to two $\tau$ leptons.


Introduction
Astrophysical evidence strongly suggests the existence of dark matter (DM) in the universe [1]. Whether the DM has a particle origin remains a mystery [2]. There are a number of wellmotivated theories beyond the standard model (SM) of particle physics that predict the existence of a particle, χ, that could serve as a DM candidate. To date, only gravitational interactions between DM and SM particles have been observed. However, the discovery of a Higgs boson by both the ATLAS and CMS Collaborations at the CERN LHC in 2012 [3][4][5] provides a new way to probe DM-SM particle interactions.
Collider experiment searches have typically looked for DM recoiling against an associated SM particle. Since any produced DM is unlikely to interact with the detector material, it creates an imbalance in the recorded momentum yielding a large amount of missing transverse momentum, p miss T . This paper presents a search for DM recoiling against an SM-like Higgs boson (h) using the h + p miss T signature. This SM-like Higgs boson can be produced from initial-or finalstate radiation, or from a new interaction between DM and SM particles. However, initial-state radiation of an SM-like Higgs boson from a quark or gluon is suppressed by Yukawa or loop processes, respectively [6][7][8].
Previous searches for h + p miss T have been performed at both the ATLAS and CMS experiments. No excesses were observed in either h → bb or h → γγ decay channels with 20.3 fb −1 of data at √ s = 8 TeV [9, 10] or with 2.3-36.1 fb −1 of data at √ s = 13 TeV [11][12][13]. This paper examines two Higgs boson decay channels: h → γγ and h → τ + τ − . This is the first search for DM produced in association with h → τ + τ − and the first to combine results from a combination of the γγ and τ + τ − decay channels.
Two simplified models for DM+h production are used as benchmarks for this search, both of which were recommended by the ATLAS-CMS Dark Matter Forum [14]. The leading order (LO) Feynman diagrams for these models are shown in Fig. 1. The first benchmark model ( Fig. 1 left) is a Z -two-Higgs-doublet model (Z -2HDM) [7]. In this scenario, the SM is extended by a U(1) Z group, with a new massive Z boson mediator, while a Type-2 2HDM framework [15,16] is used to formulate the extended Higgs sector. At LO, the Z boson is produced resonantly and decays into an SM-like Higgs boson and an intermediate heavy pseudoscalar particle (A). The A then decays into a pair of Dirac fermionic DM particles. This analysis does not consider the contribution of the decay Z → Zh which can have a h + p miss T signature if Z → νν. The second diagram ( Fig. 1 right) describes a baryonic Z model [8]. In this scenario, a new massive vector mediator Z emits a Higgs boson and then decays to a pair of Dirac fermionic DM particles. Here, the baryonic gauge boson Z arises from a new U(1) B baryon number symmetry. A baryonic Higgs boson (h B ) is introduced to spontaneously break the new symmetry and generates the Z boson mass via a coupling that is dependent on the h B vacuum expectation value. The Z couplings to quarks and DM are proportional to the U(1) B gauge couplings. There is a mixing between h B and the SM Higgs boson, allowing the Z to radiate an SM-like Higgs boson. The stable baryonic states in this model are the candidate DM particles.
In the Z -2HDM, there are several parameters that affect the predicted cross section. However, when the A is on-shell, only the Z and A masses affect the kinematic distributions of the final state particles studied in this analysis. This paper considers a Z resonance with mass between 450 and 2000 GeV and an A pseudoscalar with mass between 300 and 700 GeV. Masses of A below 300 GeV are excluded by constraints on flavor-changing neutral currents from b → sγ measurements [16] and are not considered in this analysis. The ratio of the vacuum expectation values of the Higgs doublets (tan β) in this model is fixed to 1. As given in Ref.
[12], the DM particle mass is fixed to 100 GeV, the DM-A coupling strength g DM is fixed to 1, and the Z coupling strength g Z is fixed to 0.8.
For the baryonic Z model, this paper considers a Z resonance with a mass between 10 and 2500 GeV and DM particle masses between 1 and 900 GeV. As suggested for this model [17], the mediator-DM coupling is fixed to 1 and the mediator-quark coupling (g q ) is fixed to 0.25. The mixing angle between the baryonic Higgs boson and the SM-like Higgs boson is set to 0.3 and the coupling between the Z boson and the SM-like Higgs boson is proportional to the mass of the Z boson.
For both models, values of the couplings and mixing angle are chosen to maximize the predicted cross section. Results for other values can be obtained by rescaling the cross section since these parameters do not affect the kinematic distributions of the final state particles. The SM-like Higgs boson is assumed to be the already observed 125 GeV Higgs boson, since the SM-like Higgs boson has similar properties to the SM Higgs boson. Therefore, in this paper the observed 125 GeV Higgs boson is denoted by h.
Although the SM Higgs boson branching fractions to γγ and τ + τ − are smaller than the branching fraction to bb, the analysis presented here exploits these two decay channels because they have unique advantages compared with the h → bb channel. The h → γγ channel benefits from higher precision in reconstructed invariant mass and the h → τ + τ − channel benefits from smaller SM background. Additionally, the h → γγ and h → τ + τ − channels are not dependent on p miss T trigger thresholds, as such searches in these channels are complementary to those in the h → bb channel since they can probe DM scenarios with lower p miss T . The search in the h → γγ channel uses a fit in the diphoton invariant mass spectrum to extract the signal yield. In addition to a high-p miss T category, a low-p miss T category is also considered to extend the phase space of the search. In the h → τ + τ − channel, the three decay channels of the τ lepton with the highest branching fractions are analyzed. After requiring an amount of p miss T in order to sufficiently suppress the quantum chromodynamic (QCD) multijet background, the signal is extracted by performing a simultaneous fit to the transverse mass of the p miss T and the two τ lepton candidates in the signal region (SR) and control regions (CRs).
The paper is organized as follows. Section 2 gives a brief description of the CMS detector and the event reconstruction. Section 3 details the data set and the simulated samples used in the analysis. Then Sections 4 and 5 present the event selection and analysis strategy for each decay channel, respectively. The systematic uncertainties affecting the analysis are presented in Section 6. Section 7 details the results of the analysis and their interpretations. A summary is given in Section 8.

The CMS detector and event reconstruction
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 solenoid volume are a silicon pixel and strip tracker, a lead-tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL). Extensive forward calorimetry complements the coverage provided by the barrel and endcap detectors. Charged particle trajectories are measured by the silicon pixel and strip tracker system, covering 0 ≤ φ ≤ 2π in azimuth and |η| < 2.50, where the pseudorapidity is η = − ln (tan θ/2), and θ is the polar angle with respect to the counterclockwise beam direction. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke. A more detailed description of the CMS detector can be found in Ref. [18].
Events of interest are selected using a two-tiered trigger system [19]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
Using information from all CMS subdetectors, a global event reconstruction is performed using the particle-flow (PF) algorithm [20]. The PF algorithm optimally combines all of the detector information and generates a list of stable particles (PF candidates), namely photons, electrons, muons, and charged and neutral hadrons. The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex (PV). The physics objects are the jets, clustered using the jet finding algorithm [21,22] with the tracks assigned to the vertex as inputs, and the negative vector sum of the p T of those jets. The PV is used as the reference vertex for all objects reconstructed with the PF algorithm.
Photons are reconstructed from their energy deposits in the ECAL, which can involve several crystals [23]. A photon that converts to an electron-positron pair in the tracker will yield a shower spread out in azimuth due to the deflection of the electron and positron in the strong magnetic field. In order to achieve the best photon energy resolution, corrections are applied to overcome energy losses including those from photon conversions [23]. Additional corrections, calculated from the mass distribution of Z → e + e − events, are applied to the measured energy scale of the photons in data (≤1%) and to the energy resolution in simulation (≤2%).
Electron reconstruction requires the matching of the cluster of energy deposits in the ECAL with a track in the silicon tracker. Electron identification is based on the ECAL shower shape, matching between the track and ECAL cluster, and consistency with the PV. Muons are reconstructed by combining two complementary algorithms [24]: one that matches tracks in the silicon tracker with signals in the muon system, and another in which a global track fit seeded by the muon track segment is performed.
Jets are reconstructed from PF candidates using the anti-k T clustering algorithm [21] as implemented in FASTJET [22] with a distance parameter of 0.4. Jet energy corrections are derived from simulation to bring the average measured response of jets to that of particle-level jets. Hadronically decaying τ leptons are reconstructed from jets using the hadrons-plus-strips (HPS) algorithm [25]. The HPS algorithm uses combinations of reconstructed charged hadrons and energy deposits in the ECAL to reconstruct the τ lepton's three most common hadronic decay modes: 1-prong, 1-prong + π 0 (s), and 3-prong. In the h → τ + τ − channel, events with jets originating from b quark decays are excluded in order to reduce the background from tt events. The combined secondary vertex algorithm [26] is used to identify jets originating from b quarks by their characteristic displaced vertices.
The missing transverse momentum vector ( p miss T ), with magnitude p miss T , is the negative vector sum of the p T of all PF candidates in an event. Jet energy corrections are propagated to the p miss T for a more accurate measurement [27]. Events may have anomalously large p miss T from sources such as detector noise, cosmic ray muons, and beam halo particles, which are not well modeled in simulation. Event filters [28] are applied to remove such events.

Observed and simulated data samples
The analysis is performed with pp collision data at √ s = 13 TeV collected with the CMS detector in 2016. The data correspond to an integrated luminosity of 35.9 fb −1 .
The analysis strategy and event selection were optimized using Monte Carlo (MC) simulated samples of associated DM+h production via the two benchmark models discussed in Section 1. The MADGRAPH5 aMC@NLO v2.3.0 [29] generator is used to generate both the Z -2HDM and baryonic Z signals at LO. The decay of the SM-like Higgs boson is simulated by PYTHIA 8.205 [30].
A small but irreducible background for both decay channels in this analysis comes from events in which the SM Higgs boson is produced in association with a Z boson that decays to two neutrinos. Other SM Higgs boson production mechanisms are associated with resonant but reducible backgrounds. These include gluon-gluon fusion (ggh), vector boson fusion (VBF), and production in association with a pair of top quarks (tth). The production in association with a vector boson (Vh) and other SM Higgs boson backgrounds are all generated using MAD-GRAPH5 aMC@NLO v2.2.2 at next-to-leading order (NLO) in perturbative QCD.
The dominant nonresonant backgrounds for the h → γγ channel are events with mismeasured p miss T and two photons that happen to have an invariant mass close to the mass of the SM Higgs boson. The largest contributions to this are nonresonant γγ, γ + jet, and QCD multijet production. The simulated γγ sample is generated at LO with SHERPA v2.2.2 [31] while the γ + jet and QCD multijet samples are modeled at LO with PYTHIA. Additional backgrounds originate from electroweak (EW) processes such as single top, tt, W, or Z boson production in association with one or two photons, and Drell-Yan (DY) production where the Z boson decays to pairs of electrons or muons. The DY and all other EW backgrounds considered in the analysis are generated at NLO with MADGRAPH5 aMC@NLO. These nonresonant background samples are used for optimizing the analysis selection, however, they are not used for the ultimate background estimation.
The largest backgrounds for the h → τ + τ − channel are W + jets, tt, and multiboson processes. The MADGRAPH5 aMC@NLO v2.3.0 generator is used for W + jets processes, which are generated at LO in perturbative QCD with the MLM jet matching and merging scheme [32]. A p T -dependent correction factor is applied to the W + jets sample to account for next-to-nextto-leading order QCD and NLO EW effects [33][34][35][36]. The tt process is generated at NLO with the POWHEG 2.0 [37][38][39][40] generator. Single top quark production is modeled at NLO with the POWHEG 1.0 [41] generator. The FxFx [42] merging scheme is used to generate some smaller diboson backgrounds (including WZ samples) with the MADGRAPH5 aMC@NLO generator at NLO, while the dominant diboson backgrounds, WW and ZZ in two lepton final states, are generated using POWHEG 2.0. Another reducible background considered in this analysis is Z/γ * → /ττ, where is e or µ. The Drell-Yan background which is corrected for differences in the dilepton mass m /ττ and dilepton transverse momentum p T ( /ττ) distributions using dimuon events in data [43].
All simulated samples mentioned above use the NNPDF 3.0 parton distribution function (PDF) sets [44,45] with the order matching that used in the matrix element calculations. For parton showering and hadronization, as well as for τ lepton decays, the samples are interfaced with PYTHIA using the CUETP8M1 tune [46] for all samples except tt, for which the M2 tune is used. The MC samples are processed through a full simulation of the CMS detector based on GEANT4 [47] and are reconstructed with the same algorithms that are used for the data. All samples include the simulation of additional inelastic pp interactions in the same or neighboring bunch crossings (pileup). Minimum-bias collision events generated with PYTHIA are added to the simulated samples to reproduce the pileup effects in the data. Additionally, the simulated events are weighted so that the pileup vertex distribution matches that of the data, with an average of 27 interactions per bunch crossing.

Analysis strategy in the h → γγ channel
The search for DM+h in the h → γγ channel is performed by selecting events with two photons and a large amount of p miss T . The set of requirements detailed in Section 4.1 is applied to select well-identified photons and to enhance the signal significance. A fit to the diphoton invariant mass distribution, described in Section 4.2, is performed to extract the background and signal yields.

Event selection
The events used in this analysis were selected by a diphoton trigger with asymmetric p T thresholds of 30 and 18 GeV. The trigger also has loose photon identification criteria based on the cluster shower shape, isolation requirements, and a selection on the ratio of hadronic to electromagnetic energy deposits of the photon candidates.
The photons that enter the analysis are required to fall within the fiducial range of the ECAL (|η| < 1.44 or 1.57 < |η| < 2.50) and to satisfy various preselection criteria that are slightly more stringent than the trigger requirements. An additional veto on the presence of a track pointing to the ECAL cluster is applied to reject electrons that could be reconstructed as photons. Scale factors, extracted from Z → e + e − events using the tag-and-probe method [48], are applied to the simulated samples to account for any discrepancy in identification efficiency between data and simulation.
The isolation variables that are used in the photon identification are calculated by summing the p T of PF photons, neutral hadrons, or charged hadrons associated with the PV in a cone of radius ∆R = The isolation variables are corrected by the median transverse momentum density of the event to mitigate the effects of pileup [49]. Some of the signals considered can have Lorentz-boosted topologies. For example, high-mass mediators could result in a large boost to the Higgs boson. When a boosted Higgs boson decays to two photons, the resulting photons hit the ECAL close to each other. This effect leads to large contributions from one photon to the photon isolation sum of the other. In order to maintain high efficiency for high-mass mediator signals, the photon isolation requirement is not applied to photons that are within ∆R < 0.3 of each other. that pass the preselection were used to study the discriminating power of variables such as p miss T , the p T of the diphoton system p Tγγ , and the ratio p T /m γγ for each photon. A selection on p T that scales with m γγ is chosen so that it does not distort the shape of the m γγ distribution. The p Tγγ variable is included in the selection because it has higher resolution than the event's measured p miss T and is expected to be large for signal events, since the Higgs boson is produced back-to-back with p miss T . A high-p miss T category (p miss T ≥ 130 GeV) is optimal for the two benchmark models presented in this paper. A low-p miss T category (50 < p miss T < 130 GeV), optimized using as reference the baryonic Z signal model, is also included to probe softer signals, namely signals that may not be observed in other h + p miss T searches because they rely heavily on p miss T for background rejection. The chosen requirements, found to optimize the signal sensitivity for both models in the low-and high-p miss T categories, are given in Table 1.
Further background rejection is achieved using two topological requirements. The azimuthal separation |∆φ(p γγ , p miss T )| between p miss T and the Higgs boson direction reconstructed from the two photons must be greater than 2.1 to select events in which the Higgs boson and p miss T are back-to-back. Events with highly energetic jets collinear to p miss T are removed by the requirement that the min|∆φ(p jet , p miss T )| be greater than 0.5 for any jet with p T above 50 GeV. This rejects events with a large misreconstructed p miss T arising from mismeasured jet p T . Finally, events are vetoed if they have three or more jets each with p T above 30 GeV, to reject multijet backgrounds while maintaining a high efficiency for the two benchmark signal models. The p miss T distribution of the selected events is shown in Fig. 2.

Background estimation and signal extraction
A narrow resonance search similar to the SM Higgs boson diphoton analysis of Ref. [50] is performed. The diphoton invariant mass between 105 and 180 GeV is fit with a model that is the sum of the signal and background shapes. The signal shape, taken from the simulated events, is allowed to change independently in each of the two p miss T categories. The background shape includes a smooth function, estimated from the data, to model the continuum background, and a resonant contribution from the SM Higgs boson. The fit is performed with an unbinned maximum-likelihood technique in both the low-and high-p miss T categories discussed above.
The resonant background, arising from the SM Higgs boson decays to two photons, appears as a peak under the expected signal peak. This contribution from all SM Higgs boson production modes is estimated with the simulated events by including a mass distribution template, scaled to the NLO cross section, as a resonant component in the final fitting probability density function (pdf).
The nonresonant background contribution, mostly due to γγ and to various EW processes, is estimated using data. The nonresonant diphoton m γγ distribution in the data is fit, in each p miss T category, with an analytic function. Because the exact functional form of the background is unknown, the parametric model must be flexible enough to describe a variety of potential underlying functions. Using an incorrect background model can lead to biases in the measured signal yield that can artificially modify the sensitivity of the analysis. Three functions are considered as possible models for the nonresonant background; they are analytical forms that are  Table 1. Events with p miss T below 50 GeV are not used in the analysis. The cross sections of the signals are set to 1 pb. The total simulated background is normalized to the integral of the data. The statistical uncertainty in the total background is shown by the hatched bands. The data-to-simulation ratio is shown in the lower panel.
frequently used in dijet [51] and diphoton [52] resonance searches. The best functional form found to fit the nonresonant diphoton m γγ distribution, in both p miss T categories, is a power law function f (x) = ax −b where a and b are free parameters constrained to be positive.
A detailed bias study has been performed in order to choose this function. The m γγ shape of the simulated nonresonant events is used as a template to generate 1000 pseudo-experiments for each p miss T category. For each pseudo-experiment, the number of events generated is equal to the number of events observed in data in that category. The resulting m γγ distribution is fit with each analytic function considered. The exercise is also repeated injecting a potential signal contribution. The pulls of each pseudo-experiment, defined as the difference in the number of simulated events and those predicted by the fit function divided by the statistical uncertainties of the fit, are calculated. If the bias (the median of the pulls) is five times smaller than the statistical uncertainty in the number of fitted signal events, any potential bias from the choice of background model is considered negligible. Since this criterion is satisfied for the power law function, any systematic uncertainty in the bias from the background fit function is neglected in this analysis.  (right) categories, with the sum of a power law (dashed black) fit function to describe the nonresonant contribution, and a resonant shape (dashed red), taken from simulation, to take into account the SM h → γγ contribution. The SM h contribution is fixed to the theoretical prediction in the statistical analysis. The sum of the nonresonant and resonant shapes (solid blue) is used to estimate the total background in this analysis.

Event selection
The three final states of τ lepton pairs with the highest ττ branching fractions (eτ h , µτ h , and τ h τ h ) are considered in this analysis. In the eτ h and µτ h channels, one of the τ leptons decays leptonically to an electron or a muon and two neutrinos, while the other τ lepton decays hadronically (τ h ). In the third channel, τ h τ h , both τ leptons decay hadronically. The eµ, ee, and µµ final states are not included because of the low branching fraction of the ττ pair to purely leptonic final states. The ee and µµ final states are not considered, since they are overwhelmed by DY background.
Triggers based on the presence of a single electron (muon) are used to select events in the eτ h (µτ h ) channel. In the τ h τ h channel, the triggers require the presence of two isolated τ h objects. Each τ h candidate reconstructed offline is required to match a τ h candidate at the trigger level, with a ∆R separation less than 0.5.
The electrons and muons in the eτ h and µτ h channels are required to have p T greater than 26 GeV, exceeding the trigger thresholds for the single-electron and single-muon triggers. Electrons (muons) with |η| < 2.1 (2.4) are used. An eτ h (µτ h ) event is required to have an electron (muon) passing a multivariate MVA identification discriminator (straightforward selection identification criteria) and an isolation requirement of I e rel < 0.10 (I µ rel < 0.15), where I rel is defined as in Eq. (1), with an isolation cone of size ∆R = 0.3 (0.4) surrounding the electron (muon): Here ∑ p charged T , ∑ p neutral T , and ∑ p γ T are the scalar sums of transverse momentum from charged hadrons associated with the primary vertex, neutral hadrons, and photons, respectively. The term ∑ p PU T is the sum of transverse momentum of charged hadrons not associated with the primary vertex and p T is the p T of the electron or muon.
Hadronically decaying τ leptons in all channels are required to satisfy a loose (τ h τ h channel) or a tight (eτ h and µτ h channels) working point of an MVA isolation measure. The loose (tight) working point corresponds to a 65 (50)% efficiency with a 0.8 (0.2)% misidentification probability. The τ leptons are required to be identified as decaying via one of the three modes recognized with the HPS algorithm, and also pass discriminators that reduce the rate of electrons and muons misreconstructed as τ h candidates [53]. For the eτ h and µτ h channels, the τ h candidates are required to have p T > 20 GeV and |η| < 2.3. In the τ h τ h channel, the leading (subleading) τ lepton p T is required to be greater than 55 (40) GeV, both τ h momenta exceeding the double-hadronic τ lepton trigger thresholds of 35 GeV. The selection criteria are summarized in Table 2 for all three final states.
The p miss T is further required to be greater than 105 GeV and the visible p T of the ττ system is required to be greater than 65 GeV. These stringent criteria reduce the need for tighter isolation in the τ h τ h channel. Additionally, the mass reconstructed from the visible p T of the ττ system is required to be less than 125 GeV, to ensure that the ττ system is compatible with an SM Higgs boson. In order to minimize diboson and W + jets contributions, the two τ lepton candidates must pass a loose collinearity criterion of ∆R ττ < 2.0.
Two types of event veto are employed for background reduction. Events with jets tagged as originating from hadronization of b quarks are vetoed, to reduce tt and single top processes. The working point used in the b tagging algorithm corresponds to about a 66% efficiency for a 1% misidentification probability. In addition, events with additional muons or electrons beyond those from the τ lepton candidates are discarded, to reduce the contribution of multilepton backgrounds.

Signal extraction and background estimation
The signal is extracted from a maximum-likelihood fit to the total transverse mass (M tot T ) distributions in the different channels for the SR, and for the W + jets and QCD multijet background CRs. The M tot T is defined as: where p miss The W + jets and the QCD multijet background are estimated directly from the data. The procedure to estimate these processes relies on CRs, which are included in the maximum-likelihood fit, to extract the results. The other backgrounds, tt, Z + jets, SM Higgs boson, single top quark, and diboson production processes, are extracted from simulation.
The shape of the M tot T distribution of the W + jets background is estimated from simulation by requiring the same selection as for the SR, but the isolation of the τ lepton candidates is relaxed to increase the statistical precision of the distribution. To constrain the normalization of the W + jets background, a CR enriched in W + jets events is constructed by inverting the isolation criteria on the τ h candidates while keeping a loose isolation. The CR obtained by inverting the isolation criterion is included in the maximum-likelihood fit to constrain the normalization of the W + jets background in the SR.
To estimate the QCD multijet background, a CR in data is obtained by requiring the τ lepton candidates to have the same sign. No significant amount of signal and of background with opposite-sign τ h is expected in this CR because the τ h charge misidentification is of order 1% and the charge misidentification for electrons and muons is even smaller. All simulated backgrounds are subtracted from observed events in the CR, and the remaining contribution is classified as QCD multijet background. The contribution of QCD multijet events with oppositesign τ lepton candidates in the SR is obtained by multiplying the QCD multijet background, obtained in the same-sign CR, by a scale factor. The scale factor, approximately unity with an uncertainty of 20%, is determined from events with τ h candidates failing the isolation requirement and with low p miss T , which do not overlap with events selected in the SR. To increase the statistical precision of the QCD multijet distribution, the isolation of the τ lepton candidates is relaxed for the τ h channels, while conserving the normalization obtained as detailed above. The same-sign τ lepton candidate CR constrains the QCD multijet background normalization in the SR and the other CRs in the maximum-likelihood fit.
The normalizations of the W + jets and QCD multijet background are strongly correlated since both processes contribute to both CRs. The simultaneous fit of the SR and CRs takes into account this correlation. The SR distributions included in the simultaneous maximum-likelihood fit are shown in Fig. 4. For the signal extraction, W + jets and QCD multijet background CRs are considered separately, whereas in Fig. 4, the two backgrounds are presented merged together.

Systematic uncertainties
In both analysis channels, an uncertainty of 2.5% is used for the normalization of simulated samples to reflect the uncertainty in the integrated luminosity measurement in 2016 [54]. Common to both analysis channels are systematic uncertainties related to the theoretical production cross section of the Higgs boson. The PDF, and renormalization and factorization scale uncertainties are addressed using the recommendations of PDF4LHC [55] and LHC Higgs Cross Section [56] working groups, respectively. The value of these uncertainties range from 0.3 to 9.0%. The systematic uncertainties associated with each of the analysis channels are detailed below. Uncertainties affecting normalizations are represented by log-normal pdfs in the statistical analysis.

The h → γγ channel
In the h → γγ channel, there are several sources of experimental and theoretical uncertainties that affect the signal and the SM h → γγ yields. However, the largest source of uncertainty is statistical. As mentioned in Section 4.2, no systematic uncertainties are applied to the nonresonant background, which is extracted from a fit to data, since the bias of the fit is negligible compared to the statistical uncertainty of the data set. The systematic uncertainties for the h → γγ channel are summarized in Table 3. In addition to the theoretical uncertainties mentioned above, a 20% cross section uncertainty is included for the ggh sample, based on the Observed Bkg. uncertainty CMS differential measurements of h → γγ, for diphoton p T above 70 GeV [57]. The branching fraction uncertainty [56] of 1.73% is also included.
In addition to the integrated luminosity uncertainty, several other experimental sources of systematic uncertainty are included in this analysis. The trigger efficiency uncertainty (approximately 1%) is extracted from Z → e + e − events using a tag-and-probe technique [48]. The photon identification uncertainty of 2% arises from the observed difference in efficiencies between data and simulation. A 0.5% energy scale uncertainty is assigned to take into account the knowledge of the photon energy scale at the Z boson mass peak and its extrapolation to the Higgs boson mass. Additionally, several p miss T -related uncertainties are applied. The systematic uncertainty from mismeasured p miss T is evaluated by comparing the tail of the p miss T distributions in data and simulation in a γ + jet enriched CR. The efficiencies with which data and simulated events pass the p miss T selection are compared. The difference in efficiency is 50% and is included as a systematic uncertainty associated with mismeasured p miss T . However, the contribution of simulated backgrounds with mismeasured p miss T is quite small since only the ggh and VBF SM h → γγ production modes contribute. Finally, a systematic uncertainty, which is less than 4%, is applied to take into account the difference in efficiency between data and simulation when applying the topological ∆φ requirements in the low-p miss T region. This uncertainty is evaluated using Z → e + e − events, and only affects the ggh and VBF simulated samples.

The h → τ + τ − channel
The systematic uncertainties in the h → τ + τ − channel are related to the normalization of signal and background processes and, in several instances, the shapes of the signal and background distributions. As mentioned earlier, the simultaneous maximum-likelihood fit is performed in the SR and CRs, where the shape and normalization uncertainties are represented by nuisance parameters in the likelihood. Uncertainties affecting the distribution of M tot T (shape uncertainties) are represented by Gaussian pdfs, whereas log-normal pdfs are used for normalization, as stated above. The largest overall uncertainty is statistical. Table 4 summarizes the different sources of systematic uncertainty in this channel.
An uncertainty of 2% is assigned to simulated events containing an electron or muon candidate. In simulated events with a τ h candidate, an additional uncertainty of 5% per τ h is applied. These uncertainties account for the observed differences in the performance of electron, muon, and τ h identification, isolation, and trigger algorithms, between data and simulation. The hadronic τ h efficiency is not fully correlated across all ττ final states because there are different discriminators used in each channel. The τ h τ h channel has a 9% τ h uncertainty due to a correlation with the τ h τ h trigger systematic uncertainty. An uncertainty of 12% is assigned to simulated events containing an electron misidentified as a τ h candidate, and 25% for a muon misidentified as a τ h candidate [43]. A 2 (4)% uncertainty is assigned to the yield of multiboson and single top (tt) processes to account for changes in overall normalization arising from uncertainties in the b tagging performance. Similarly, a 5% b tagging uncertainty is assigned to Z + jets and SM Higgs boson processes, while all other processes, including signal, receive a 2% uncertainty. A systematic uncertainty of up to 20% is applied to QCD multijet background to account for yield differences in the same-sign CR. All of the background systematic uncertainties in the same-sign region are propagated to the total QCD multijet background uncertainty, which is taken to be 40%.
The W + jets background has a p T -dependent uncertainty, which approaches 10%, from predicted NLO EW K-factors where the full EW correction is treated as the systematic uncertainty [33][34][35][36]. Cross section uncertainties of the order of 5% are applied to the tt (6%), top quark (5%), and diboson (5%) processes [58][59][60][61]. In simulated Z + jets samples, a shape uncertainty of 10% of the Z boson p T reweighting correction, to account for higher-order effects, is used. The uncertainty in the Z + jets background contribution is about 12% in the SR. The tt contribution includes a shape systematic uncertainty equivalent to 5% related to the top quark p T spectrum, since there is evidence that the spectrum is softer in data than in simulation [61]. A 1.2% uncertainty in the τ lepton energy scale [43] is propagated through to the final signal extraction variables. The τ lepton energy scale depends on the τ h decay mode and is correlated across all channels. A shape uncertainty is used for the uncertainty in the double τ h trigger. A shift of 3% of the p T of the trigger-level τ h candidate leads to a 12% normalization difference at 40 GeV, and a 2% difference at 60 GeV. For p T > 60 GeV, a constant 2% systematic uncertainty is applied.
To account for potentially different rates of jets misidentified as τ h candidates between data and simulation, an uncertainty, applied as a function of the p T of the τ h candidate, is used for background events where the reconstructed τ h candidate is matched to a jet at generator-level. The uncertainty increases to about 20% near a τ h candidate p T = 200 GeV, and acts to change the shape of the M tot T distribution. An asymmetric uncertainty related to the identification of τ h with a high p T is applied to signal and background simulation. The high-p T τ h efficiency measurement uses selected highly virtual W bosons and has limited statistical precision in comparison to the lower p T Z → ττ and tt τ h efficiency studies. Therefore the asymmetric uncertainty is used in combination with a constant scale factor. It is proportional to p T and has a value of +5% and −35% per 1 TeV. For the application of all of the aforementioned τ lepton uncertainties, simulated backgrounds are separated depending on whether the reconstructed τ h candidates are matched to generated τ leptons.
In all simulated samples, uncertainties in the p miss T calculation related to unclustered energy deposits are taken into account. Uncertainties in the jet energy scale are included on an eventby-event basis and propagated to the p miss T calculation. Lastly, an uncertainty in the statistical precision of each process in each bin of the distribution is also included.

Results
The results of the analysis are derived from the maximum-likelihood fits presented in Sections 4 and 5 for the h → γγ and h → τ + τ − channels, respectively.

Observed yields
For the h → γγ channel, from the signal plus background fit to m γγ , the number of expected events from background processes are determined. The background yields and the observed number of events within 3 GeV of the SM Higgs boson mass are listed in Table 5 for both the low-and high-p miss T categories. The excess at low-p miss T has negligible effect on the results when combined with the high-p miss T category for the benchmark signals considered in this paper. Table 5: Expected background yields and observed numbers of events for the h → γγ channel in the m γγ range of 122-128 GeV are shown for the low-and high-p miss T categories. The nonresonant background includes QCD multijet, γγ, γ + jet, and EW backgrounds and is estimated from the analytic function fit to data. The SM Higgs boson background is presented separately for the irreducible Vh production and for the other production modes. For the resonant background contributions, both the statistical and the systematic uncertainties are listed. As detailed in Section 4.2, the systematic uncertainty associated with the nonresonant background is negligible.

Expected background
Low-p miss T category High-p miss T category SM h → γγ (Vh) 2.9 ± 0.1 (stat) ± 0.2 (syst) 1.26 ± 0.05 (stat) ± 0.09 (syst) SM h → γγ (ggh, tth, VBF) 5.3 ± 0.3 (stat) ± 1.2 (syst) 0.11 ± 0.01 (stat) ± 0.01 (syst) Nonresonant background 125.1 ± 11.2 (stat) 4.5 ± 2.1 (stat) Total background 133 ± 11 (stat) ± 1 (syst) 5.9 ± 2.1 (stat) ± 0.1 (syst) Observed events 159 6 In the h → τ + τ − channel, the final simultaneous fit to the M tot T distributions for the SR, and W + jets and QCD multijet CRs is performed in each of the three considered τ decay channels (eτ h , µτ h , and τ h τ h ). The extracted post-fit yields for the expected number of background events and the number of events observed in data are shown in Table 6. The number of events observed is in good agreement with the number of events predicted by the SM backgrounds. Aside from the small excess in the low-p miss T category of the h → γγ channel, the observed numbers of events are consistent with SM expectations. All of the results presented here are interpreted in terms of the two benchmark models of DM production mentioned earlier. Expected signal yields and the product of the predicted signal acceptances and their efficiencies (A ) are summarized in Table 7 for selected mass points, in both h → γγ and h → τ + τ − channels. Table 7: The expected signal yields and the product of acceptance and efficiency (A ) for the two benchmark models. The Z -2HDM signal is shown for the parameters m A = 300 GeV and m Z = 1000 GeV, and the baryonic Z signal, for the parameters m DM = 1 GeV and m Z = 10 GeV.

Interpretation in the Z -2HDM model
For the event selection given in Sections 4 and 5, the results interpreted in terms of the Z -2HDM associated production of DM and a Higgs boson are presented here. The expected and observed yields are used to calculate an upper limit on the production cross section of DM+h production via the Z -2HDM mechanism. Upper limits are computed [62] at 95% confidence level (CL) using a profile likelihood ratio and the modified frequentist criterion [63,64] with an asymptotic approximation [65]. The upper limits are obtained for each Higgs boson decay channel separately and for the statistical combination of the two. The two decay channels are combined using the Higgs boson branching fractions predicted by the SM [56]. In the combination of the two analyses, the theoretical uncertainties in the Higgs boson cross section and the systematic uncertainty in the integrated luminosity are assumed to be fully correlated between the two decay channels. Figure 5 shows the 95% CL expected and observed upper limits on the DM production cross section (σ 95% CL ) as a function of Z mass. Both the h → γγ and h → τ + τ − channels, as well as the combination of the two, are shown for m A = 300 GeV. These upper limits, although obtained with a DM mass of 100 GeV, can be considered valid for any DM mass below 100 GeV since the branching fraction for decays of A to DM particles decreases as the dark matter mass increases. The theoretical cross section (σ th ) is calculated with m DM = 100 GeV, g Z = 0.8, and g DM = tan β = 1, as mentioned in Section 1.
To produce exclusion limits in the two-dimensional plane of Z mass and A mass, an interpolation is performed. Fully simulated signal samples (mentioned in Section 3) were generated in a coarse grid of m A and m Z . For the h → γγ channel, the m γγ shape does not depend on the mass of these particles, only the expected yield is affected by these masses. Therefore, the product A of the fully simulated samples is parametrized and used to extract the expected number of events for intermediate mass points. In the h → τ + τ − channel, this is not sufficient because the M tot T shape does depend on the particle masses. A reweighting technique is used to extract the yields for the intermediate mass points. Simulation samples were produced at generator-level for m Z between 450 and 2000 GeV in steps of 50 GeV and for m A between 300 and 700 GeV in steps of 25 GeV. These are compared with the full-simulation samples at generator-level. The bin-by-bin ratio of the SM-like Higgs boson p T between the two samples is used to weight the full-simulation samples. This method was validated by applying the same procedure at the generator-level among the samples for which full-simulation is available.
The interpolation between mass points is improved using kernel algorithms to display smooth, continuous exclusion contours. The resulting two-dimensional exclusion for the Z -2HDM signal is shown in Fig. 6. The 95% CL expected and observed upper limits on signal strength (σ 95% CL /σ th ) are shown. Regions of the parameter space with σ 95% CL /σ th < 1 are excluded at 95% CL under the nominal σ th hypothesis. For m A = 300 GeV, the h → γγ channel alone excludes at 95% CL Z masses up to 860 GeV, while the h → τ + τ − channel excludes m Z up to 1200 GeV. The combination of these two decay channels excludes Z masses up to 1265 GeV for m A = 300 GeV. The Z mass range considered is extended from previous CMS searches to include 450 ≤ m Z < 600 GeV.

Baryonic Z model interpretation
Here the results presented in Section 7.1 are interpreted in the context of the baryonic Z model. This paper presents the first baryonic Z model interpretation of h + p miss T searches with the CMS detector. The 95% CL upper limits on DM+h cross section are calculated for the baryonic Z production mechanism. The upper limits for each decay channel and the combination of the two channels are shown in Fig. 7. The σ th is calculated assuming the choice of parameters detailed in Section 1. Results in the two-dimensional plane of m DM and m Z are produced using an interpolated grid produced in the same way as described in Section 7.2. The twodimensional exclusion for this model is shown in Fig. 8, where the 95% CL upper limits on the signal strength are shown for each decay channel and for the combination of the h → γγ and h → τ + τ − channels. For m DM = 1 GeV, the h → γγ channel excludes m Z masses up to 574 GeV. The h → τ + τ − channel similarly excludes m Z masses up to 450 GeV. The combination of the two decay channels excludes m Z up to 615 GeV for m DM = 1 GeV.       Figure 9: The 90% CL exclusion limits on the DM-nucleon SI scattering cross section as a function of m DM . Results obtained in this analysis are compared with those from a selection of direct detection (DD) experiments. The latter exclude the regions above the curves. Limits from CDMSLite [66], LUX [67], XENON-1T [68], PandaX-II [69], and CRESST-II [70] are shown.

Simplified DM model interpretation
Limits from the baryonic Z model are reinterpreted to infer limits on the s-channel simplified DM models that were proposed by the ATLAS-CMS Dark Matter Forum [14] for comparison with direct detection experiments. In the model considered in this analysis, Dirac DM particles couple to a vector mediator, which in turn couples to the SM quarks. A point in the parameter space of this model is determined by four variables: the DM particle mass m DM , the mediator mass m med , the mediator-DM coupling g DM , and the universal mediator-quark coupling g q . The couplings for this analysis are fixed to g DM = 1.0 and g q = 0.25, following the recommendation of Ref. [17].
The results are interpreted in the spin-independent (SI) cross section σ SI for DM scattering off a nucleus. The value of σ SI for a given point in the s-channel simplified DM model is determined by the equation [17]: where µ nDM is the reduced mass of the DM-nucleon system and f (g q ) is the mediator-nucleon coupling, which is dependent on g q . The resulting σ SI limits as a function of DM mass are shown in Fig. 9. In the same plot, exclusions from several direct detection experiments are shown. For the baryonic Z model, the limits are more stringent than direct detection experiments for m DM < 2.5 GeV.

Summary
A search for dark matter particles produced in association with a Higgs boson has been performed. The study focuses on the case where the 125 GeV Higgs boson decays to either two photons or two τ leptons. This analysis is based on proton-proton collision data collected with the CMS detector during 2016 at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 . The results of the search are interpreted in terms of a Z -two-Higgs-doublet model (Z -2HDM) and a baryonic Z simplified model of dark matter production.
A statistical combination of the two channels was performed and these results were used to produce upper limits on dark matter production. Limits on the signal production cross section are calculated for both simplified models. For the Z -2HDM signal, with an intermediate pseudoscalar of mass m A = 300 GeV and m DM = 100 GeV, Z masses up to 1265 GeV are excluded at 95% confidence level. For the baryonic Z model, with m DM = 1 GeV, Z masses up to 615 GeV are excluded. This is the first search for dark matter produced in association with a Higgs boson decaying to two τ leptons and the first to combine results from the γγ and τ + τ − decay channels. The Z -2HDM interpretation extended the Z mass range compared with previous CMS searches. The interpretation of the results include the first baryonic Z model interpretation for CMS and an extrapolation to limits on the spin-independent cross section for the dark matter-nucleon interaction.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies:  [9] ATLAS Collaboration, "Search for dark matter produced in association with a Higgs boson decaying to two bottom quarks in pp collisions at √ s = 8 TeV with the ATLAS detector", Phys. Rev. D 93 (2016) 47, doi:10.1103/PhysRevD.93.072007, arXiv:1510.06218.