Search for dark matter produced in association with a leptonically decaying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Z}} $$\end{document}Z boson in proton–proton collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\,\text {Te}\text {V} $$\end{document}s=13Te

A search for dark matter particles is performed using events with a Z boson candidate and large missing transverse momentum. The analysis is based on proton–proton collision data at a center-of-mass energy of 13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {Te}\text {V}$$\end{document}Te, collected by the CMS experiment at the LHC in 2016–2018, corresponding to an integrated luminosity of 137\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {fb}^{-1}$$\end{document}fb-1. The search uses the decay channels \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Z}} \rightarrow {\mathrm{e}} {\mathrm{e}} $$\end{document}Z→ee and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Z}} \rightarrow {{\upmu }{}{}} {{\upmu }{}{}} $$\end{document}Z→μμ. No significant excess of events is observed over the background expected from the standard model. Limits are set on dark matter particle production in the context of simplified models with vector, axial-vector, scalar, and pseudoscalar mediators, as well as on a two-Higgs-doublet model with an additional pseudoscalar mediator. In addition, limits are provided for spin-dependent and spin-independent scattering cross sections and are compared to those from direct-detection experiments. The results are also interpreted in the context of models of invisible Higgs boson decays, unparticles, and large extra dimensions.


Introduction
The existence of dark matter (DM) is well established from astrophysical observations [1], where the evidence relies entirely on gravitational interactions. According to fits based on the Lambda cold dark matter model of cosmology [2] to observational data, DM comprises 26.4% of the current matter-energy density of the universe, while baryonic matter accounts for only 4.8% [3]. In spite of the abundance of DM, its nature remains unknown. This mystery is the subject of an active experimental program to search for dark matter particles, including direct-detection experiments that search for interactions of ambient DM with ordinary matter, indirect-detection experiments that search for the products of self-annihilation of DM in outer space, and searches at accelerators and colliders that attempt to create DM in the laboratory. e-mail: cms-publication-committee-chair@cern.ch The search presented here considers a "mono-Z" scenario where a Z boson, produced in proton-proton (pp) collisions, recoils against DM or other beyond the standard model (BSM) invisible particles. The Z boson subsequently decays into two charged leptons ( + − , where = e or μ) yielding a dilepton signature, and the accompanying undetected particles contribute to missing transverse momentum. The analysis is based on a data set of pp collisions at a center-ofmass energy of 13 TeV produced at the CERN LHC. The data were recorded with the CMS detector in the years 2016-2018, and correspond to an integrated luminosity of 137 fb −1 . The results are interpreted in the context of several models for DM production, as well as for two other scenarios of BSM physics that also predict invisible particles.
These results extend and supersede a previous search by CMS in the mono-Z channel based on a data set collected at √ s = 13 TeV corresponding to an integrated luminosity of 36 fb −1 [4]. The ATLAS experiment has published searches in this channel as well with the latest result based on a data set corresponding to an integrated luminosity of 36 fb −1 [5]. Similar searches for DM use other "mono-X" signatures with missing transverse momentum recoiling against a hadronic jet [6,7], a photon [8], a heavy-flavor (bottom or top) quark [9-11], a W or Z boson decaying to hadrons [5,7,12], or a Higgs boson [13][14][15][16][17][18]]. An additional DM interpretation is explored in searches for Higgs boson decays to invisible particles [19,20].
The paper is organized as follows. The DM and other BSM models explored are introduced along with their relevant parameters in Sect. 2. Section 3 gives a brief description of the CMS detector. The data and simulated samples are described in Sect. 4, along with the event reconstruction. The event selection procedures and background estimation methods are described in Sects. 5 and 6, respectively. Section 7 details the fitting method implemented for the different models presented, while Sect. 8 discusses the systematic uncertainties. The results are given in Sect. 9, and the paper is summarized in Sect. 10.

