Search for invisible decays of a Higgs boson produced through vector boson fusion in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for invisible decays of a Higgs boson is performed using proton-proton collision data collected with the CMS detector at the LHC in 2016 at a center-of-mass energy $\sqrt{s} =$ 13 TeV, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The search targets the production of a Higgs boson via vector boson fusion. The data are found to be in agreement with the background contributions from standard model processes. An observed (expected) upper limit of 0.33 (0.25), at 95% confidence level, is placed on the branching fraction of the Higgs boson decay to invisible particles, assuming standard model production rates and a Higgs boson mass of 125.09 GeV. Results from a combination of this analysis and other direct searches for invisible decays of the Higgs boson, performed using data collected at $\sqrt{s}=$ 7, 8, and 13 TeV, are presented. An observed (expected) upper limit of 0.19 (0.15), at 95% confidence level, is set on the branching fraction of invisible decays of the Higgs boson. The combined limit represents the most stringent bound on the invisible branching fraction of the Higgs boson reported to date. This result is also interpreted in the context of Higgs-portal dark matter models, in which upper bounds are placed on the spin-independent dark-matter-nucleon scattering cross section.


Introduction
Since the discovery of the Higgs boson at the CERN LHC [1][2][3], the ATLAS and CMS Collaborations have pursued a wide-ranging program to study its properties and interactions. Precise measurements of the couplings of the Higgs boson to standard model (SM) particles indicate that the properties of the new particle are consistent with the SM predictions [4]. These measurements also provide indirect constraints on additional contributions to the Higgs boson width from beyond the SM (BSM) decays. Based on the results presented in Ref. [4], an indirect upper limit on the Higgs boson branching fraction to BSM particles of 0.34 is set at 95% confidence level (CL).
In the SM, the Higgs boson decays invisibly (H → inv) only through the H → ZZ → 4ν process, with a branching fraction, B(H → inv), of about 10 −3 . The rate for invisible decays of the Higgs boson may be significantly enhanced in the context of several BSM scenarios [5][6][7][8], including those in which the Higgs boson acts as a portal to dark matter (DM) [9][10][11][12]. Direct searches for H → inv decays increase the sensitivity to B(H → inv) beyond the indirect constraints. The ATLAS Collaboration [13] presented a combination of direct searches using √ s = 7 and 8 TeV data from proton-proton (pp) collisions, yielding an observed (expected) upper limit of 0.25 (0.27) on B(H → inv) at 95% CL [14]. The CMS Collaboration [15] performed a similar combination based on √ s = 7, 8, and 13 TeV pp collision data collected up to the end of 2015, setting an observed (expected) upper limit of 0.24 (0.23) on B(H → inv) at 95% CL [16].
This Letter presents a search for invisible decays of a Higgs boson, using pp collision data at √ s = 13 TeV collected with the CMS detector in 2016, corresponding to an integrated luminosity of 35.9 fb −1 . The search targets events in which a Higgs boson is produced in association with jets from vector boson fusion (VBF), as illustrated by the Feynman diagram in Fig. 1 (left). In these events, a Higgs boson is produced along with two jets that exhibit a large separation in pseudorapidity (|∆η jj |) and a large dijet invariant mass (m jj ). This characteristic signature allows for the suppression of SM backgrounds, making the VBF channel the most sensitive mode for invisible decays of a Higgs boson at hadron colliders [16,17]. The invisible particles produced by the Higgs boson decay can recoil with high transverse momentum (p T ) against the visible VBF-jet system, resulting in an event with large p T imbalance, which can be used to select signal enriched regions. In this phase space, the main expected backgrounds originate from Z(νν)+jets and W( ν)+jets processes. They are estimated from data using dedicated control regions (CRs), which consist of high purity samples of Z or W bosons decaying leptonically ( = µ, e). While earlier searches probing this final state at the LHC were based on counting experiments, the analysis presented in this Letter more optimally exploits the distinctive kinematic features of the VBF topology by fitting the shape of the m jj distribution. This approach is referred to as the "shape analysis". The shape analysis has been designed to provide a substantially improved sensitivity to invisible decays of the SM Higgs boson, resulting in the most sensitive VBF H → inv search reported to date. In addition, a simple but less sensitive counting approach, referred to as the "cut-and-count analysis", allows for an easier interpretation of the results of this search in the context of other phenomenological models predicting the same final-state signature. Upper limits on the product of the cross section and branching fraction for an additional Higgs boson with SM-like couplings, which does not mix with the 125 GeV Higgs boson, are also reported.
To further improve the sensitivity, results from a combination of searches for invisible decays of the Higgs boson, using data collected at √ s = 7, 8, and 13 TeV, are also presented. The searches considered in this combination target the VBF, the associated production (denoted by VH, where V denotes a W or Z boson), and the gluon fusion modes, whose representative Feynman diagrams are shown in Fig. 1. The VH-tag includes both a search for ZH production, in which the Z boson decays to a pair of leptons (e, µ) or b quarks, and one where a Lorentz-boosted W or Z boson decays to light-flavor quarks, whose corresponding hadronization products are reconstructed as a single large-radius jet. Additional sensitivity is achieved by including a search for gg → gH production (hereafter referred to as ggH), where a high-p T Higgs boson candidate is produced in association with jets from initial-state radiation. When these searches are combined to set an upper limit on B(H → inv), SM production cross sections are assumed. The result of this combination is also interpreted in the context of Higgs-portal models of DM interactions [9][10][11][12], in which the 125 GeV Higgs boson plays the role of a mediator between the SM and DM particles, thereby allowing for the possibility of producing DM candidates. This Letter is organized as follows: after a brief description of the CMS detector in Section 2, the event reconstruction in Section 3, and the simulated signal and background processes in Section 4, Section 5 is dedicated to the event selection requirements followed by a detailed description of the analysis strategy in Section 6. Section 7 reports the results of the VBF search in terms of upper limits on B(H → inv). Section 8 reports the upper limit on B(H → inv) from a combination of the aforementioned searches for invisible decays of the Higgs boson based on 13 TeV data collected in 2016 while, in Section 9, results from a more complete combination, involving also similar analyses performed on the 7 and 8 TeV data sets, are presented. The Letter is summarized in Section 10.

