Search for dark matter and unparticles in events with a Z boson and missing transverse momentum in proton-proton collisions at sqrt(s) = 13 TeV

A search for dark matter and unparticle production at the LHC has been performed using events containing two charged leptons (electrons or muons), consistent with the decay of a Z boson, and large missing transverse momentum. This study is based on data collected with the CMS detector in 2015, corresponding to an integrated luminosity of 2.3 inverse femtobarns of proton-proton collisions at the LHC, at a center-of-mass energy of 13 TeV. No excess over the standard model expectation is observed. Compared to previous searches in this topology, which exclusively relied on effective field theories, the results are interpreted in terms of a simplified model of dark matter production for both vector and axial vector couplings between a mediator and dark matter particles. The first study of this class of models using CMS data at sqrt(s) = 13 TeV is presented. Additionally, effective field theories of dark matter and unparticle production are used to interpret the data.


Introduction
According to the well-established Λ CDM model of cosmology, known matter only comprises about 5% of the total energy content of the universe, with 27% contributed by dark matter (DM) and the rest being dark energy [1]. Although strong astrophysical evidence indicates the existence of DM, there is no evidence yet for nongravitational interactions between DM and standard model (SM) particles. DM searches exploit a number of methods including direct detection [2] and indirect detection [3]. If there are DM particles that can be observed in direct detection experiments, they could have substantial couplings to nucleons, and therefore could be produced at the CERN LHC. A theoretically promising possibility is that DM may take the form of weakly interacting massive particles. Searches for production of such particles at colliders typically consider the case of DM recoiling against a standard model particle ("tag") to obtain a defined signature [4]. Such searches have been performed using various standard model signatures as tags [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In models where DM production is mediated by an interaction involving SM quarks, the monojet signature is typically the most sensitive. If DM particles are instead produced via radiation emitted by a standard model boson, searches in the Z/W/γ + E miss T channels are advantageous.
The study presented here considers the case of a Z boson recoiling against a pair of DM particles, χχ. The Z boson subsequently decays into two charged leptons ( + − , where = e or µ) producing a well-defined signature together with missing transverse momentum due to the undetected DM particles. A simplified tree-level ultraviolet-complete model [4] that contains a massive spin-1 mediator exchanged in the s-channel is considered here. In this model, the spin-1 mediator A could have either vector or axial-vector couplings to the SM and DM particles. The DM particle χ is assumed to be a Dirac fermion. The interaction Lagrangian of the s-channel vector mediated DM model can be written as: where the mediator is labeled as A, and its coupling to DM particles is labeled as g χ . The coupling between the mediator and SM quarks is labeled as g q , and is assumed to be universal to all quarks. The Lagrangian for an axial-vector mediator is obtained by making the replacement γ µ → γ µ γ 5 in all terms. As a benchmark model for DM production via a scalar coupling, an effective field theory (EFT) with dimension-7 operators is also considered [4]. It contains SU(2) L × U(1) Y gauge invariant couplings between a DM pair and two SM gauge bosons in a four-particle contact interaction. The corresponding interaction Lagrangian is: in which B µν and F i µν are the U(1) Y and SU(2) L field tensors, and Λ denotes the cutoff scale. The coupling parameter c 1 controls the relative importance of the U(1) Y and SU(2) L fields for DM production. Any multiplicative factor for the U(1) Y and SU(2) L couplings is absorbed into Λ. Note that the choice of Λ modifies the signal cross section, but not the expected kinematic properties of events. The model is nonrenormalizable and should be considered as a benchmark of the sensitivity to this class of interaction. It should be used with caution when making comparisons with other sources of DM constraints, such as direct detection experiments. Figure 1 shows the Feynman diagrams for production of DM pairs (χχ) in association with a Z boson in these two types of models. The signature for DM production considered in this paper is the production of a pair of leptons (e + e − or µ + µ − ) consistent with a Z boson decay, together with a large missing transverse momentum. This same signature is sensitive to other models of physics beyond the SM (BSM), e.g. "unparticles"(U).
The unparticle physics concept [21-24] is particularly interesting because it is based on scale invariance, which is anticipated in many BSM physics scenarios [25][26][27]. The effects of the scale invariant sector (unparticles) appear as a noninteger number of invisible massless particles. In this scenario, the SM is extended by introducing a scale invariant Banks-Zaks (BZ) field, which has a nontrivial infrared fixed point [28]. This field can interact with SM particles by exchanging heavy particles with a high mass scale M U . Below this mass scale, the coupling is nonrenormalizable and the interaction is suppressed by powers of M U . The EFT Lagrangian can be expressed as: in which C U is a normalization factor, d U represents the possible noninteger scaling dimension of the unparticle operator O U , O SM is an operator composed of SM fields with dimension d SM , k = d SM + d BZ − 4 > 0 is the scaling dimension, Λ U is the energy scale of the interaction, and d BZ denotes the scaling dimension of the BZ operator at energy scales above Λ U . The parameter λ = C U Λ d BZ U /M k U is a measure of the coupling between SM particles and unparticles. The scaling dimension d U ≥ 1 is constrained by the unitarity condition. Additional details regarding this unparticle model are available in Ref. [17].
In this paper, real emission of scalar unparticles is considered. The unparticles are assumed to couple to the standard model quarks in an effective three-particle interaction. In the scalar unparticle case, O SM = qq, which yields numerically identical results to the pseudo-scalar operator choice O SM = qiγ 5 q [29]. Figure 2 shows the corresponding tree-level diagram for the production of unparticles associated with a Z boson.
The analysis is based on a data set recorded with the CMS detector in 2015 in pp collisions at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 2.3 ± 0.1 fb −1 . A previous CMS search in the same final state [17], based on data collected at a center-ofmass energy of 8 TeV, found no evidence of new physics and set limits on DM and unparticle production using an EFT description. A CMS analysis of the 8 TeV data set in the combined monojet and hadronic mono-V (where V = W or Z) channels [19] has previously set limits on the simplified model parameters considered here. Dark matter particle masses of up to 500 q q − + U Z Figure 2: Leading order Feynman diagram for unparticle (denoted by U) production in association with a Z boson. The hatched circle indicates the interaction modeled with an EFT operator.
GeV (400 GeV) and mediator masses of up to 1.6 TeV have been excluded in the vector (axialvector) coupling scenarios for g q = g DM = 1. A search performed by the ATLAS Collaboration using √ s = 13 TeV data corresponding to an integrated luminosity of 3.2 fb −1 in events with a hadronically decaying V boson and E miss T has recently reported exclusion of the dimension-7 EFT scenario up to Λ = 700 GeV (460 GeV) for DM particle masses of 1 GeV (1 TeV) [20].

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 (η) [30] coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. 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. [30].
Variables of particular relevance to the present analysis are the missing transverse momentum vector p miss T and the magnitude of this quantity, E miss T . The quantity p miss T is defined as the projection on the plane perpendicular to the beams of the negative vector sum of the momenta of all reconstructed particles in an event.

Simulation
Samples of simulated DM particle events for both the simplified model and EFT interpretations are generated using MADGRAPH5 aMC@NLO 2.2.2 [31] at leading order (LO) and matched to PYTHIA 8.205 [32] using tune CUETP8M1 for parton showering and hadronization [33,34]. The factorization and renormalization scales are set to the geometric mean of √ p 2 T + m 2 for all final-state particles [4,31], where p T and m are the transverse momentum and mass of each particle.
For the simplified model of DM production, couplings are chosen according to the recommendations in Ref. [35]. The coupling g χ is set to unity. For g q , values of 1.0 and 0.25 are considered. The width of the mediator is assumed to be determined exclusively by the contributions from the couplings to quarks and the DM particle χ. Under this assumption, the width is in the range 1-5% (30-50%) of the mediator mass for g q = 0.25 (g q = 1.00). The signal simulation samples with g q = 1.0 are processed using the detector simulation described below. Signal predictions for g q = 0.25 are obtained by applying event weights based on the E miss T distribution at the generator level to the fully simulated samples with g q = 1.0. This procedure allows to take into account any effect of the coupling dependent mediator width on the E miss T distribution [35]. The exact dependence of the width on the model parameters is reported in [35].
Samples for the EFT DM benchmark are generated with Λ = 3 TeV and c 1 = 1. Signal predictions for other values of Λ are obtained by rescaling the signal cross section accordingly, while other values of c 1 are evaluated using the same reweighting method as for the simplified model case.
The events for the unparticle model are generated at LO with PYTHIA 8 [29,36] assuming a cutoff scale Λ U = 15 TeV, using tune CUETP8M1 for parton showering and hadronization. We evaluate other values of Λ U by rescaling the cross sections as needed. The parameter Λ U acts solely as a scaling factor for the cross section and does not influence the kinematic distributions of unparticle production [29].
The POWHEG 2.0 [37][38][39][40][41] event generator is used to produce samples of events for the tt, tW, qq → ZZ, and WZ background processes, which are simulated at next-to-leading order (NLO). The gg → ZZ process is simulated using MCFM 7.0.1 [42] at NLO. The Drell-Yan (DY, Z/γ * → + − ) process is generated using the MADGRAPH5 aMC@NLO event generator at LO and normalized to the next-to-next-to-leading order (NNLO) cross section as calculated using FEWZ 3.1 [43,44]. Triboson events (WZZ, WWZ and ZZZ) are simulated using MAD-GRAPH5 aMC@NLO at NLO. Samples of quantum chromodynamics (QCD) production of multijet events are generated using PYTHIA 8 at LO. For all SM simulation samples, parton showering and hadronization are performed with PYTHIA 8 with tune CUETP8M1.
The parton distribution function (PDF) set NNPDF3.0 [45] is used for Monte Carlo (MC) samples, and the detector response is simulated using a detailed description of the CMS detector, based on the GEANT4 package [46,47]. Minimum bias events are superimposed on the simulated events to emulate the effect of additional pp interactions in the same or nearby bunch crossings (pileup). All MC samples are corrected to reproduce the pileup distribution as measured in the data. The average number of pileup interactions per proton bunch crossing is about 12 for the 2015 data sample.
The upper left panel of Fig. 3 shows the distribution of E miss  event vertex. Simulation studies show that this requirement correctly selects the event vertex in more than 99% of both signal and background events. The lepton candidate tracks are required to be compatible with the event vertex.
A particle-flow (PF) event algorithm [49,50] reconstructs and identifies each individual particle with an optimized combination of information from the various elements of the CMS detector. Photon energies are directly obtained from the ECAL measurement, corrected for zerosuppression effects [30]. Electron energies are determined from a combination of the electron momentum at the event vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. Muon momenta are obtained from the curvature of the corresponding track. Charged hadron energies are determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, cor-rected for zero-suppression effects and for the response function of the calorimeters to hadronic showers [50]. Finally, neutral hadron energies are obtained from the corresponding corrected ECAL and HCAL energies.
Electron candidates are reconstructed using an algorithm that combines information from the ECAL and the tracker [51]. To reduce the electron misidentification rate, the candidates have to satisfy additional identification criteria that are based on the shape of the electromagnetic shower in the ECAL. In addition, the electron track is required to originate from the event vertex and to match the shower cluster in the ECAL. Electron candidates with an ECAL cluster in the transition region between ECAL barrel and endcap (1.44 < |η| < 1.57) are rejected because the reconstruction of an electron candidate in this region is not optimal. Candidates that are identified as coming from photon conversions [51] in the detector material are explicitly removed.
Muon candidate reconstruction is based on two algorithms: in the first, tracks in the silicon tracker are matched with at least one muon segment in any detector plane of the muon system, and in the second algorithm, a combined fit is performed to hits in both the silicon tracker and the muon system [52]. The muon candidates in this analysis are required to be reconstructed with at least one of the two algorithms and to be further identified as muons by the PF algorithm. To reduce the muon misidentification rate, additional identification criteria are applied based on the number of spatial points measured in the tracker and in the muon system, the fit quality of the muon track, and its consistency with the event vertex location.
Leptons produced in the decay of Z bosons are expected to be isolated from hadronic activity in the event. Therefore, an isolation requirement is applied based on the sum of the momenta of the PF candidates found in a cone of radius ∆R = 0.4 around each lepton. The isolation sum is required to be smaller than 15% (20%) of the p T of the electron (muon). For each electron, the mean energy deposit in the isolation cone of the electron, coming from other pp collisions in the same bunch crossing, is estimated following the method described in Ref. [51], and subtracted from the isolation sum. For muon candidates, only charged tracks associated with the event vertex are included. The sum of the p T for charged particles not associated with the event vertex in the cone of interest is rescaled by a factor of 0.5, corresponding to the average neutral to charged energy density ratio in jets, and subtracted from the isolation sum.
For the purpose of rejecting events containing τ leptons, hadronically decaying τ leptons (τ h ) are identified using the "hadron-plus-strips" algorithm. The algorithm identifies a jet as a τ h candidate if a subset of the particles assigned to the jet is consistent with the decay products of a τ h [53]. In addition, τ h candidates are required to be isolated from other activity in the event.
Jets are reconstructed from PF candidates by using the anti-k T clustering algorithm [54] with a distance parameter of 0.4, as implemented in the FASTJET package [55,56]. Jets are identified over the full calorimeter acceptance, |η| < 5. The jet momentum is defined as the vector sum of all particle momenta assigned to the jet, and is found in simulation to be within 5 to 10% of the true hadron-level momentum over the whole p T range and detector acceptance. An overall energy subtraction is applied to correct for the extra energy clustered in jets due to pileup, following the procedure described in Ref. [57]. Additional corrections to the jet energy scale and resolution are derived from simulation, and are complemented by measurements of the energy balance in dijet and γ+jets events [57].

Event selection
A preselection with a large yield is used to validate the background model and is followed by a final selection that is designed to give maximal sensitivity to the signal, as quantified by the expected limits achieved. Preselected events are required to have exactly two wellidentified, isolated leptons with the same flavor and opposite charge (e + e − or µ + µ − ), each with p T > 20 GeV. The invariant mass of the lepton pair is required to be within ±10 GeV of the nominal mass of the Z boson [58]. Only electrons (muons) within the range of |η| < 2.5 (2.4) are considered. To reduce the background from the WZ process where the W boson decays leptonically, events are removed if an additional electron or muon is reconstructed with p T > 10 GeV. The event is also removed from the final selection if a τ h candidate is reconstructed with p T > 20 GeV. As a loose preselection requirement, the dilepton transverse momentum (p T ) is required to be larger than 50 GeV to reject the bulk of DY background events.
Since only a small amount of hadronic activity is expected in the final state of both DM and unparticle events, any event having two or more jets with p T > 30 GeV is rejected. Processes involving top quarks are further suppressed with the use of techniques based on soft-muon and secondary-vertex b jet tagging, aimed at identifying the b quarks produced in top quark decays. Soft muons are identified using a specialised low-p T set of identification criteria focused on the muon candidate track quality. The rejection of events with soft muons having p T > 3 GeV reduces the background from semileptonic decays of B mesons. The b jet tagging technique employed is based on the "combined secondary vertex" algorithm [59,60]. The algorithm is calibrated to provide, on average, 80% efficiency for tagging jets originating from b quarks, and 10% probability of light-flavor jet misidentification. Events are rejected if at least one b-tagged jet is reconstructed with p T > 20 GeV within the tracker acceptance (|η| < 2.5).
For the final selection, further kinematic requirements are set in order to achieve the best possible signal extraction. A minimal E miss T of 80 GeV is required. The angle between the Z boson and the missing transverse momentum in the transverse plane ∆φ , p miss T is required to be larger than 2.7 radians. The momentum balance of the event defined by |E miss T − p T |/p T is required to be smaller than 0.2. These variables suppress background processes such as DY and top quark production. The event selection criteria used for the electron and muon channels are the same. They are summarized in Table 1.    after preselection for the Z → e + e − (left) and Z → µ + µ − (right) channels. Representative expected signal distributions are shown for the simplified model of DM production with vector couplings, the EFT scenario of DM production, and unparticles. The SM expectation is based on simulation only. The total statistical uncertainty in the overall background prediction is shown as a hatched region. Overflow events are included in the rightmost bins. The upper error bars on data points are shown for bins with zero entries (Garwood procedure) in the region up to the last non-zero entry. In the lower panels, the ratio between data and predicted background is shown.

Background estimation
The ZZ and WZ backgrounds are modeled using MC simulation, and normalized to their respective NLO cross sections. Other backgrounds, including tt, tW, WW, Z → ττ, single top quark, and DY production are estimated from data for the final selection.
The simulation of the ZZ process includes the qq-and gg-induced production modes. In order to correct the ZZ differential cross section from NLO to NNLO in QCD, ∆φ(Z, Z)-dependent K-factors are applied [61]. We apply NLO electroweak (EW) K-factors as a function of the p T of the trailing boson, following the calculations in Refs. [62][63][64]. Electroweak corrections to WZ production are also available, but considered small [64] and not applied.
The background processes involving ee or µµ pairs not directly resulting from the decay of a Z boson are referred to as nonresonant backgrounds. These backgrounds arise mainly from leptonic W boson decays in tt, tW, and WW events. There are also small contributions from the sand t-channel single top quark events, W+jets events, and Z → ττ events in which τ lepton decays result in electrons or muons and E miss T . We estimate these nonresonant backgrounds using a data control sample, consisting of events with opposite-charge different-flavor dilepton pairs (e ± µ ∓ ) that otherwise pass the full selection. As the decay rates for Z → e + e − and Z → µ + µ − are almost equal, by equating the ratio of observed dilepton counts to the square of the ratio of efficiencies, the nonresonant backgrounds in the ee and µµ channels can be estimated from the eµ channel: in which the coefficient of 1/2 in the transfer factors k ee and k µµ comes from the dilepton decay ratios for ee, µµ, and eµ in these nonresonant backgrounds, and N data ee and N data µµ are the numbers of selected ee and µµ events from data with masses in the Z boson mass window. The ratio √ N data ee /N data µµ and the reciprocal quantity take into account the difference between the electron and muon selection efficiencies. The term N data, corr eµ is the number of eµ events observed in data corrected by subtracting the estimated ZZ, WZ, and DY background contributions. The kinematic distributions of the estimated nonresonant backgrounds are obtained from simulation with the overall normalization determined by the method described above. The validity of this procedure for predicting nonresonant backgrounds is checked with simulated events containing tt, tW, WW, W+jets, and Z → ττ processes. We assign a systematic uncertainty of 26% for this background estimation in both the electron and muon channels for E miss T > 80 GeV, based on closure tests that compare the predictions obtained from the control sample with those from the simulated events.
The DY process is dominant in the region of low E miss T . This process does not produce undetectable particles, and therefore the measured E miss T arises from limited detector acceptance and mismeasurement of particle momenta. The estimation of this background uses simulated DY events, which are normalized to data with scale factors obtained by measuring the number of DY events in a background-dominated control region, after subtracting other processes. These scale factors are of order 1.0-1.2. The control region is defined by applying the full selection with the E miss T requirement inverted. The reliability of this approach in the high-E miss T regime has been studied by considering variables sensitive to E miss T mismeasurement, such as the angular separation between the E miss T direction and any jet. A normalization uncertainty of 100%, which accommodates any differences observed in these control regions, is assigned for the DY background estimate. The assigned uncertainty has little impact on the overall signal sensitivity because of the small overall contribution from the DY background prediction.
Contributions from QCD production of multijet events is estimated using simulation and found to be negligible after final selection.

Efficiencies and systematic uncertainties
The efficiencies for selecting, reconstructing, and identifying isolated leptons are determined from simulation, and corrected with scale factors determined from applying a "tag-and-probe" technique [65] to Z → + − events in data. The trigger efficiencies for the electron and muon channels are found to be above 90%, varying as a function of p T and η of the lepton. The identification efficiency, when applying the selection criteria described in Section 4, is found to be about 80-86% for electrons and 95% for muons, depending on the p T and η of the corresponding lepton. The corresponding data-to-MC scale factors are typically in the range 0.96-1.00 for the electron and 0.96-0.98 for the muon channel, depending on the p T and |η| of the lepton candidate. The lepton momentum scale uncertainty is computed by varying the momentum of the leptons by its uncertainties. The lepton momentum uncertainty is 1% for the muons, while the uncertainty for the electrons is 2% in the barrel and 5% in the endcaps. For both channels, the overall uncertainty in the efficiency of selecting and reconstructing leptons in an event is about 3%.
In the treatment of systematic uncertainties, both normalization effects, which only affect the overall size of individual contributions, as well as shape uncertainties, which also affect their distribution, are taken into account. The systematic uncertainties are summarized in Table 2. Where applicable, the symbol V is used to refer to both Z and W bosons. The impact of each source of uncertainty on the observed strength of a potential signal is also reported. The signal strength is defined as the ratio of the observed or excluded signal cross-section to the signal cross-section predicted by theory. To calculate the impact, a maximum likelihood fit of the combined background and signal model to the expected distribution for unity signal strength is performed. The fit is repeated with each individual nuisance parameter varied by its uncertainty. The impact of the uncertainty is then defined as the relative change induced in the expected best fit signal strength by the variation of the respective parameter. In the table, the reference signal is the simplified model DM scenario with a vector mediator of mass 200 GeV, a DM particle mass of 50 GeV, and coupling g q = 1.0. Table 2: Summary of systematic uncertainties. Each background uncertainty represents the variation of the relative yields of the particular background components. The signal uncertainties represent the relative variations in the signal acceptance, and the ranges quoted cover both signals of DM and unparticles with different DM masses or scaling dimensions. For shape uncertainties, the numbers correspond to the overall effect of the shape variation on the yield or acceptance. The symbol "-" indicates that the systematic uncertainty is not applicable. The impact of each group of systematic uncertainties is calculated by performing a maximum likelihood fit to obtain the signal strength with each parameter separately varied by its uncertainty. The number given in the impact column is the relative change of the expected best fit signal strength that is introduced by the variation for the simplified model signal scenario with a vector mediator of mass 200 GeV, DM of mass 50 GeV, and coupling g q = 1.0. The normalization uncertainties in the background estimates from data have been described in Section 6. The PDF and α S uncertainties (referred to as PDF+α S in the following) for signal and background processes are estimated from the standard deviation (s.d.) of weights according to the replicas provided in the NNPDF3.0 parton distribution set [66]. While the influence on the estimated signal acceptance arising from theory-related uncertainties is included in the limit calculation, the corresponding effect on the normalization of the signal process is not. For the simplified model of DM production, the effect of the signal normalization uncertainty is treated separately from the experimental uncertainty and is shown as a dashed band around the observed limit. Since the EFT benchmark and unparticle scenarios are extremely simplified, theory-related cross-section uncertainties are not considered to be realistic for these models and are thus neglected. The efficiencies for signal, ZZ, and WZ processes are estimated using simulation, and the uncertainties in the corresponding yields are derived by varying the renormalization and factorization scales, α S , and choice of PDFs. The factorization and renormalization scale uncertainties are assessed by varying the original scales by factors of 0.5 or 2.0, and amount to 2-3% for ZZ and WZ processes. The effect of variations in α S and choice of PDFs is 2% for the ZZ and WZ backgrounds. A 3% normalization uncertainty is assigned to the WZ background to account for higher-order EW corrections [64]. The uncertainty assigned to the integrated luminosity measurement is 2.7% [67].

Source of uncertainty
Experimental sources of shape uncertainty are the lepton momentum scale, the jet energy scale and resolution, the b tagging efficiency, and the pileup modeling. The effect of each uncertainty is estimated by varying the respective variable of interest by its uncertainties, and propagating the variations to the distribution of E miss T after the final selection. In the case of the lepton momentum scale, the uncertainty is computed by varying the momentum of the leptons by their uncertainties. The uncertainty due to the lepton momentum scale is evaluated to be less than 1% (1-7%) for signal (background).
The uncertainties in the calibration of the jet energy scale and resolution directly affect the assignments of jets to jet categories, the E miss T computation, and all the selections related to jets. The effect of the jet energy scale uncertainty is estimated by varying the energy scale by ±1 s.d. A similar strategy is used to evaluate the systematic uncertainty related to the jet energy resolution. The effect of the shifts is propagated to E miss T . The uncertainties in the final yields are found to be less than 1% for signal and less than 4% for background.
In order to reproduce b tagging efficiencies observed in data, an event-by-event reweighting using data-to-simulation scale factors is applied to simulated events. The uncertainty associated with this procedure is obtained by varying the event-by-event weight by ±1 s.d. The total uncertainty in the final yields due to b tagging is less than 1% for both signal and background. All simulated events are reweighted to reproduce the pileup conditions observed in data. To compute the uncertainty related to pileup modeling, we shift the mean of the distribution in simulation by 5% [68]. The variation of the final yields induced by this procedure is 0.5-1% for signal and 1-2% for background. For the processes estimated from simulation, the sizes of the MC samples limit the precision of the modeling, and the resulting statistical uncertainty is incorporated into the shape uncertainty. A similar treatment is applied to the backgrounds estimated from control samples in data, based on the statistical uncertainties in the corresponding control samples. and expected events are shown in Table 3, which also includes the expectation for a selected parameter point for each type of signal. Figure 5 shows the E miss T distributions after the final selection. The observed distributions agree with the SM background predictions and no excess of events is observed.

The DM interpretation
The results are interpreted in the context of a simplified model of DM production. Figure 6 shows 95% confidence level (CL) expected and observed limits on the signal strength σ obs /σ th in the case of vector and axial-vector mediators and for two possible values of the quarkmediator coupling constant, g q = 0.25 or 1. Independent of the type of coupling, production of DM particles via an on-shell mediator (2m χ < M med ) can be excluded up to mediator masses of ≈400 GeV for g q = 1.0 and up to ≈300 GeV for g q = 0.25. Dark matter particle masses are probed up to 100-150 GeV for vector and up to 50-100 GeV for axial-vector couplings. For g q = 1.0, a small region of off-shell parameter space can also be excluded. In the case of g q = 0.25, sensitivity is limited to the on-shell region.
The simplified model allows a calculation of the DM relic abundance in the universe for each parameter point [72,73]. Parameter combinations consistent with measurements of the DM relic abundance in the universe are indicated in Fig. 6. For these parameter combinations, no BSM phenomena other than the simplified model are needed to account for the relic abundance in the universe. For other parameter values, additional phenomena, such as an extended dark sector, are necessary. for the final selection in the e + e − (left) and µ + µ − (right) channels. Expected signal distributions are shown for the simplified model of DM production with vector couplings, the EFT DM production benchmark, and unparticle model. The total uncertainty (stat. ⊕ sys.) in the overall background is shown as a hatched region. Overflow events are included in the rightmost bins. In the lower panels, the ratio between data and predicted background is shown.
The exclusion limits in the M med -m χ plane are translated into limits on the DM-nucleon scattering cross section using the prescription of Ref. [35]. The limits are set at 90% CL, assuming g q = 0.25. The resulting exclusion curves for both spin-independent (vector) and spindependent (axial-vector) cases are shown in Fig. 7, which compares them to the results from direct detection experiments. The comparison of collider and direct detection experiments highlights the complementarity of the two approaches. Especially in the case of lower DM masses and axial-vector couplings, a collider-based search can exclude parameter space not covered by direct detection experiments. In all cases, the DM-mediator coupling g χ is set to one. Figure 8 shows 95% CL expected limits on the cutoff scale Λ of the EFT benchmark model with DM pair coupling to gauge bosons. The limits are derived as a function of the DM particle mass. At low masses, cutoff scales up to ≈ 480 GeV can be excluded. With increasing DM particle mass, sensitivity decreases with Λ < 250 GeV excluded for m χ = 1.3 TeV. The 95% CL expected limits on the cutoff scale Λ and signal strength σ obs /σ th as a function of coupling c 1 and DM mass m χ are shown in Fig. 9. At c 1 ≈ 1, the interaction is dominated by the ZZχχvertex. With increasing c 1 , the γZχχ-vertex begins to contribute, yielding an improvement in the sensitivity.

Unparticle interpretation
In the unparticle scenario, 95% CL lower limits are set on the effective cutoff scale Λ U . A fixed coupling λ = 1 is assumed. The limits on Λ U are shown in Fig. 10 as a function of the scaling dimension d U . The result is compared with the limits obtained from previous CMS searches in the monojet [15] and mono-Z [17] channels, as well as with a reinterpretation of LEP searches [83]. Comparable sensitivity to the previous CMS mono-Z search is achieved owing to the increase in collision energy, which offsets the larger size of the previous dataset.  Figure 6: The 95% CL observed limits on the signal strength σ obs /σ theo in both vector (left) and axial-vector (right) mediator scenarios, for mediator-quark coupling constant values g q = 0.25 (upper) and 1 (lower). In all cases, the DM-mediator coupling g χ is set to one. The expected exclusion curves for unity signal strength are shown as a reference, with black dashed lines indicating the expected ±1 s.d. interval due to experimental uncertainties. The red dashed lines show the influence of theory-related signal normalization uncertainties on the observed limits, which are estimated to be 15%. The solid line labeled "Ω c × h 2 = 0.12" identifies the parameter region where no additional new physics beyond the simplified model is necessary to reproduce the observed DM relic abundance in the universe [1, 35, 72-74].

Model-independent limits
As an alternative to the interpretation of the results in specific models, a simple counting experiment is performed to obtain model-independent expected and observed 95% CL upper limits on the visible cross section σ BSM vis = σ A for BSM physics processes, where A is the acceptance and is the identification efficiency for a hypothetical signal. The limits as a function of E miss T thresholds are shown in Fig. 11. Table 4 shows the total SM background predictions for the numbers of events passing the selection requirements, for different E miss T thresholds, compared with the observed numbers of events. The 95% CL expected and observed upper limits for the contribution of events from BSM sources are also shown. Since the efficiency of reconstructing potential signal events depends on the characteristics of the signal, the model-independent limits are not corrected for the efficiency. For the models considered in this analysis, typical efficiencies are in the range 50-70% (simplified DM model), 60-70% (EFT DM model), and  Figure 7: Observed 90% CL limits on the DM-nucleon scattering cross sections in both spinindependent (left) and spin-dependent (right) cases, assuming a mediator-quark coupling constant g q = 0.25 and mediator-DM coupling constant g χ = 1. The line shading indicates the excluded region. Limits from the LUX [75], CDMSLite [76], PandaX-II [77], and CRESST-II [78] experiments are shown for the spin-independent case. Limits from the Super-Kamiokande [79], PICO-2L [80], PICO-60 [81], and IceCube [82] experiments are shown for the spin-dependent case.

Summary
A search for physics beyond the standard model has been performed in events with a Z boson and missing transverse momentum, using a data set corresponding to an integrated luminosity of 2.3 fb −1 of pp collisions at a center-of-mass energy of 13 TeV. The observed data are consistent with the expected standard model processes. The results are analyzed to obtain limits in three different scenarios of physics beyond the standard model. In a simplified model of DM production via a vector or axial vector mediator, 95% confidence level limits are obtained on the masses of the DM particles and the mediator. Limits on the DM-nucleon scattering cross section are set at 90% confidence level in spin-dependent and spin-independent coupling scenarios. In an effective field theory approach, limits are set on the DM coupling parameters to U(1) and SU(2) gauge fields and on the scale of new physics. For an unparticle model, 95% confidence level limits are obtained on the effective cutoff scale as a function of the scaling dimension. In addition, model-independent limits on the contribution to the visible Z + E miss T cross section from non-standard-model sources are presented as a function of the minimum requirement on E miss T . These results are the first in this signal topology to be interpreted in terms of a simplified model. Furthermore, the limits on unparticle production are the first of their kind to be presented at √ s = 13 TeV.

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.      [20] ATLAS Collaboration, "Search for dark matter produced in association with a hadronically decaying vector boson in pp collisions at √ s = 13 TeV with the ATLAS detector", Phys. Lett