Signal models
Several models of BSM physics can lead to a signature of a Z boson subsequently decaying into a lepton pair and missing transverse momentum. The goal of this paper is to explore a set of benchmark models for the production of DM that can contribute to this final state. In all DM models we consider, the DM particles are produced in pairs, χχ, where χ is assumed to be a Dirac fermion.
First, we consider a set of simplified models for DM production [21,22]. These models describe the phenomenology of DM production at the LHC with a small number of parameters and provide a standard for comparing and combining results from different search channels. Each model contains a massive mediator exchanged in the s-channel, where the mediator (either a vector, axial-vector, scalar, or pseudoscalar particle) couples directly to quarks and to the DM particle χ. An example tree-level diagram is shown in Fig. 1 (upper left). The free parameters of each model are the mass of the DM particle m χ , the mass of the mediator m med , the mediator-quark coupling g q , and the mediator-DM coupling g χ . Following the suggestions in Ref. [22], for the vector and axial-vector studies, we fix the couplings to values of g q = 0.25 and g χ = 1 and vary the values of m χ and m med , and for the scalar and pseudoscalar studies, we fix the couplings g q = 1 and g χ = 1, set the dark matter particle mass to m χ = 1 GeV, and vary the values of m med . The comparison with data is carried out separately for each of the four spin-parity choices for the mediator.
We also explore a two-Higgs-doublet model (2HDM) with an additional pseudoscalar boson, a, that serves as the mediator between DM and ordinary matter. This "2HDM+ a" model [23,24] is a gauge-invariant and renormalizable model that contains a Higgs scalar (h), which we take to be the observed 125 GeV Higgs boson, a heavy neutral Higgs scalar (H), a charged Higgs scalar (H±), and two pseudoscalars (A, a), where the pseudoscalar bosons couple to the DM particles. For the process studied in this paper, the H boson is produced via gluon fusion and decays into a standard model (SM) Z boson and the pseudoscalar a. These subsequently decay into a pair of leptons and a pair of DM particles, respectively, as shown in Fig. 1 (upper right). The sizable couplings of the Z boson to the Higgs bosons makes the mono-Z channel more sensitive to this model than the mono-jet or monophoton channels. Among the parameters of this model are the Higgs boson masses, the ratio tan β of the vacuum expectation values of the two Higgs doublets, and the mixing angle θ of the pseudoscalars. We consider only configurations in which m H = m H± = m A , and fix the values tan β = 1 and sin θ = 0.35, following the recommendations of Ref. [24].
We also examine the case where the h boson acts as a mediator for DM production, as discussed in "Higgs portal" models [25][26][27][28]. If m χ < m h /2, the Higgs boson could decay invisibly into a pair of DM particles. The mechanism for such decays can be found, for example, in many supersymmetric theoretical models that contain a stable neutral lightest supersymmetric particle, e.g., a neutralino [29], that is sufficiently light. An illustrative Feynman diagram for such a case is shown in Fig. 1 (lower left), while additional gluon-induced diagrams are also considered.
In addition to the DM paradigm, we consider a model where unparticles are responsible for the missing transverse momentum in the final state. The unparticle physics concept [30,31] is based on scale invariance, which is anticipated in many BSM physics scenarios [32][33][34]. The effects of the scale-invariant sector ("unparticles") appear as a non-integral production. Here A represents the DM mediator, χ represents a DM particle, while (H, h) and a represent the scalar and pseudoscalar Higgs bosons, respectively. Here h is identified with the 125 GeV scalar boson. The dotted line represents either an unparticle or a graviton number of invisible massless particles. In this scenario, the SM is extended by introducing a scale-invariant Banks-Zaks field, which has a nontrivial infrared fixed point [35]. This field can interact with the SM particles by exchanging heavy particles with a high mass scale M U [36]. Below this mass scale, where the coupling is nonrenormalizable, the interaction is suppressed by powers of M U and can be treated within an effective field theory (EFT). The parameters that characterize the unparticle model are the possible noninteger scaling dimension of the unparticle operator d U , the coupling of the unparticles to SM fields λ, and the cutoff scale of the EFT Λ U . In order to remain in the EFT regime, the cutoff scale is set to Λ U = 15 TeV and to maintain unitarity, only d U > 1 is considered. Figure 1 (lower right) shows the treelevel diagram considered in this paper for the production of unparticles associated with a Z boson.
The final SM extension considered in this paper is the Arkani-Hamed-Dimopoulos-Dvali (ADD) model of large extra dimensions [37,38], which is motivated by the disparity between the electroweak (EW) unification scale (M EW ∼ 100 GeV) and the Planck scale (M Pl ∼ 10 19 GeV). This model predicts graviton (G) production via the process qq → Z + G, as shown in Fig. 1 (lower right). The graviton escapes detection, leading to a mono-Z signature. In the ADD model, the apparent Planck scale in four spacetime dimensions is given by where M D is the fundamental Planck scale in the full (n+4)-dimensional spacetime and R is the compactification length scale of the extra dimensions. Assuming M D is of the same order as M EW , the observed large value of M Pl suggests values of R much larger than the Planck length. These values are on the order of nm for n = 3, decreasing with larger values of n. The consequence of the large compactification scale is that the mass spectrum of the Kaluza-Klein graviton states becomes nearly continuous [37,38], resulting in a broadened spectrum for the transverse momentum ( p T ) of the Z boson.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system [39]. The first level (L1), 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 (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [40].

Data samples and event reconstruction
This search uses pp collision events collected with the CMS detector during 2016, 2017, and 2018 corresponding to a total integrated luminosity of 137 fb −1 . The data sets from the three different years are analyzed independently with appropriate calibrations and corrections to take into account the different LHC running conditions and CMS detector performance.
Several SM processes can contribute to the mono-Z signature. The most important backgrounds come from diboson processes: WZ → ν where one lepton escapes detection, ZZ → νν, and WW → νν. There can also be contributions where energetic leptons are produced by decays of top quarks in tt or tW events. Smaller contributions may come from triple vector boson processes (WWZ, WZZ, and ZZZ), ttW → WWbbW, ttZ → WWbbZ, and ttγ → WWbbγ, referred to collectively as VVV due to the similar decay products. Drell-Yan (DY) production of lepton pairs, Z/γ * → , has no intrinsic source of missing transverse momentum but can still mimic a mono-Z signature when the momentum of the recoiling system is poorly measured. A minor source of background is from events with a vector boson and a misreconstructed photon, referred to as Vγ.
Monte Carlo simulated events are used to model the expected signal and background yields. Three sets of simulated events for each process are used in order to match the different data taking conditions. The samples for DM production are generated using the dmsimp package [41,42] interfaced with MadGraph5_amc@nlo 2.4.2 [43][44][45][46]. The pseudoscalar and scalar model samples are generated at leading order (LO) in quantum chromodynamics (QCD), while the vector and axial-vector model samples are generated at next-to-leading-order (NLO) in QCD. The powhegv2 [47][48][49][50][51] generator is used to simulate the Zh signal process of the invisible Higgs boson at NLO in QCD, as well as the tt, tW, and diboson processes. The BSM Higgs boson production cross sections, as a function of the Higgs boson mass for the Zh process are taken from Ref. [52]. Samples for the 2HDM+ a model are generated at NLO with MadGraph5_amc@nlo 2.6.0. Events for both the ADD and unparticle models are generated at LO using an EFT implementation in pythia 8.205 in 2016 and 8.230 in 2017 and 2018 [53,54]. In order to ensure the validity of the effective theory used in the ADD model, a truncation method, described in Ref. [55], is applied. Perturbative calculations are only valid in cases where the square of the center-of-mass energy (ŝ) of the incoming partons is smaller than the fundamental scale of the theory (M 2 D ). As such, this truncation method suppresses the cross section for events withŝ  [57] and the simplified DM (2HDM+ a) samples, which uses CP3 [58] (CP5) tunes for all years. All events are processed through a simulation of the CMS detector based on Geant4 [59] and are reconstructed with the same algorithms as used for data. Simultaneous pp collisions in the same or nearby bunch crossings, referred to as pileup, are also simulated. The distribution of the number of such interactions in the simulation is chosen to match the data, with periodic adjustments to take account of changes in LHC operating conditions [60]. The average number of pileup interactions was 23 for the 2016 data and 32 for the 2017 and 2018 data.
Information from all subdetectors is combined and used by the CMS particle-flow (PF) algorithm [61] for particle reconstruction and identification. The PF algorithm aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energies of photons are obtained from the ECAL measurement. The energies of electrons are determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy from the corresponding ECAL cluster, and the energy sum from all bremsstrahlung photons spatially compatible with originating from the electron track. The momentum of muons is obtained from the curvature of the corresponding track in the tracker detector in combination with information from the muon stations. The energies of charged hadrons are determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers. Finally, the energies of neutral hadrons are obtained from the corresponding corrected ECAL and HCAL energies.
The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [62,63] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
Both electron and muon candidates must pass certain identification criteria to be further selected in the analysis. They must satisfy requirements on the transverse momentum and pseudorapidity: p T > 10 GeV and |η| < 2.5 (2.4) for electrons (muons). At the final level, a medium working point [64,65] is chosen for the identification criteria, including requirements on the impact parameter of the candidates with respect to the primary vertex and their isolation with respect to other particles in the event. The efficiencies for these selections are about 85 and 90% for each electron and muon, respectively.
In the signal models considered in this paper, the amount of hadronic activity tends to be small, so events with multiple clustered jets are vetoed. For each event, hadronic jets are clustered from reconstructed particle candidates using the infrared and collinear safe anti-k T algorithm [62,63] with a distance parameter of 0.4. Jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, on average, within 5 to 10% of the true momentum over the entire spectrum and detector acceptance. Pileup interactions can contribute additional tracks and calorimetric energy depositions to the jet momentum. To mitigate this effect, charged particles identified to be originating from pileup vertices are discarded and an offset is applied to correct for remaining contributions [66]. Jet energy corrections are derived from simulation to bring the measured response of jets to the average of simulated jets clustered from the generated final-state particles. In situ measurements of the momentum balance in dijet, photon+jet, Z+jet, and multijet events are used to determine corrections for residual differences between jet energy scale in data and simulation [66]. The jet energy resolution amounts typically to 15% at 10 GeV, 8% at 100 GeV, and 4% at 1 TeV. Additional selection criteria are applied to each jet to remove jets potentially dominated by anomalous contributions from some subdetector components or reconstruction failures [67]. Jets with p T > 30 GeV and |η| < 4.7 are considered for the analysis.
To identify jets that originated from b quarks, we use the medium working point of the DeepCSV algorithm [68]. This selection was chosen to remove events from top quark decays originating specifically from tt production, without causing a significant loss of signal. For this working point, the effi-ciency to select b quark jets is about 70% and the probability for mistagging jets originating from the hadronization of gluons or u/d/s quarks is about 1% in simulated tt events.
To identify hadronically decaying τ leptons (τ h ), we use the hadron-plus-strips algorithm [69]. This algorithm constructs candidates seeded by PF jets that are consistent with either a single or triple charged pion decay of the τ lepton. In the single charged pion decay mode, the presence of neutral pions is detected by reconstructing their photonic decays. Mistagged jets originating from non-τ decays are rejected by a discriminator that takes into account the pileup contribution to the neutral component of the τ h decay [69]. The efficiency to select real hadronically decaying τ leptons is about 75% and the probability for mistagging jets is about 1%.
The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [70]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event. Events with anomalously high p miss T can originate from a variety of reconstruction failures, detector malfunctions, or noncollision backgrounds. Such events are rejected by event filters that are designed to identify more than 85-90% of the spurious highp miss T events with a misidentification rate of less than 0.1% [70].

Event selection
Events with electrons (muons) are collected using dielectron (dimuon) triggers, with thresholds of p T > 23 (17) GeV and p T > 12 (8) GeV for the electron (muon) with the highest and second-highest measured p T , respectively. Single-electron and single-muon triggers with p T thresholds of 25 (27) and 20 (24) GeV for 2016 (2017-2018) are used to recover residual inefficiencies, ensuring a trigger efficiency above 99% for events passing the offline selection.
In the signal region (SR), events are required to have two (N = 2) well-identified, isolated electrons or muons with the same flavor and opposite charge (e+e− or μ+μ−). At least one electron or muon of the pair must have p T > 25 GeV, while the second must have p T > 20 GeV. In order to reduce nonresonant background, the dilepton invariant mass is required to be within 15 GeV of the world-average Z boson mass m Z [71]. Additionally, we require the p T of the dilepton system p T to be larger than 60 GeV to reject the bulk of the DY background. Since little hadronic activity is expected for the signal, we reject events having more than one jet with p T > 30 GeV within |η| < 4.7. The top quark background is further suppressed by rejecting events containing any b-tagged jet with p T > 30 GeV reconstructed within the tracker acceptance of |η| < 2.4. To reduce the WZ background in which both bosons decay leptonically, we remove events containing additional electrons or muons with loose identification and with p T > 10 GeV. Events containing a loosely identified τ h candidate with p T > 18 GeV and |η| < 2.3 are also rejected. Decays that are consistent with production of muons or electrons are rejected by an overlap veto.
In addition to the above criteria, there are several selections designed to further reduce the SM background. The main discriminating variables are: the missing transverse momentum, p miss T ; the azimuthal angle formed between the dilepton p T and the p miss T , Δφ( p T , p miss T ); and the balance ratio, | p miss T − p T |/ p T . The latter two variables are especially powerful in rejecting DY and top quark processes. Selection criteria are optimized to obtain the best signal sensitivity for the range of DM processes considered. The final selection requirements are: p miss T > 100 GeV, Δφ( p T , p miss T ) > 2.6 radians, and | p miss For the 2HDM+ a model, the selection differs slightly. We make a less stringent requirement on the missing transverse momentum, p miss T > 80 GeV, and require the trans- ] to be greater than 200 GeV. The kinematic properties of the 2HDM+ a production yield a peak in the m T spectrum near the neutral Higgs scalar (H) mass that is advantageous for background discrimination.
In order to avoid biases in the p miss T calculation due to jet mismeasurement, events with one jet are required to have the azimuthal angle between this jet and the missing transverse momentum, Δφ( p j T , p miss T ), larger than 0.5 radians. To reduce the contribution from backgrounds such as WW and tt, we apply a requirement on the distance between the two leptons in the (η, φ) plane, ΔR < 1.8, where ΔR = √ (Δφ ) 2 + (Δη ) 2 . A summary of the selection criteria for the SR is given in Table 1.