The CMS detector
The CMS detector is a multi-purpose apparatus designed to study a wide range of physics processes in both pp and heavy ion collisions. The central feature of the experiment is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T parallel to the beam direction. Within the solenoid volume a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL) are installed, each composed of a barrel and two endcap sections. The tracker system measures the momentum of charged particles up to |η| = 2.5, while the ECAL and HCAL provide coverage up to |η| = 3.0. In addition, the steel and quartz-fiber Cherenkov hadron forward calorimeter extends the coverage to |η| = 5.0. Muons are detected in gasionization chambers embedded in the steel flux-return yoke outside the solenoid, which cover up to |η| = 2.4. reconstruction optimized for fast processing, reducing the event rate to about 1 kHz.
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. [15].

Event reconstruction
The particle-flow (PF) algorithm [19] aims to reconstruct and identify each 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 momentum of the associated track at the primary interaction vertex, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The momentum of muons is obtained from the curvature of the corresponding tracks. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matched ECAL and HCAL energy deposits, corrected 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 energy.
The missing transverse momentum vector ( p miss T ) is computed as the negative vector p T sum of all the PF candidates in an event, and its magnitude is denoted as p miss T . Hadronic jets are reconstructed by clustering PF candidates using the anti-k T algorithm [20,21], with a distance parameter of 0.4. The reconstructed vertex with the largest value of summed physics object p 2 T is taken to be the primary pp interaction vertex, where physics objects correspond to the jets and the p miss T measured in the event. The charged PF candidates originating from any other vertex are ignored during the jet finding procedure. Jet momentum is determined as the vector sum of all particle momenta inside the jet, and is found from simulation to vary, on average, between 5 and 10% of the true momentum over the whole p T spectrum and detector acceptance. An offset correction is applied to jet energies to take into account the contribution from additional pp interactions within the same or adjacent bunch crossings (pileup) [22]. Jet energy corrections are derived from simulation and are confirmed with in situ measurements of the energy balance in dijet, multijet, γ+jets, and leptonically decaying Z+jets events [23]. These energy corrections are also propagated to the p miss T calculation [24].
Muon candidates, within the geometrical acceptance of the silicon tracker and muon subdetectors (|η| < 2.4), are reconstructed by combining the information from the tracker and the muon chambers [25]. These candidates are required to satisfy a set of quality criteria based on the number of hits measured in the tracker and the muon system, the properties of the fitted muon track, as well as the impact parameters of the track with respect to the primary vertex of the event.
Electron candidates within |η| < 2.5 are reconstructed using an algorithm that associates fitted tracks in the silicon tracker with electromagnetic energy clusters in the ECAL [26]. To reduce the misidentification rate, these candidates are required to satisfy identification criteria based on the shower shape of the energy deposit, the matching of the electron track to the ECAL energy cluster, the relative amount of energy deposited in the HCAL detector, and the consistency of the electron track with the primary vertex. Because of non-optimal reconstruction performance, electron candidates in the transition region between the ECAL barrel and endcap, 1.44 < |η| < 1.57, are not considered in the analysis. Electrons identified as coming from photon conversions in the detector are discarded [27].
Identified electrons or muons are required to be isolated from hadronic activity in the event. The isolation is defined by summing the p T of all the PF candidates within a cone of radius R = √ (∆η) 2 + (∆φ) 2 = 0.4 (0.3) around the muon (electron) track, and is corrected for the contribution of neutral hadrons from pileup interactions [25,26].
Hadronically decaying τ leptons (τ h ) are identified from reconstructed jets via the hadron-plusstrip algorithm [28], that requires a subset of particles inside the jet to be consistent with the decay products of a τ lepton. In addition, the τ h candidate must be isolated from other activity in the detector. The isolation is computed by summing the p T of all the charged PF candidates and PF photons within a cone of radius R = 0.3 around the jet axis. Hadronic τ leptons are selected with an average efficiency between 60 and 65%.

Simulated samples
The signal and background processes are simulated using several Monte Carlo (MC) generators. Higgs boson signal events, produced through ggH and VBF, are generated with POWHEG v2.0 [29][30][31][32][33] at next-to-leading order (NLO) approximation in perturbative quantum chromodynamics (QCD). Signal events are normalized to the inclusive Higgs boson production cross sections taken from the recommendations of Ref. [34]. The ggH production cross section is computed at next-to-next-to-NLO (N 3 LO) precision in QCD, and at NLO in electroweak (EW) theory [35]. The cross section for Higgs boson production through VBF is calculated at nextto-NLO (NNLO) in QCD, including also NLO EW corrections. The ggH process is simulated using calculations in which the top quark loop is fully resolved. The p T distribution of the Higgs boson produced via ggH is reweighted to match the NNLO plus next-to-next-toleading-logarithmic (NNLL) prediction from HRES v2.1 [36,37]. When upper limits are set on B(H → inv) for the SM Higgs boson, both ggH and VBF signal events are generated assuming a Higgs boson mass of 125.09 GeV, which is consistent with the combined ATLAS and CMS measurement [38] based on 7+8 TeV data, as well as the recent CMS measurement at 13 TeV in the H → ZZ → 4 channel [39].
The Z/γ * ( + − )+jets, Z(νν)+jets, and W( ν)+jets backgrounds are simulated at leading order (LO) using MADGRAPH5 aMC@NLO v2.2.2 [40], where up to four partons in the final state are included in the matrix element calculation. The background processes involving the production of a vector boson (V) in association with two jets exclusively through EW interactions, i.e. of order α 4 , are simulated at LO via MADGRAPH5 aMC@NLO. In addition, the QCD multijet background is also simulated at LO using MADGRAPH5 aMC@NLO. The tt and single top quark background samples are produced at NLO QCD using POWHEG v2.0 and v1.0, respectively [41][42][43]. Finally, the WZ and ZZ diboson productions are simulated at LO with PYTHIA v8.205 [44], while the Vγ and WW processes are simulated at NLO QCD using MAD-GRAPH5 aMC@NLO and POWHEG [45], respectively.
In all cases, generated events are interfaced with PYTHIA v8.205 or higher for the simulation of fragmentation, parton showering, and the underlying event description, using the parameters from the CUETP8M1 tune [46]. In the case of LO (NLO) MADGRAPH5 aMC@NLO samples, partons from the matrix elements are matched to the parton shower description via the MLM [47] (FxFx [48]) scheme. The NNPDF v3.0 [49] parton distribution functions (PDFs) are used for all the matrix element calculations. Interactions of the final-state particles with the CMS detector are simulated with GEANT4 [50]. Simulated events include the effects of pileup, and are weighted to reproduce the observed pileup distribution.

Event selection
Events in the signal region (SR) are selected initially by the L1 trigger exploting the p miss T information, whose threshold varies between 60 and 90 GeV depending on the instantaneous luminosity. The p miss T at the L1 trigger is computed from the vector p T sum of all the energy depositions in the calorimeters with |η| < 3. Partial mistiming of signals in the forward region of the ECAL endcaps (2.5 < |η| < 3.0) led to a reduction in the L1 trigger efficiency. A correction for this effect was determined using an unbiased data sample. This correction was found to be about 1% for m jj of 200 GeV and it increases to about 20% for m jj larger than 3.5 TeV.
At the HLT level, events of interest are collected using triggers with thresholds of 110 or 120 GeV, depending on the data taking period, applied equally to both the missing transverse momentum computed at the trigger level (p miss T, trig ) and the H miss T, trig variable. The H miss T, trig is defined as the magnitude of the vector p T sum of the reconstructed jets at the trigger level in the event with p T > 20 GeV and |η| < 5. The energy fraction attributed to neutral hadrons in jets with |η| < 3.0 is required to be less than 90%, in order to remove spurious jets originating from detector noise. Both p miss T, trig and H miss T, trig are calculated without including muon candidates, allowing the same triggers to be used also for selecting events in the muon CRs, which are used in the background estimation procedure described in Section 6.
Offline, events considered in the VBF search are required to have at least two jets with p T larger than 80 (40) GeV for the leading (subleading) jet. Since the L1 trigger decision does not use information from the hadronic activity in the forward region, at least one of the two leading jets in the event is required to have |η| < 3. To ensure a high and stable trigger efficiency, events are further required to have p miss T > 250 GeV. The trigger efficiency is measured as a function of H miss T , computed from jets with p T > 30 GeV and |η| < 3. After correcting for the L1 mistiming inefficiency, these triggers are found to be fully efficient for events passing the analysis selection with H miss T > 250 GeV. In addition, if the leading jet is within the geometrical acceptance of the tracker (|η| < 2.4), its energy fraction attributed to charged hadrons is required to be greater than 10%, while the energy fraction attributed to neutral hadrons is required to be smaller than 80%. These requirements, along with quality filters applied to tracks, muon candidates, and other physics objects, reduce the contamination arising from large misreconstructed p miss T from noncollision backgrounds [51]. To further suppress the contamination from QCD multijet events, in which a large p miss T may arise from a severe mismeasurement of the jet momentum, the jets in the event, with p T > 30 GeV and |η| < 4.7, are required to not be aligned with the p miss T . The minimum value of the azimuthal angle between the p miss T vector and each jet (min∆φ( p miss T , p jet T )) is required to be larger than 0.5 rad, where only the first four leading jets are included in the min∆φ( p miss T , p jet T ) definition. This selection reduces the QCD multijet contamination to less than 1% of the total background.
The two leading jets in VBF signal events typically show a large separation in η, large m jj and a small azimuthal separation (|∆φ jj |). The discriminatory power of |∆φ jj | results from a combination of the spin-parity properties of the Higgs boson and the high-p T regime explored by this search [17], in which the two VBF jets tend to recoil against the invisible system. The Z(νν)+jets and W( ν)+jets processes constitute the largest backgrounds in this search. The shape analysis primarily employs the large separation power of m jj to discriminate between VBF signal and V+jets backgrounds. Therefore, in this scenario, a set of loose requirements is applied on both m jj and |∆η jj |, i.e. |∆η jj | > 1.0 and m jj > 200 GeV. To further reduce the V+jets contamination, |∆φ jj | is required to be smaller than 1.5 rad and the two jets must lie in opposite hemispheres, η j1 η j2 < 0. In the cut-and-count approach, Z(νν)+jets and W( ν)+jets processes are suppressed by a more stringent event selection requiring |∆η jj | > 4.0 and m jj > 1.3 TeV, while the require-ment applied on |∆φ jj | remains unchanged.
The W( ν)+jets background is further suppressed by rejecting events that contain at least one isolated electron or muon with p T > 10 GeV, or a τ h candidate with p T > 18 GeV and |η| < 2.3, where the isolation is required to be less than 25 (16)% of the muon (electron) p T . With this strategy, prompt muons (electrons) are selected with an average efficiency of about 98 (95)%. In order to further reduce the contribution from γ+jets and Vγ processes, events containing an isolated photon with p T > 15 GeV and |η| < 2.5, passing identification criteria based on its ECAL shower shape [27], are vetoed.
Top quark backgrounds (tt and single top quark processes) are suppressed by rejecting events in which at least one jet, with p T > 20 GeV and |η| < 2.4, is identified as a b quark jet using the combined secondary vertex (CSVv2) algorithm [52]. A working point that yields a 60% efficiency for tagging a b quark jet and a 1 (10)% probability of misidentifying a light-flavor (c quark) jet as a b quark jet is used.
A summary of the selection criteria for the SR for both the shape and the cut-and-count analyses is shown in Table 1. Table 1: Summary of the kinematic selections used to define the SR for both the shape and the cut-and-count analyses.

Observable
Shape analysis Cut-and-count analysis Target background Leading (subleading) jet

Analysis strategy
The search exploits the large m jj and |∆η jj | that characterize events from VBF Higgs boson production. In the shape analysis case, the signal is extracted by fitting the sum of the signal and background shapes to the binned m jj distribution observed in data. The signal is expected to accumulate as an excess of events over the background at large values of m jj . This strategy necessitates a precise estimation of the shape of the background m jj distribution.
About 95% of the total expected background in this search is due to the V+jets processes, namely Z(νν)+jets and W( ν)+jets. A fraction of the V+jets background, referred to as V+jets (EW), can be attributed to the EW production of a Z or a W boson in association with two jets. A representative Feynman diagram contributing to V+jets (EW) production is shown in Fig. 2 (left). The remaining V+jets contribution arises from the production of a vector boson in association with QCD radiation, as shown in Fig. 2 (right). This is referred to as the V+jets (QCD) background. For both EW and QCD productions, the expected Z(νν)+jets rate in the SR is about two times larger than the W( ν)+jets contribution.
A comparison of the shapes of the key discriminating observables used in this analysis, obtained after applying the requirements listed in Table 1 except for those on m jj , |∆η jj | and |∆φ jj |, Figure 2: Representative leading order Feynman diagrams for the production of a Z boson in association with two partons arising from EW (left) and QCD (right) interactions. The left diagram contributes to the Z(νν)+jets (EW) production cross section, while the diagram on the right to the Z(νν)+jets (QCD) one. Diagrams for EW and QCD production of a W boson in association with two jets are similar to those reported above for the Z(νν)+jets process.
is shown in Fig. 3 for simulated signal and V+jets background events. From these distributions, it can be seen that the V+jets (EW) background is kinematically similar to the VBF Higgs boson signal. Therefore, its contribution to the total V+jets background rate increases when the two leading jets have large m jj and |∆η jj |. The V+jets (EW) process constitutes about 2% of the total V+jets background for m jj around 200 GeV. Its contribution increases to about 20% for m jj ≈ 1.5 TeV, and is more than 50% for m jj > 3 TeV.  Figure 3: Comparison between the shapes of the m jj (left), |∆η jj | (middle) and |∆φ jj | (right) distributions of signal events, produced by VBF (solid black) and ggH (dashed black) mechanisms, and V+jets backgrounds from both QCD (solid red) and EW (solid blue) production. Both signal and background distributions are scaled in order to have unit area. Distributions are obtained from simulated events passed through the CMS event reconstruction.