Background estimation
We estimate the background contributions using combined information from simulation and control regions (CRs) in data. A simultaneous maximum likelihood fit to the p miss T or m T distributions in the SR and CRs constrains the background normalizations and their uncertainties. Specific CRs target different categories of background processes, as described below.

The three-lepton control region
The WZ → ν decay mode can contribute to the SR when the third lepton ( = e or μ) escapes detection, and this same process can be monitored in an orthogonal CR, where  the third lepton is identified and then removed. The construction of the three-lepton (3 ) CR is based on events with three well-reconstructed charged leptons. A Z boson candidate is selected in the same manner as for the SR , while an additional electron or muon with identical quality and isolation is required. In cases where there are multiple Z boson candidates, the candidate with invariant mass closest to the Z boson mass is selected. To enhance the purity of the WZ selection, p miss T of at least 30 GeV is required and the invariant mass of three leptons is required to be larger than 100 GeV. The backgrounds in this CR are similar to those in the SR, with a sizable nonprompt background from DY events where a jet is misidentified as a lepton [72]. An additional minor source of background is from events with a vector boson and a misreconstructed photon (Vγ). All background estimates for this CR are taken from simulation.
To simulate the consequences of not detecting the third lepton, the "emulated p miss T " is estimated from the vectorial sum of p miss T and the transverse momentum ( p T ) of the additional lepton. The emulated p miss T is then used in place of the reconstructed p miss T and the same selection is applied as for the SR. Since there is negligible contamination from WZ → τν and top quark backgrounds in this CR, no veto is applied on additional τ h or b jet candidates. The resulting emulated p miss T spectrum is shown in Fig. 2 (upper). For the 2HDM+ a case, the "emulated m T " is used instead of "emulated p miss T " with the same selections.