Overview of the V+jets background estimation
The Z+jets and W+jets backgrounds are estimated using four mutually exclusive CRs. These include a dimuon and a dielectron CR consisting mostly of Z( )+jets events that are kinematically similar to Z(νν)+jets background if the presence of the two leptons in the event is ignored. The W+jets background is estimated using CRs consisting of single-muon and singleelectron events stemming mainly from leptonic decays of a W boson. In contrast to the W+jets background in the SR, the single-lepton CRs consist of leptons that fall within the detector acceptance and pass the identification requirements. The p miss T in all the CRs is calculated by excluding the contribution of the identified leptons. Therefore, it corresponds to the p T of the hadronic recoil system, which resembles the p miss T expected from the V+jets backgrounds in the SR.
The event yield in the dilepton CRs is considerably smaller than the Z(νν)+jets contribution in the SR because the Z( ) branching fraction, where = µ or e, is six times smaller than the Z(νν) branching fraction. Consequently, the dilepton CRs have a limited statistical power to constrain the Z(νν)+jets background by themselves. In contrast, the yield of the single-lepton CRs is comparable to the Z(νν)+jets background. Furthermore, the Z(νν)+jets and W( ν)+jets processes are kinematically similar if the presence of the charged lepton is ignored. The theoretical uncertainties involved in the prediction of the Z+jets and W+jets cross sections largely cancel out in their ratio. Therefore, this ratio is predicted very reliably by the simulation and can be used as a constraint to connect the statistically rich single-lepton CRs to the Z(νν)+jets background in the SR.
The predictions for the V+jets processes obtained from simulation are referred to as "pre-fit" expectations, and are considered to be the initial estimates for the V+jets yields in the CRs and SR. These V+jets yields are then treated as freely floating parameters, and are fit to the data in all CRs and the SR. The V+jets yields obtained from this fit are referred to as "post-fit" estimates, and serve as the final V+jets background predictions in the analysis.

Definition of control regions
Dimuon and single-muon CRs are selected using the same L1 and HLT p miss T -based triggers that are used to collect events in the SR. Dimuon events are required to contain exactly two oppositely charged muons with p T > 10 GeV that form an invariant mass (m µµ ) between 60 and 120 GeV, which is compatible with a Z boson decay. Events with additional electrons or photons are rejected. At least one of the two muons must have p T > 20 GeV, and is required to pass tighter identification criteria based on the number of measurements in the tracker and the muon systems, the quality of the muon track fit, and the consistency of the muon track with the primary vertex. The isolation, as defined in Section 3, is required to be smaller than 15% of the muon p T . These tightly identified muons are selected with an average efficiency of 90%.
In the single-muon CR, events are required to contain exactly one muon with p T > 20 GeV, passing both tight identification and isolation requirements. The transverse mass (m T ) of the muon-p miss T is the p T of the muon, and ∆φ is the angle between p µ T and p miss T in the transverse plane. The m T is required to be smaller than 160 GeV, and no additional electrons or photons are allowed in the event.
Events in the dielectron and single-electron CRs are collected mainly using a single-electron trigger with a p T threshold of 27 GeV. In the case of dielectron events where the Z boson has p T > 600 GeV, the two electrons have a small angular separation, and are likely to get included in each other's isolation cones. This results in an inefficiency for the chosen trigger that imposes isolation requirements on electron candidates. This inefficiency is mitigated by including events collected by a single-electron trigger with a p T threshold of 105 GeV and no isolation requirements on the electron candidate.
The dielectron events are required to contain exactly two oppositely charged electrons with p T > 10 GeV and no additional muons or photons. As in the case of the dimuon events, the invariant mass of the dielectron system is required to be between 60 and 120 GeV. At least one of the two electrons must have p T > 40 GeV, and is required to pass a tight identification criterion based on the shower shape of its ECAL energy deposit, the matching of the electron track to the ECAL energy cluster, and the consistency of the electron track with the primary vertex. Furthermore, the isolation is required to be smaller than 6% of the electron p T . These selection requirements for electrons have an average efficiency of 70%.
Events in the single-electron CR are required to contain exactly one tightly identified and isolated electron with p T > 40 GeV; no additional muons or photons are allowed. The contamination from QCD multijet events is reduced by requiring p miss T > 60 GeV and m T < 160 GeV.
Events in the CRs must also satisfy the requirements imposed on events in the SR. When doing so, the negative p T of the hadronic recoil system is used instead of the p miss T in the event.

Estimation of V+jets backgrounds
The V+jets yields in the CRs are translated to the background estimates in the SR using transfer factors that are derived from simulation. The transfer factors are defined as the ratio of the yields of a given V+jets background in the SR and the corresponding process measured in each CR.
The transfer factors for the dilepton CRs account for the difference in the branching fractions of the Z(νν) and Z( ) decays, and the γ * ( ) contribution, as well as the impact of lepton acceptance and selection efficiencies. In the case of dielectron events, the transfer factors also account for the difference in trigger efficiencies. Transfer factors between the W( ν)+jets event yields in the single-lepton CRs and the W+jets background estimate in the SR take into account the effect of lepton acceptance, selection efficiencies, and lepton and τ h veto efficiencies, as well as the difference in trigger efficiencies in the case of the single-electron CR.
The constraint on the ratio of the cross sections of the Z+jets and W+jets processes, which is used to connect the single-lepton CRs to the Z(νν)+jets in the SR, is also implemented as a transfer factor, and is computed as the ratio of the Z(νν)+jets and W( ν)+jets yields in the SR. In order to have the most precise estimate of this constraint, the LO simulations for the Z+jets (QCD) and the W+jets (QCD) processes are corrected using boson p T and m jj -dependent NLO QCD K-factors derived with MADGRAPH5 aMC@NLO. The Z+jets and W+jets simulations are also corrected as a function of boson p T with NLO EW K-factors derived from theoretical calculations [53]. Similarly, Z+jets (EW) and W+jets (EW) processes are corrected with NLO QCD K-factors derived using the VBFNLO event generator [54,55] as a function of boson p T and m jj .
The V+jets background yields are determined using a maximum-likelihood fit, performed simultaneously across all CRs and the SR. The likelihood function is defined as: where P(x|y) = y x e −y /x!. The symbol i denotes each bin of the m jj distribution in the shape analysis, while, in the cut-and-count case, i stands for a single bin that represents the event yields obtained at the end of the event selection. The symbols d , and d i denote the observed number of events in each bin i of the dimuon, dielectron, single-muon, single-electron CRs, and the SR, respectively. The symbols f i (θ) Q and f i (θ) E indicate the ratios between the W( ν)+jets and Z(νν)+jets backgrounds in the SR from QCD and EW production, respectively.
and R e i (θ) Q are the transfer factors relating the dimuon, dielectron, single-muon, and single-electron CRs, respectively, to the SR for the V+jets (QCD) processes. Similarly, and R e i (θ) E indicate the transfer factors for the V+jets (EW) processes. The parameters κ i νν represent the yield of the Z(νν)+jets (QCD) background in each bin i of the SR, and are left to float freely in the fit. In a given bin, the Z(νν)+jets (EW) background yield is obtained from κ i νν through the transfer factor R Z i that represents the ratio between the Z(νν)+jets (QCD) and Z(νν)+jets (EW) processes. The contributions from subleading backgrounds in each region are estimated directly from simulation and they are denoted by B Systematic uncertainties are modeled as constrained nuisance parameters (θ), for which lognormal or Gaussian priors, indicated by P j (θ) in previous equation, are considered. The systematic uncertainties in the V+jets background estimates are introduced in the likelihood as variations of the transfer factors. These include theoretical uncertainties in the Z+jets to W+jets differential cross section ratio for both the QCD and EW processes due to the choice of the renormalization and the factorization scales, as well as the choice of the PDFs. The QCD scale variations are assumed to be uncorrelated between the Z+jets and W+jets processes, and therefore they do not cancel in the Z+jets to W+jets cross section ratio. This results in larger uncertainties compared to those from NLO calculations recommended in Ref. [53]. The uncertainty due to the choice of the renormalization scale varies between 8 and 12% as a function of m jj for both Z+jets/W+jets (QCD) and (EW) ratios. Similarly, the uncertainty due to the choice of the factorization scale varies between 2 and 7%. This also covers the uncertainty in the Z+jets/W+jets cross section ratio due to the interference between the V+jets (QCD) and V+jets (EW) processes, which is not included in the simulation. The PDF uncertainties are assumed to be correlated across V+jets processes, resulting in a residual uncertainty smaller than 1% on the Z+jets/W+jets cross section ratio. The uncertainties related to NLO EW corrections to the V+jets (QCD) processes are estimated according to the recommendations in Ref. [53], and are found to be about 1-2% across the entire m jj spectrum. Additional uncertainties included in the transfer factors include uncertainties in the reconstruction efficiencies of leptons (around 1% per muon or electron), the selection efficiencies of leptons (about 1% per muon, 1.5% per electron), the veto efficiency of leptons (around 2% per muon, 1% per electron) and τ h candidates (about 3.5% per τ h ), the knowledge of the jet energy scale (1-2%), and the efficiency of the electron (around 1%) and p miss T triggers (about 2%).
The full set of systematic uncertainties related to the V+jets transfer factors are listed in Table 2. Before any fit is performed, the total uncertainty in the expected background in the SR ranges between 4.5 and 6% as a function of m jj , dominated by the theoretical uncertainties in the Z+jets to W+jets cross section ratio for both QCD and EW production. The impact of each source of systematic uncertainty, as reported in Table 2 in the context of the shape analysis, is defined as the maximum difference in the fitted value of the signal strength, (σ/σ SM ) B(H → inv), obtained by varying the associated nuisance parameter within one standard deviation of its maximum likelihood estimate. In this procedure, the per-bin κ i νν parameters are profiled when a given nuisance parameter is shifted from its best fit estimate.
The m jj distributions in the dilepton and single-lepton CRs are shown in Fig. 4. The pre-fit predictions for the V+jets processes are shown in red. These indicate the level of agreement between data and simulation before a fit is performed. An estimate of the V+jets backgrounds is then obtained by fitting the data across all the CRs. This is depicted by the blue line in Fig. 4. This fit is referred to as the "CR-only" fit, since it does not impose any constraint on the V+jets yields due to the data in the SR.
To assess the level of agreement between data and simulation obtained through the application of p T -m jj dependent NLO corrections to both the V+jets (QCD) and V+jets (EW) processes, the ratio between the number of Z+jets and W+jets events in the CRs in bins of m jj is used as a benchmark. Figure 5 shows the ratio of the Z+jets and W+jets event yields in the muon (left) and electron (right) CRs, respectively. A good agreement is observed between data and simulation and local differences are covered by the systematic uncertainties listed in Table 2.

Other backgrounds
In addition to the V+jets processes, several other minor sources of background contribute to the total event yield in the SR. These include QCD multijet events that typically have small genuine p miss T . However, jet momentum mismeasurements and instrumental effects may give rise to large p miss   Table 2, as well as the statistical uncertainty in the simulation.
30% in the estimated QCD multijet background in the SR. The MC statistical uncertainty of the QCD multijet samples, which affects the transfer factor prediction, is also considered and is found to vary between 40 and 100% as a function of m jj . Lastly, a validation of the ∆φ method is performed using a purer sample of QCD multijet events that pass the analysis requirements, but have p miss T in the range of 100-175 GeV. In this validation region, the predicted QCD background is found to agree with the observation within 50%, which is taken as a conservative estimate of an additional uncertainty.
The remaining background sources include top quark production and diboson processes, which are estimated from simulation. The p T distribution of the top quark in simulation is corrected to match the observed p T distribution in data [57]. An uncertainty of about 10% is assigned to the overall top quark background normalization, while an additional 10% uncertainty is added to account for the modeling of the top quark p T distribution in simulation. The overall normalization of the diboson background has an uncertainty of about 15% [58,59]. The uncertainties in the top quark and diboson backgrounds are correlated across the SR and the CRs. Several experimental sources of uncertainty are also assigned to these backgrounds. An uncertainty of 2.5% in the integrated luminosity measurement [60] is propagated to the background yields. The uncertainty in the efficiency of the b quark jet veto is estimated to be around 3% for the top quark background and of about 1% for the other simulated processes. The uncertainty related to the jet energy scale varies between 8 and 15%, depending on both the process and the CR.