The four-lepton control region
The ZZ process contributes to the SR through the ZZ → νν decay mode, and the same production process can be monitored via the decay mode ZZ → 4 . The 4 CR is based on events with two pairs of charged leptons. Each pair comprises two leptons of opposite charge and the same flavor and corresponds to a Z candidate. Two of the four leptons must fulfill the same requirements on the leptons as in the SR, while, in order to increase the yield, the other two leptons need only pass relaxed lepton quality requirements. The highest p T Z boson candidate is required to have an invariant mass within 35 GeV of the Z boson mass m Z [71]. Additionally, we require the transverse momentum of this Z boson candidate to be larger than 60 GeV. Additional backgrounds to the ZZ final state are events from triboson processes, events with a vector boson and a higgs boson (Vh) and from nonprompt events. These backgrounds are almost negligible. All background estimates for this CR are taken from simulation. For these four-lepton events, the emulated p miss T is calculated as the vectorial sum of the p miss T and the p T of the Z boson candidate with the larger absolute mass difference to m Z . The choice of which Z boson to use as a proxy for an invisibly decaying boson negligibly alters the emulated p miss T spectrum. The same selection as the SR is then applied using the emulated p miss T in place of the reconstructed p miss T , with the exception of the τ h and b jet candidate vetoes. The resulting emulated p miss T spectrum is shown in Fig. 2 (lower). Similarly to the 3 CR, the "emulated m T " is used instead of "emulated p miss T " for the 2HDM+ a case and the distribution is well described by the SM background estimations.

The electron-muon control region
We estimate the contribution of the flavor-symmetric backgrounds from an eμ CR based on events with two leptons of different flavor and opposite charge (e ± μ ∓ ) that pass all other analysis selections. This CR is largely populated by nonresonant backgrounds (NRB) consisting mainly of leptonic W boson decays in tt, tW, and WW events, where the dilepton mass happens to fall inside the Z boson mass window. Small contributions from single top quark events produced via s-and t-channel processes, and Z → ττ events in which τ leptons decay into light leptons and neutrinos, are also considered in the NRB estimation.

The DY control region
The DY background is dominant in the region of low p miss T . This process does not produce undetectable particles. Therefore, any nonzero p miss T arises from mismeasurement or limitations in the detector acceptance. The estimation of this background uses simulated DY events, for which the normalization is taken from data in a sideband CR of 80 < p miss T < 100 GeV where the signal contamination is negligible, with all other selections applied. For the 2HDM+ a analysis, a similar approach is taken with relaxed p miss T selection of 50 < p miss T < 100 GeV and an additional selection of m T < 200 GeV applied. The sideband CR is included in the maximum likelihood fit and a 100% uncertainty is assigned to the extrapolation from this CR to the SR. This uncertainty has little effect on the results because of the smallness of the overall contribution from the DY process in the SR.