Results
This section presents the results obtained from the shape and the cut-and-count analyses. These include 95% CL upper limits on B(H → inv), and an interpretation of the search in the context of a BSM model which allows for the presence of a SM-like Higgs boson with a mass between 110 and 1000 GeV.

The shape analysis
The observed and the expected m jj distributions in the SR, obtained after applying the full event selection, are shown in Fig. 6. The background prediction shown in Fig. 6 (left) is obtained from a fit to the data in the CRs. Signal distributions for the SM Higgs boson produced via the ggH and VBF modes are overlaid, assuming B(H → inv) = 1. The estimated background yields from the CR-only fit are listed in Table 3, along with the observed event yield in the SR. The large contamination from ggH production arises from the low-m jj bins, which represent the least sensitive region to H → inv decays. Systematic uncertainties in the V+jets transfer factors and in the minor backgrounds introduce correlations across the m jj bins used in the fit. The correlations in the predicted background yields in each m jj bin are reported in Table 4.   Figure 6: The observed m jj distribution of the shape analysis SR compared to the post-fit backgrounds from various SM processes. On the left, the predicted backgrounds are obtained from a combined fit to the data in all the CRs, but excluding the SR. On the right, the predicted backgrounds are obtained from a combined fit to the data in all the CRs, as well as in the SR, assuming the absence of any signal. Expected signal distributions for a 125 GeV Higgs boson produced through ggH and VBF modes, and decaying to invisible particles with a branching fraction B(H → inv) = 1, are overlaid. The last bin includes all events with m jj > 3.5 TeV. The description of the ratio panels is the same as in Fig. 4.
An excess of about 4-10% is observed in the SR data when compared to the estimated backgrounds. The discrepancy resides mainly in the bulk of the m jj distribution. The shape of the excess is inconsistent with the characteristic features of a VBF signal, whose presence is expected to produce an increasing discrepancy between data and backgrounds as m jj increases. A goodness of fit test, based on a saturated χ 2 test statistic [61,62], yields a p-value of about 6% indicating that the data are compatible with the SM prediction. Figure 6 (right) shows the background prediction obtained after including events from the SR in the fit, but assuming the absence of a signal. Such a fit is referred to as the "b-only fit". Table 3: Expected event yields in each m jj bin for various background processes in the SR of the shape analysis. The background yields and the corresponding uncertainties are obtained after performing a combined fit across all the CRs, but excluding data in the SR. The "other backgrounds" includes QCD multijet and Z( )+jets processes. The expected total signal contribution for the 125 GeV Higgs boson, decaying to invisible particles with a branching fraction B(H → inv) = 1, and the observed event yields are also reported. The comparison between the results of the b-only fit with that allowing for the presence of the signal is used to set an upper limit on B(H → inv). In the b-only fit, the V+jets estimate in the SR can vary with respect to the prediction from the CRs within the systematic uncertainties assigned to the transfer factors. Therefore, the additional constraint due to the data in the SR mitigates the excess shown in Fig. 6 (left), yielding a p-value for the b-only fit of about 65%.
The results of this search are interpreted in terms of an upper limit on the product of the Higgs boson production cross section and its branching fraction to invisible particles, σB(H → inv), relative to the predicted cross section assuming SM interactions, σ SM . Observed and expected 95% CL upper limits are computed using an asymptotic approximation of the CL s method [63,64] with a profile likelihood ratio test statistic [65] in which systematic uncertainties are modeled as nuisance parameters following a frequentist approach [66]. The profile likelihood ratio is defined as: where (σ/σ SM )B(H → inv) represents the value of the signal strength that maximizes the likelihood L for the data, whileθ (κ) andθ a (κ a ) denote the best fit estimates for the nuisance parameters (Z(νν)+jets rate in each bin) and the estimates for a given fixed value of (σ/σ SM ) B(H → inv), respectively.
The relative contributions of the VBF and ggH production modes are fixed to the SM prediction within their uncertainties. The uncertainties in the predictions of the inclusive VBF and ggH production cross sections due to PDF uncertainties, renormalization and factorization scale variations are taken from Ref. [34]. An additional uncertainty of 40% is assigned to the expected ggH contribution. This accounts for both the limited knowledge of the ggH cross section in association with two or more jets, as well as the uncertainty in the prediction of the ggH differential cross section for large Higgs boson transverse momentum, p H T > 250 GeV. The former contribution is obtained by following the recipe outlined in Ref. [34] and is found to be about 30%, while the latter uncertainty is estimated by comparing the prediction from POWHEG+MINLO [67] with the one from MADGRAPH5 aMC@NLO and ranges between 20 and 25%. Furthermore, the uncertainties in the signal acceptance due to the choice of the PDF set are also evaluated independently for the different signal processes, and are treated as independent nuisance parameters in the fit.
The observed (expected) 95% CL upper limit on B(H → inv) is measured to be 0.33 (0.25), and the regions containing 68% and 95% of the distribution of upper limits, expected in absence of a signal, are found to be 0.18-0.35 and 0.14-0.47, respectively.
The background estimates reported in Table 3, along with the correlation matrix presented in Table 4, can be used in the simplified likelihood approach detailed in Ref. [68] to reinterpret these results in theoretical models different from those presented in this Letter. Table 4: Correlation between the uncertainties in predicted background yields across the m jj bins of the shape analysis SR. The backgrounds are estimated by fitting the data in the CRs. Bin ranges are expressed in TeV.

The cut-and-count analysis
The cut-and-count analysis is presented because it allows for an easier reinterpretation of this search in the context of other theoretical models that predict p miss T plus VBF dijet signatures. The observed event yield after the cut-and-count selection is reported in Table 5, along with the predicted backgrounds in the SR. The backgrounds are estimated by fitting the data in the CRs. An excess, characterized by a significance of about 2.5 standard deviations, is observed in the SR compared to the background prediction obtained from the CRs. As for the shape analysis, this excess is mostly due to low m jj events. The excess is incompatible with a VBF Higgs boson signal and, upon detailed scrutiny, is ascribed to a statistical fluctuation.
The results of the cut-and-count analysis are presented in terms of a 95% CL upper limit on B(H → inv) using the statistical procedure outlined in Section 7.1. The observed (expected) upper limit is found to be 0.58 (0.30), and the regions containing 68% and 95% of the distribution of upper limits, expected in absence of a signal, are found to be 0.22-0.43 and 0.17-0.58, respectively.

Constraints on a SM-like Higgs boson
The results presented are also interpreted in the context of an additional SM-like Higgs boson that does not mix with the 125 GeV boson and decays to invisible particles [69]. Such a boson may be produced via both the ggH and VBF mechanisms. This model has already been studied in earlier CMS publications [70][71][72]. Upper limits, computed at 95% CL on (σ/σ SM ) B(H → inv), are shown in Fig. 7 as a function of the SM-like Higgs boson mass hypothesis (m H ) for both the shape and the cut-and-count analyses. Table 5: Expected event yields in the SR and in the CRs of the cut-and-count analysis for various SM processes. The background yields and the corresponding uncertainties are obtained from a combined fit to data in all the CRs, but excluding data in the SR. The expected total signal contribution for the 125 GeV Higgs boson, decaying to invisible particles with a branching fraction B(H → inv) = 1, and the observed event yields are also reported.

Process
Signal region Dimuon CR Dielectron CR Single-muon CR Single-electron CR Z(νν) (QCD)   On the left, observed (solid black) and expected (dashed black) upper limits are obtained from the shape analysis while, on the right, results from the cut-and-count analysis are reported. The 68% (green) and 95% (yellow) CL intervals around the expected upper limits are also shown for both the shape and the cut-and-count analyses.

Combined limits on H → inv from 2016 data
The common feature of all the searches included in this combination is a large p miss T , where at least one high-p T jet or a weak boson recoils against the invisible particles produced by the Higgs boson decay. Specific topological selections are designed to reduce the contamination from large SM backgrounds, targeting a particular Higgs boson production mode. The analyses included in this combination are listed in Table 6, together with their expected signal composition and their individual upper limits on B(H → inv). The results quoted for the VBF channel come from the shape analysis described earlier in this Letter. The Z( + − )H analysis is identical to the one described in Ref. [72], where the expected signal comes entirely from invisible decays of the SM Higgs boson produced in association with a leptonically decaying Z boson, via either qq → ZH or gg → ZH production. In contrast, the V(qq')H and the ggH-tagged searches are similar to those described in Ref. [73], but events which overlap with the VBF analysis have been removed to avoid double counting. In both the ggH and V(qq')H searches, overlapping events represent about 6 (15)% of the total background for a p miss T of about 250 (1000) GeV. The overlap removal introduces a 5% loss in the expected exclusion sensitivity compared to that of Ref. [73]. Both the V(qq')H and the ggH searches target events with at least one high-p T central jet, and their SRs contain a mixture of different production modes. This mixture results from the limited discrimination power of the substructure observables exploited to select boosted V(qq')H candidates. No significant deviations from the SM expectations are observed in any of the searches. The results are interpreted as an upper limit on (σ/σ SM ) B(H → inv). These limits are calculated following the same approach described in Section 7.1. The combined likelihood fit accounts for correlations between the nuisance parameters in each search. The uncertainties in the diboson backgrounds (except for those considered in the Z( )H channel), tt and single top quark cross sections, lepton efficiencies, momentum scales, integrated luminosity, b quark jet and τ h vetoes are correlated among all the searches. In addition, the uncertainties in the inclusive signal production cross sections, due to renormalization and factorization scale variations, and the PDF uncertainties are also correlated across the channels. In contrast, since the jet kinematics in the VBF search differ from that in the other analyses, jet energy scale and resolution uncertainties are correlated only across the ggH and VH-tagged categories. The theoretical uncertainties applied to the V+jets (QCD) ratios are assumed to be uncorrelated between the VBF analysis and the other searches.
Observed and expected upper limits on (σ/σ SM ) B(H → inv) are computed at 95% CL and are presented in Fig. 8 (left). Assuming SM cross sections for each production mode, the combination yields an observed (expected) upper limit of B(H → inv) < 0.26 (0.20). The profile likelihood ratios as a function of B(H → inv), for both the combined fit and each individual search channel, are reported in Fig. 8 (right). Results are shown for both data and an Asimov dataset [65], defined by fixing the nuisances parameters to their maximum likelihood estimate obtained from a fit to the data in which B(H → inv) = 0 is assumed.