Fitting method
After applying the selection, we perform a binned maximum likelihood fit to discriminate between the potential signal and the remaining background processes. The data sets for each data-taking year are kept separate in the fit. This yields a better expected significance than combining them into a single set because the signal-to-background ratios are different for the three years due to the different data-taking conditions. The electron and muon channels have comparable signalto-background ratios, and are combined in the fit, while the contributions, corrections and systematic uncertainties are calculated individually.
The p miss T distribution of events passing the selection is used as the discriminating variable in the fit for all of the signal hypotheses except for the 2HDM+ a model. For this model, the m T distribution is used since a Jacobian peak around the pseudoscalar Higgs boson mass is expected. Events in the SR are split into 0-jet and 1-jet categories to take into account the different signal-to-background ratios. In addition, for the CRs defined in Sect. 6, events with 0-jet and 1-jet are included as a single category in the fit. The eμ and DY CRs are each included as a single bin corresponding to the total yield. The p miss T or m T spectra in the 3 and 4 CRs are included in the fit with the same binning as in the SR, where these spectra are based upon the emulated p miss T . To allow for further freedom in the ZZ and WZ background estimation, the p miss T and emulated p miss T distributions are split into three regions with independent normalization parameters: low (< 200 GeV), medium (200-400 GeV), and high (> 400 GeV), with uncertainties of 10, 20, and 30%, respectively. These values are based on the magnitudes of the theoretical uncertainties as described in Sect. 8. For fits to the 2HDM+ a model, three similar m T regions are cho-sen with the same uncertainties: low (< 400 GeV), medium (400-800 GeV), and high (> 800 GeV). To make the best use of the statistical power in the CRs and to take advantage of the similarities of the production processes, we take the normalization factors to be correlated for the ZZ and the WZ backgrounds in each p miss T region. For each individual bin, a Poisson likelihood term describes the fluctuation of the data around the expected central value, which is given by the sum of the contributions from signal and background processes. Systematic uncertainties are represented by nuisance parameters θ with log-normal probability density functions used for normalization uncertainties and Gaussian functions used for shape-based uncertainties, with the functions centered on their nominal valueŝ θ . The uncertainties affect the overall normalizations of the signal and background templates, as well as the shapes of the predictions across the distributions of observables. Correlations among systematic uncertainties in different categories are taken into account as discussed in Sect. 8. The total likelihood is defined as the product of the likelihoods of the individual bins and the probability density functions for the nuisance parameters: The factors of the likelihood can be written more explicitly as The purpose of the fit is to determine the confidence interval for the signal strengths μ. Here P(N | λ) is the Poisson probability to observe N events for an expected value of λ, and f NP (θ |θ) describes the nuisance parameters with lognormal probability density functions used for normalization uncertainties and Gaussian functions used for shape-based uncertainties. The index i indicates the bin of the p miss T or m T distribution, r (i) corresponds to the region (low, medium, high) of bin i, and the index j indicates either the 0-jet or 1-jet selection. The diboson process normalization in the region r (i) is μ VV,r (i) , while μ DY is the DY background normalization and μ NRB is the normalization for the nonresonant background. The yield prediction from simulation for process x in region y is noted as N y x . The smaller backgrounds in each region are merged together and are indicated collectively as "other". The method above for constructing likelihood functions follows that of Ref. [73], where a more detailed mathematical description may be found.

Systematic uncertainties
In the following, we describe all of the uncertainties that are taken into account in the maximum likelihood fit. We consider the systematic effects on both the overall normalization and on the shape of the distribution of p miss T or m T for all applicable uncertainties. We evaluate the impacts by performing the full analysis with the value of the relevant parameters shifted up and down by one standard deviation. The final varied distributions of p miss T or m T are used for signal extraction and as input to the fit. For each source of uncertainty, variations in the distributions are thus treated as fully correlated, while independent sources of uncertainty are treated as uncorrelated. Except where noted otherwise, the systematic uncertainties for the three different years of data taking are treated as correlated.
The assigned uncertainties in the integrated luminosity are 2.5, 2.3, and 2.5% for the 2016, 2017, and 2018 data samples [74-76], respectively, and are treated as uncorrelated across the different years.
We apply scale factors to all simulated samples to correct for discrepancies in the lepton reconstruction and identification efficiencies between data and simulation. These factors are measured using DY events in the Z boson peak region [65,77,78] that are recorded with unbiased triggers. The factors depend on the lepton p T and η and are within a few percent of unity for electrons and muons. The uncertainty in the determination of the trigger efficiency leads to an uncertainty smaller than 1% in the expected signal yield.
For the kinematic regions used in this analysis, the lepton momentum scale uncertainty for both electrons and muons is well represented by a constant value of 0.5%. The uncertainty in the calibration of the jet energy scale (JES) and resolution directly affects the p miss T computation and all the selection requirements related to jets. The estimate of the JES uncertainty is performed by varying the JES. The variation corresponds to a re-scaling of the jet four-momentum as p → p(1 ± δp JES T / p T ), where δp JES T is the absolute uncer-tainty in the JES, which is parameterized as function of the p T and η of the jet. In order to account for the systematic uncertainty from the jet resolution smearing procedure, the resolution scale factors are varied within their uncertainties. Since the uncertainties in the JES are derived independently for the three data sets, they are treated as uncorrelated across the three data sets. The signal processes are expected to produce very few events containing b jets, and we reject events with any jets that satisfy the b tagging algorithm working point used. In order to account for the b tagging efficiencies observed in data, an event-by-event reweighting using b tagging scale factors and efficiencies is applied to simulated events. The uncertainty is obtained by varying the event-by-event weight by ±1 standard deviation. Since the uncertainties in the b tagging are derived independently for the three data sets, they are treated as uncorrelated across the three data sets. The variation of the final yields induced by this procedure is less than 1%.
Simulated samples are reweighted to reproduce the pileup conditions observed in data. We evaluate the uncertainty related to pileup by recalculating these weights for variations in the total inelastic cross section by 5% around the nominal value [79]. The resulting shift in weights is propagated through the analysis and the corresponding p miss T and m T spectra are used as input to the maximum likelihood fit. The variation of the final yields induced by this procedure is less than 1%.
Shape-based uncertainties for the ZZ and WZ backgrounds, referred to jointly as VV, and signal processes are derived from variations of the renormalization and factorization scales, the strong coupling constant α S , and PDFs [80][81][82]. The scales are varied up and down by a factor of two. Variations of the PDF set and α S are used to estimate the corresponding uncertainties in the yields of the signal and background processes following Ref. [56]. The missing higherorder EW terms in the event generation for the VV processes yield another source of theoretical uncertainty [83,84]. The following additional higher-order corrections are applied: a constant (approximately 10%) correction for the WZ cross section from NLO to NNLO in QCD calculations [85]; a constant (approximately 3%) correction for the WZ cross section from LO to NLO in EW calculations, according to Ref. [86]; a Δφ(Z, Z)-dependent correction to the ZZ production cross section from NLO to next-to-next-to-leading order (NNLO) in QCD calculations [87]; a p T -dependent correction to the ZZ cross section from LO to NLO in EW calculations, following Refs. [83,84,86], which is the dominant correction in the signal region. We use the product of the above NLO EW corrections and the inclusive NLO QCD corrections [88] as an estimate of the missing NLO EW×NLO QCD contribution, which is not used as a correction, but rather assigned as an uncertainty. The resulting variations in the p miss T and m T and m T distributions are needed for each of the background processes. For the DY and nonresonant processes, we take the shape directly from simulation. The distributions for the ZZ and WZ processes are obtained by taking the shapes from the simulation and normalizing them to the yield seen in the data in the CR. The gluon-induced and the quark-induced ZZ processes have different acceptances and their uncertainties are treated separately, while the normalization factors are taken to be correlated. In all cases, the limited number of simulated events in any given bin gives rise to a systematic uncertainty. This uncertainty is treated as fully uncorrelated across the bins and processes.
A summary of the impact on the signal strength of the systematic uncertainties is shown in Table 2. The Zh(invisible) model is used as an example to illustrate the size of the uncertainties, both for the presence (B(h → invisible) = 1) and absence (B(h → invisible) = 0) of a signal. These two paradigms are used to generate Asimov data sets that are then fit to give the uncertainty estimates shown in Table 2. The systematic uncertainties are dominated by the theoretical uncertainty in the ZZ and WZ background contributions.

Results
The number of observed and expected events in the SR after the final selection is given in Table 3, where the values of the  expected yields and their uncertainties are obtained from the maximum likelihood fit. The observed numbers of events are compatible with the background predictions. The expected yields and the product of acceptance and efficiency for several signal models used in the analysis are shown in Table 4. The post-fit p miss T distributions for events in the signal region in the 0-jet and 1-jet categories are shown in Fig. 3. The final m T distributions used for the 2HDM+ a model are shown in Fig. 4.
For each of the models considered, simulated signal samples are generated for relevant sets of model parameters. The observed p miss T and m T spectra are used to set limits on theories of new physics using the modified frequentist construction CL s [73, 89,90] used in the asymptotic approximation [91].