Combination of 7, 8, and 13 TeV searches for H → inv decays
The analyses previously described and listed in Table 6 are further combined with earlier searches performed using data collected at √ s = 7, 8, and 13 TeV up to the end of 2015, as reported in Refs. [16,70,74]. The 7 and 8 TeV data, collected in 2011 and 2012, correspond to integrated luminosities of up to 4.9 and 19.7 fb −1 [75,76], respectively. The 13 TeV data collected in 2015 correspond to an integrated luminosity of 2.3 fb −1 [77]. Systematic uncertainties in the inclusive ggH, VBF, and VH production cross sections are fully correlated across the 7, 8, and 13 TeV analyses. The uncertainty in the prediction of the Higgs boson p T distribution in ggH production, included in both the ggH and VH channels, and those arising from the limited knowledge of the PDFs, are also correlated among all searches. The uncertainties in the lepton and photon reconstruction and identification efficiencies, in the lepton momentum scales, and in the veto efficiency of leptons, τ h candidates, and b quark jets are uncorrelated between 7+8 and 13 TeV searches. Similarly, uncertainties in the jet energy scale and resolution, mistag rate of leptons, and modeling of the unclustered particles are also uncorrelated between 7+8 and 13 TeV searches. The b jet energy scale and resolution uncertainties for the Z(bb)H analysis are estimated using different techniques, and therefore are treated as uncorrelated with other channels. Theoretical uncertainties affecting the ratio of Z(νν) and W( ν) predictions in the VBF searches, for both QCD and EW V+jets processes, are uncorrelated across data sets because different strategies are followed to quantify and assign these uncertainties. In contrast, those affecting the Z(νν)/W( ν) and Z(νν)/γ+jets ratios in the ggH and V(qq)H channels are correlated across the 7+8 and the 13 TeV (2015 data) searches, as described in Ref. [16]. The uncertainties in the tune of the underlying event simulation and in the pileup modeling are uncorrelated between 7+8 and 13 TeV searches. Finally, theoretical uncertainties affecting diboson and top quark cross sections, except for those considered on the ZZ/WZ ratio in the Z( )H channel, are correlated across all data sets.
Observed and expected upper limits on (σ/σ SM ) B(H → inv) at 95% CL are presented in Fig. 9 (left). Limits are computed for the combination of all data sets, as well as for partial combinations based either on 7+8 or 13 TeV data. The relative contributions from different Higgs production mechanisms are constrained to their SM values within the theoretical uncertainties. The combination yields an observed (expected) upper limit of B(H → inv) < 0.19 (0.15) at 95% CL. The corresponding profile likelihood ratios as a function of B(H → inv) are shown in Fig. 9 (right). The measured value of the invisible branching fraction and an approximate 68% CL interval, obtained from the profile likelihood, are B(H → inv) = 0.05 ± 0.03 (stat) ± 0.07 (syst).  Figure 9: On the left, observed and expected 95% CL upper limits on (σ/σ SM ) B(H → inv) for partial combinations based either on 7+8 or 13 TeV data as well as their combination, assuming SM production cross sections for the Higgs boson with mass of 125.09 GeV. On the right, the corresponding profile likelihood ratios as a function of B(H → inv) are presented. The solid curves represent the observations in data, while the dashed lines represent the expected result obtained from the background-only hypothesis.
The relative sensitivity of each search considered in the combination depends on the assumed SM production rates. The cross sections for the ggH, VBF and VH production modes are parametrized in terms of coupling strength modifiers κ V and κ F , which directly scale the coupling of the Higgs boson to vector bosons and fermions, respectively [69]. The contribution from the gg → ZH production is scaled to account for the interference between the tH and ZH diagrams, as described in Ref. [34]. In this context, SM production rates are obtained for κ V = κ F = 1. Figure 10 (left) shows the observed 95% CL upper limits on (σ/σ SM ) B(H → inv) evaluated as a function of κ V and κ F . The LHC best estimates for κ V and κ F from Ref. [4] are superimposed, along with the 68% and 95% CL limit contours. Within the 95% CL region, the observed (expected) upper limit on B(H → inv) varies between 0.14 (0.11) and 0.24 (0.19).
The upper limit on B(H → inv), obtained from the combination of √ s = 7, 8, and 13 TeV searches, is interpreted in the context of Higgs-portal models of DM interactions, in which a stable DM particle couples to the SM Higgs boson. Direct-detection experiments are sensitive to the interaction between a DM particle and an atomic nucleus, which may be mediated by the exchange of a Higgs boson, producing nuclear recoil signatures that can be interpreted in terms of the DM-nucleon scattering cross section. The sensitivity of these experiments depends mainly on the DM particle mass (m χ ). If m χ is smaller than half of the Higgs boson mass, the Higgs boson invisible width (Γ inv ) can be translated, within an effective field theory approach, into a spin-independent DM-nucleon elastic scattering cross section, as outlined in Ref. [9]. This translation is performed assuming that the DM candidate is either a scalar or a Majorana fermion, and both the central value and the uncertainty of the dimensionless nuclear formfactor f N are taken from the recommendations of Ref. [78]. The conversion from B(H → inv) to Γ inv uses the relation B(H → inv) = Γ inv /(Γ SM + Γ inv ), where Γ SM is set to 4.07 MeV [69]. Since renormalizable models predicting a vectorial DM candidate require an extended dark Higgs sector, which may lead to modifications of kinematic distributions assumed for the invisible Higgs boson signal, such interpretation is not provided in the context of this Letter. Figure 10 (right) shows the 90% CL upper limits on the spin-independent DM-nucleon scattering cross section as a function of m χ , for both the scalar and the fermion DM scenarios. These limits are computed at 90% CL so that they can be compared with those from direct detection experiments such as XENON1T [79], LUX [80], PandaX-II [81], CDMSlite [82], CRESST-II [83], and CDEX-10 [84] which provide the strongest constraints in the m χ range probed by this search. In the context of Higgs-portal models, the result presented in this Letter provides the most stringent limits for m χ smaller than 18 (7)    GeV, whose production cross section varies as a function of the coupling modifiers κ V and κ F . Their best estimate, along with the 68% and 95% CL contours from Ref. [4], are also reported. The SM prediction corresponds to κ V = κ F = 1. On the right, 90% CL upper limits on the spin-independent DM-nucleon scattering cross section in Higgsportal models, assuming a scalar (solid orange) or fermion (dashed red) DM candidate. Limits are computed as a function of m χ and are compared to those from the XENON1T [79], LUX [80], PandaX-II [81], CDMSlite [82], CRESST-II [83], and CDEX-10 [84] experiments.

Summary
A search for invisible decays of a Higgs boson is presented using proton-proton (pp) collision data at a center-of-mass energy √ s = 13 TeV, collected by the CMS experiment in 2016 and corresponding to an integrated luminosity of 35.9 fb −1 . The search targets events in which a Higgs boson is produced through vector boson fusion (VBF). The data are found to be consistent with the predicted standard model (SM) backgrounds. An observed (expected) upper limit of 0.33 (0.25) is set, at 95% confidence level (CL), on the branching fraction of the Higgs boson decay to invisible particles, B(H → inv), by means of a binned likelihood fit to the dijet mass distribution. In addition, upper limits are set on the product of the cross section and branching fraction of an SM-like Higgs boson, with mass ranging between 110 and 1000 GeV.
A combination of CMS searches for the Higgs boson decaying to invisible particles, using pp collision data collected at √ s = 7, 8, and 13 TeV (2015 and 2016), is also presented. The combination includes searches targeting Higgs boson production via VBF, in association with a vector boson (with hadronic decays of the W boson and hadronic or leptonic decays of the Z boson) and via gluon fusion with initial state radiation. The VBF search is the most sensitive channel involved in the combination. No significant deviations from the SM predictions are observed in any of these searches. The combination yields an observed (expected) upper limit on B(H → inv) of 0.19 (0.15) at 95% CL, assuming SM production rates for the Higgs boson and a Higgs boson mass of 125.09 GeV. The observed 90% CL upper limit of B(H → inv) < 0.16 is interpreted in terms of Higgs-portal models of dark matter (DM) interactions. Constraints are placed on the spin-independent DM-nucleon interaction cross section. When compared to the upper bounds from direct detection experiments, this limit provides the strongest constraints on fermion (scalar) DM particles with masses smaller than about 18 (7) GeV.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus [24] CMS Collaboration, "Performance of the CMS missing transverse momentum reconstruction in pp data at √ s = 8 TeV", JINST 10 (2015) P02006, doi:10.1088/1748-0221/10/02/P02006, arXiv:1411.0511.