Simplified dark matter model interpretation
In the framework of the simplified models of DM, the signal production is sensitive to the mass, spin, and parity of the tor mediators, we observe a limit around m med > 870 GeV for most values of m χ less than m med /2. For axial-vector mediators the highest limit reached in the allowed region is about m med > 800 GeV. In both cases, the previous limits from this channel are extended by about 150 GeV, but the limits are still less restrictive than those from published mono-jet results [7] because weakly coupled Z bosons are radiated from the initial state quarks much less frequently than gluons. Figure 6 shows the 90% CL limits on the DM- Fig. 5 The 95% CL exclusion limits for the vector (upper) and the axial-vector (lower) simplified models. The limits are shown as a function of the mediator and DM particle masses. The coupling to quarks is fixed to g q = 0.25 and the coupling to DM is set to g χ = 1 nucleon cross sections calculated following the suggestions in Ref. [22]. Limits are shown as a function of the DM particle mass for both the spin-independent and spin-dependent cases and compared to selected results from direct-detection experiments.
In addition to vector and axial-vector mediators, scalar and pseudoscalar mediators are also tested. For these models, we fix both couplings to quarks and to DM particles: g q = 1 and g χ = 1 as suggested in Ref. [22]. Since the choice of DM particle mass is shown to have negligible effects on the kinematic distributions of the detected particles, we set it to the constant value of m χ = 1 GeV. Figure 7 gives the 95% CL exclusion limits on the production cross section over the predicted cross section as a Fig. 6 The 90% CL DM-nucleon upper limits on the cross section for simplified DM in the spin-independent (upper) and spin-dependent (lower) cases. The coupling to quarks is set to g q = 0.25 and the coupling to DM is set to g χ = 1. Limits from the XENON1T [93], LUX [94], PandaX-ll [95], CRESST-III [96], and DarkSide-50 [97] experiments are shown for the spin-independent case with vector couplings. Limits from the PICO-60 [98], PICO-2L [99], IceCube [100], and Super-Kamiokande [101] experiments are shown for the spindependent case with axial-vector couplings function the mediator mass m med . The expected limits are about 25% better than the previous results in this channel [4], but are not yet sensitive enough to exclude any value of m med . The best limits obtained on the cross section are about 1.5 times larger than the predicted values for low values of m med . Fig. 7 The 95% CL upper limits on the cross section for simplified DM models with scalar (upper) and pseudoscalar (lower) mediators. The coupling to quarks is set to g q = 1, the coupling to DM is set to g χ = 1 and the DM mass is m χ = 1 GeV

Two-Higgs-doublet model interpretation
For the 2HDM+ a model, the signal production is sensitive to the heavy Higgs boson and the pseudoscalar a masses. As discussed in Sect. 7, the m T distribution is used in the fit rather than p miss T . The limits on both the heavy Higgs boson and the additional pseudoscalar mediator a are shown in Fig. 8. The mixing angles are set to tan β = 1 and sin θ = 0.35 with a DM particle mass of m χ = 10 GeV. The mediator mass with the most sensitivity is m H = 1000 GeV, where the observed (expected) limit on m a is 440 (340) GeV. For small values

Invisible Higgs boson interpretation
For the search for invisible decays of the Higgs boson, we use the p miss T distribution as input to the fit. We obtain upper limits on the product of the Higgs boson production cross section and branching fraction to invisible particles σ Zh B(h → invisible). This can be interpreted as an upper limit on B(h → invisible) by assuming the production rate [52,103,104] for an SM Higgs boson at m h = 125 GeV. The observed (expected) 95% CL upper limit at m h = 125 GeV on B(h → invisible) is 29% (25 +9 −7 %) as shown in Fig. 9. The observed (expected) limit from the previous CMS result in this channel was B(h → invisible) < 45 (44) Table 5 gives a summary of the limits expected and observed for a selection of relevant parameters in all of the models considered.

Summary
Events with a Z boson recoiling against missing transverse momentum in proton-proton collisions at the LHC are used Table 5 Observed and expected 95% CL limits on parameters for the simplified DM models, invisible decays of the Higgs boson, two-Higgsdoublet model, large extra dimensions in the ADD scenario, and unparticle model. For the scalar and pseudoscalar mediators, the limits are dependent on the mediator mass, so the lowest values for the ratio of observed to theoretical cross sections are presented. For the vector and axial-vector mediators, the limits are dependent on the DM particle mass, so the limits are shown for m χ < 300 GeV for the vector mediator and m χ = 240 GeV for the axial-vector mediator

Model
Parameter  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .