Search for high-mass resonances in dilepton final states in proton-proton collisions at $\sqrt{s}=$ 13 TeV

A search is presented for new high-mass resonances decaying into electron or muon pairs. The search uses proton-proton collision data at a centre-of-mass energy of 13 TeV collected by the CMS experiment at the LHC in 2016, corresponding to an integrated luminosity of 36 fb$^{-1}$. Observations are in agreement with standard model expectations. Upper limits on the product of a new resonance production cross section and branching fraction to dileptons are calculated in a model-independent manner. This permits the interpretation of the limits in models predicting a narrow dielectron or dimuon resonance. A scan of different intrinsic width hypotheses is performed. Limits are set on the masses of various hypothetical particles. For the Z$'_\mathrm{SSM}$ (Z$'_{\psi}$) particle, which arises in the sequential standard model (superstring-inspired model), a lower mass limit of 4.50 (3.90) TeV is set at 95% confidence level. The lightest Kaluza-Klein graviton arising in the Randall-Sundrum model of extra dimensions, with coupling parameters $k/\overline{M}_\mathrm{Pl}$ of 0.01, 0.05, and 0.10, is excluded at 95% confidence level below 2.10, 3.65, and 4.25 TeV, respectively. In a simplified model of dark matter production via a vector or axial vector mediator, limits at 95% confidence level are obtained on the masses of the dark matter particle and its mediator.


Introduction
Neutral resonances decaying to lepton pairs occur in a variety of theoretical models that attempt to extend the standard model (SM) of particle physics. Grand unified theories (GUTs), including superstring and left-right-symmetric models (LR), achieve unification of the three forces at a high energy scale and predict the existence of new neutral gauge bosons [1,2]. These bosons might be light enough to be produced at current or future colliders. Theories that allow the gravitational force to propagate into extra spatial dimensions [3] could explain the large separation between the electroweak symmetry breaking energy scale and the gravitational energy scale. In such models, graviton excitations could be observed as spin-2 high-mass resonances. Moreover, high mass neutral resonances are predicted in models where dark matter (DM), whose existence is suggested by many astrophysical and cosmological observations [4], has a particle explanation. In these theories, interactions between the DM and SM particles could be mediated by high-mass, weakly coupled particles. Indirect evidence for DM could potentially be seen at the CERN LHC by searching for a high-mass DM mediator decaying to a dilepton final state.
Typical models predicting extra Z-like bosons extend the gauge group of the SM by additional U (1) gauge groups. The U (1) gauge groups, or a linear combination of them, can be broken near the TeV scale, giving rise to new massive gauge bosons denoted as Z . A generalized version of these models uses a continuously varying angle to describe the mixing of the U (1) generators. The cross section for charged lepton pair production via a Z vector boson can, in the narrow-width approximation (NWA), be expressed in terms of the quantity c u w u + c d w d [5,6]. The parameter c u (c d ) contains information about the model-dependent Z boson couplings to the up-type (down-type) quarks, while w u (w d ) depends on the up-type (down-type) quark parton distribution functions (PDFs). The parameterization of the linear mixing of the relevant U (1) generators produces a contour in the (c d , c u ) plane that represents each class of models. Commonly considered models are the generalized sequential model (GSM) [6], containing the Z SSM boson that has SM-like couplings to SM fermions [7]; GUT models based on the E 6 gauge group, containing the Z ψ boson [1,8]; and high-mass neutral bosons of the left (L)-right (R) symmetric extensions of the SM based on the SU(2) L ⊗ SU(2) R ⊗ U(1) B-L gauge group, where B-L refers to the difference between baryon and lepton numbers. Specific choices of the mixing angle produce different models such as those shown in Table 1, with their exact definition given in Ref. [6]. The (c d , c u ) plane parameterization provides a model-independent way to create a direct correspondence between the experimental bounds on Z production cross sections and the parameters of the Lagrangian. The translation of the experimental limits into the (c d , c u ) plane is studied both in the context of the NWA and by taking finite widths into account. The two procedures, NWA and finite widths, have been shown to give the same results [6]. A further study including the effects of interference [9] has demonstrated that the two procedures can still be used with an appropriate choice of the invariant mass window, within which the cross section is calculated.
Searches for high-mass Z gauge bosons have been performed by the CMS Collaboration at the LHC with proton-proton collision data collected at √ s = 7 TeV [10,11] and with data collected at 8 TeV [12,13]. More recently CMS performed a search using the combination of 2015 data collected at 13 TeV with data collected at 8 TeV [14]. Searches for high-mass Z gauge bosons have also been performed by the ATLAS Collaboration with data collected at 7 TeV [15,16], with data collected at 8 TeV [17], and with data collected at 13 TeV [18].
Kaluza-Klein graviton (G KK ) excitations arising in the Randall-Sundrum (RS) model of extra spatial dimensions [3,19] involve a finite five-dimensional bulk that is warped as a function We also consider a simplified model with a single DM particle that has sizeable interactions with the SM fermions through an additional spin-1 high-mass particle mediating the SM-DM interaction. In this model, the mediator is a vector or an axial-vector boson and is exchanged in the s channel [26,27]. There are five free parameters in this model: the DM mass m DM , the mediator mass m Med , the coupling g DM between the mediator and the DM particle, and the universal couplings g and g q between the mediator and the SM charged leptons and quarks, respectively. These five parameters define the production rate of the mediator, its DM and leptonic/hadronic decay rates, and the kinematic distributions of the signal events. We investigate two sets of benchmark coupling values that illustrate the complementary strengths of dijet and dilepton searches and the typical impact of searches for dilepton resonances in this model [26]: The results presented in this paper are obtained from an analysis of the data sample collected in 2016 at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 for the dielectron channel and 36.3 fb −1 for the dimuon channel. The invariant mass spectra of the observed dilepton final states are scrutinized for possible deviations from the SM background predictions. Background yields are estimated with simulated samples and normalized to their relative cross sections using the highest order calculations available. The sum of these backgrounds is normalized to the observed yield in the dilepton invariant mass region of 60-120 GeV. An exception to this is the background from multijet events, which is estimated from the data using control regions. Additionally, the tt background prediction from simulation is also crosschecked via eµ events in data, as discussed in Sect. 5.
Limits are set on the ratio of the cross section for dilepton production via a new boson to the cross section for dilepton production via the SM Z boson. This is done in order to remove the dependence on the CMS integrated luminosity measurement and to suppress the correlated uncertainties between the low-and high-mass regions. The computation of the observed limit and significance involves an arbitrary choice of the intrinsic width of the new high-mass resonance. The choice of the intrinsic width can potentially affect the statistical interpretation of the result, therefore we provide limits by scanning different width hypotheses, as discussed in Sect. 6. The analysis is designed to minimize the effect of the specific model assumptions on the results, allowing the results to be interpreted in the framework of any high-mass resonance decaying into lepton pairs with a width and pseudorapidity distribution similar to the reference model used. The couplings of the Z model will not only impact the width of the peak region of the high-mass resonance but also the tails due to PDFs and interference effects between the SM electroweak bosons and the new resonance in DY Z/γ * → e + e − /µ + µ − processes. In the NWA chosen for this analysis, the contributions to the signal cross section from PDFs and interference off-shell effects, which are highly model dependent, are removed. Therefore, in this work we consider only the resonant peak, taking into account the effects of the intrinsic width of the high-mass resonance on the experimentally observed mass peak region. This paper is structured as follows. Section 2 contains a short description of the CMS detector. In Section 3, we discuss the data sample and the Monte Carlo (MC) event generators used for the signal and background simulation. The lepton reconstruction and event selection relevant for the Z boson search are discussed in Section 4. The background estimation methods are described in Section 5. In Section 6 the statistical method used to extract the results and the statistical treatment of the systematic uncertainties are explained. Results are summarized in Section 7.

The CMS detector
The central feature of the CMS detector is a superconducting solenoid providing an axial magnetic field of 3.8 T and enclosing an inner tracker, an electromagnetic calorimeter (ECAL), and a hadron calorimeter (HCAL). The inner tracker is composed of a silicon pixel detector and a silicon strip tracker, and measures charged particle trajectories in the pseudorapidity range |η| < 2.5. The ECAL and HCAL, each composed of a barrel and two endcap sections, extend over the range |η| < 3.0. The finely segmented ECAL consists of nearly 76 000 lead tungstate crystals, while the HCAL is constructed from alternating layers of brass and scintillator. Forward hadron calorimeters encompass 3.0 < |η| < 5.0. The muon detection system covers |η| < 2.4 with up to four layers of gas-ionization detectors installed outside the solenoid and sandwiched between the layers of the steel flux-return yoke. Additional detectors and upgrades of electronics were installed before the beginning of the 13 TeV data collection period in 2015, yielding improved reconstruction performance for muons relative to the 8 TeV data collection period in 2012. The efficiency to reconstruct and select muons that result from Z boson decays and pass specific identification selection criteria, has increased by approximately 2% between Run 1 and Run 2 as a result of these upgrades [28]. 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. [29].
The CMS experiment has a two-level trigger system. The level-1 (L1) trigger [30], composed of custom hardware processors, selects events of interest using information from the calorimeters and muon detectors and reduces the readout rate from the 40 MHz bunch crossing frequency to a maximum of 100 kHz. The software based high-level trigger (HLT) [30] uses the full event information, including that from the inner tracker, to reduce the event rate to around the 1 kHz that is retained for further processing.

Simulated data samples
The dominant background in this search is the DY process. The simulated DY background is generated with POWHEG v2 [31][32][33][34][35][36] from next-to-leading order (NLO) matrix elements using the NNPDF3.0 [37] PDF set, and with PYTHIA 8.205 [38] for parton showering and hadronization. For all simulated SM samples, the default tune for PYTHIA, CUETP8M1 [39], is used. The DY cross section at NLO is corrected to next-to-next-to-leading order (NNLO) in perturbative quantum chromodynamics (QCD) by using a dilepton invariant mass dependent K-factor according to the predictions of the FEWZ 3.1.b2 program [40]. In addition, these predictions incorporate missing EW corrections at NLO. For the FEWZ calculations, the LUXqed plus -PDF4LHC15 nnlo 100 [41] PDF set is used in which the QCD PDFs based on the PDF4LHC [42] set are combined with the photon PDFs to account for pure quantum electrodynamics effects. Another nonresonant background arises from a γγ initial state via t and u channel processes. The photon-induced (PI) process produces two leptons in the final state [43,44]. This contribution is included in the K-factor that corrects the DY NLO cross section.
The tt, tW and WW backgrounds are simulated using POWHEG v2, with parton showering and hadronization described by PYTHIA 8.205. The NNPDF3.0 PDF set is used for all these samples. The tt cross section is calculated at NNLO with TOP++ [45] assuming a top quark mass of 172.5 GeV. The inclusive diboson processes WZ, and ZZ are simulated at leading order (LO) using the PYTHIA 8.205 program along with the NNPDF3.0 PDFs. The production of DY τ + τ − and W+jets is simulated at LO with the MADGRAPH5 aMC@NLO version 2.2.2 [46] program. The PDFs are evaluated using the LHAPDF library [47][48][49].
We use a sample of events in which a Z ψ boson is generated with a mass of 3000 GeV, and RS samples with the graviton generated at different mass values from 250 to 4000 GeV. In all samples the high-mass resonances decay to electron and muon pairs. Both signal samples are generated using the PYTHIA 8.205 program with the NNPDF3.0 PDFs. The Z ψ sample is used to create simulated peaks in the dilepton mass plots of Fig. 1. However, we use DY samples to model the Z at high masses, since the dilepton behaviour in this region is identical in the two cases.
The presence of additional pp interactions in the same or adjacent bunch crossing observed in data (pileup) is incorporated in simulated events by including overlapping pp interactions. MC samples are corrected to reproduce the pileup distribution as measured in data (pileup reweighting), with an average number of pileup interactions per proton bunch crossing of approximately 22 for the 2016 data sample. The detector response is simulated using the GEANT4 [50] package.

Lepton reconstruction and event selection
The electron and muon reconstruction algorithms and event selection criteria used in this highmass dilepton search are mostly unchanged from the previous analysis [14]. However, the muon selection criteria were modified in both the online and offline selection in order to increase the efficiency in the high mass region, above 1 TeV.
Energy deposits in the ECAL are combined into clusters under the assumption that each local maximum represents a single particle. Any clusters consistent with originating from a single particle that may have undergone bremsstrahlung emission are grouped together. If a track from the nominal interaction point is geometrically associated with a cluster, this track together with the cluster form an electron candidate. The angular information of the electron candidate is taken from the track. The energy of the electron uses only the ECAL deposits and is not combined with the track momentum. Electron candidates are required to have a transverse momentum p T > 35 GeV and satisfy |η C | < 1.44 (ECAL barrel region) or 1.57 < |η C | < 2.50 (ECAL endcap region), where η C is the pseudorapidity of the cluster of ECAL deposits comprising the electron with respect to the nominal centre of the CMS detector. The transition region 1.44 < |η C | < 1.57 is excluded as it leads to lower-quality reconstructed clusters, owing mainly to services and cables exiting between the barrel and endcap calorimeters.
The electron candidates are also required to pass a set of dedicated high energy electron selection criteria [51]. This selection requires that the lateral spread of deposits in the ECAL be consistent with that of a single electron, that the track be matched to the ECAL deposits and be consistent with a particle originating from the nominal interaction point, and that the associated energy in the HCAL around the electron direction be less than 5% of the reconstructed energy of the electron, once noise and pileup are taken into account. The selection also requires that the electron be isolated in a cone of radius ∆R = √ (∆η) 2 + (∆φ) 2 = 0.3 in both the calorimeter and tracker [14]. Only well-measured tracks that are consistent with originating from the same vertex as the electron are included in the isolation sum.
The efficiency of the trigger to select events with two electrons passing the analysis selection requirements is 98.5%, when the barrel (endcap) electrons satisfy p T > 36 (38) GeV. This primary trigger is monitored by a suite of higher-threshold triggers which have progressively fewer selection requirements, culminating in a trigger that simply requires 800 GeV of p T in the ECAL. These triggers are included to minimize the chance that unexpected reconstruction problems cause a lower than expected efficiency in the primary trigger.
In the selection of dielectron pairs, at least one of the two electrons must be in the ECAL barrel region in order to reduce the background from multijet events. This also allows the endcapendcap events to be used as a control sample for the QCD background estimate. Dielectron pairs are not required to be oppositely charged as this leads to a significant efficiency loss at high invariant mass [10]. As the electron energy is solely obtained from the calorimeter, an incorrectly measured charge does not impact the measured mass. If there are multiple possible dielectron pairs, the pair with the two highest p T electrons is selected.
The efficiency to trigger, reconstruct, and select an electron pair with invariant mass equal to 1 TeV within the detector acceptance is 69 (65)% for barrel-barrel (barrel-endcap) events. The trigger efficiency is measured in data and the total efficiency is estimated using simulated DY events and validated using data measurements at the Z boson mass peak. Using Z bosons, it is possible to probe the efficiencies up to electron p T of 500 GeV within a few percent precision. The uncertainty in the efficiency is ±3 (±5)% in the barrel (endcap) for electrons coming from a dielectron pair with a mass m > 120 GeV. The simulated efficiency reproduces the energy evolution of the observed efficiency in the measurable region from 40 to 500 GeV.
A candidate muon pair at the L1 trigger is required to have at least one muon reconstructed with segments in the muon detectors and transverse momentum p T above 22 GeV. These muons are then required to have p T > 50 GeV and |η| < 2.4 at the HLT. A prescaled HLT path is used to extract the observed yield in the invariant mass region of the Z boson peak (60 < m < 120 GeV) in order to construct the normalization factor to that region. This prescaled trigger has a p T threshold of 27 GeV with the same L1 requirements as the main trigger.
The trigger efficiency for dimuon events, where both muons have p T > 53 GeV, is parameterized using simulated DY events. It is measured to be around 99.5% if both muons are in the barrel (|η| < 1.2), and 99.0% for events with one muon in the endcap (|η| > 1.2). These efficiencies are validated as a function of muon p T and η, using data events that pass the full offline selection criteria. This is done using muons present in high mass dilepton or high-p T Z boson events (free from background contributions) and other data sets (selected by electron, missing transverse momentum, or jet requirements). The measurements are found to be in agreement for muons with p T up to 1.5 TeV and an uncertainty of ±0.3 (±0.7)% for barrel-barrel (barrelendcap and endcap-endcap) dimuon events over the full mass range is assigned.
High-p T (>200 GeV) muon offline reconstruction uses dedicated algorithms [52] to take into account the effects of radiative processes of high-energy muon interactions with the detector material. Muon candidates are required to have p T > 53 GeV, to be within the region of |η| < 2.4, and to pass dedicated high-momentum identification selection criteria [13]. Muon candidates are also required to pass isolation requirements. The summed p T of tracks within a cone of radius ∆R = 0.3 around the candidate direction should be less than 10% of the p T of the candidate, excluding from summation the muon candidate under consideration.
Oppositely charged muon candidates passing the selection are combined to form dimuon candidates. A χ 2 fit is performed to the common vertex between the two muons to ensure that they originate from the same vertex. This fit is required to have reduced χ 2 < 20. The angle between the directions of the two muon candidates is required to be less than π − 0.02 in order to suppress cosmic ray backgrounds. If there are multiple possible dimuon pairs, the pair with the two highest p T muons is selected.
Standalone muon reconstruction and identification efficiencies as a function of the muon momentum, both in data and simulation, were probed using high quality isolated inner tracks extrapolated to the muon system. In the barrel region MC efficiencies are validated up to 1.2 TeV, while in the endcaps we observed a small discrepancy with respect to data at the order of 10% for muon momentum of 3 TeV. The efficiency ratios between data and MC obtained as a function of the single muon momentum are then propagated as a function of the dimuon mass. It leads to −1.5 (−6.5)% one-sided uncertainties on a 4 TeV mass for barrel-barrel (barrel-endcap and endcap-endcap) events. This represents the dominant uncertainty in the signal acceptance times efficiency.
The efficiency to trigger, reconstruct, and select a muon pair with invariant mass equal to 1 TeV within the detector acceptance is 92.7 +0.3 −0.5 % (92.5 +0.7 −2.7 %) for barrel-barrel (barrel-endcap and endcap-endcap) events. The efficiency for each of the two muons is correlated as is the case for electron pairs. The experimental dilepton mass resolution is determined from simulation as a function of the generated dilepton mass. The simulated mass resolution measured in simulations is smeared in order to be comparable to Z boson events selected in data. This smearing is performed for electrons in the entire p T range while for the muon channel this smearing is performed as a function of the leading muon p T up to 800 (450) GeV in the barrel-barrel (barrel-endcap and endcap-endcap) events. The experimental mass resolution is 1.0 (1.5)% for barrel-barrel (barrel-endcap) electron pairs with a mass of 1 TeV. No uncertainty is assigned to the electron pair mass resolution as we use the resolution measured in the low-mass region, which also represents the most pessimistic value for higher masses. The resolution for muon pairs with a mass of 1 TeV is 3.0 ± 0.5 (4.0 ± 0.6)% for barrel-barrel (barrel-endcap and endcap-endcap) events. We assign a 15% systematic uncertainty in these values, based on measurements with different functions to parameterize the resolution for muon pairs.
The response of the detector to leptons might depend on the increasing dilepton invariant mass. For electrons this could reflect a nonlinear response of the readout electronics. There is no evidence for such effects in the current data. The energy scale uncertainty at high p T (400 GeV) is validated at the 2 (1)% level for electrons in the barrel (endcaps). Since the calorimeter is expected to behave linearly with energy we therefore assume this is true for a 2 TeV mass. As the muon p T increases, it becomes increasingly sensitive to the detector alignment. New methods have been developed for the 2016 data to determine a potential bias from this source. The usage of muon alignment position uncertainties has been added in the HLT and offline muon reconstruction algorithms in order to compensate for misalignment. The curvature distributions of positive and negative muons (q/p T ) in data are compared to those obtained in simulation for different η and φ ranges. The mass scale is within 1 (3)% for barrel-barrel (barrel-endcap and endcap-endcap) dimuon events for the entire mass range up to 2.3 TeV, covering the region where dimuon events are observed.

Backgrounds
The dominant SM background arises from the DY process. In addition to the DY process, e + e − and µ + µ − pairs can be produced in the PI process γγ → + − [43] from photons radiated by the incoming protons. The PI process was studied through investigations of proton PDFs that include photon contributions [43,44]. Although the relative contribution of PI processes increases with dilepton mass, the effect on the statistical analysis of the data was found to be negligible. Uncertainties in the PDFs, in the contributions from PI processes, and in the NNLO corrections to the cross sections dominate the systematic uncertainty in the amount of estimated background. This in turn dictates the uncertainty in the background shape. The uncertainty due to the PDFs is assessed using the PDF4LHC prescription [42] and is found to vary from 1.4 to 20% as the dilepton mass increases from 0.2 to 6 TeV.
Other sources of background are real leptons from the top quark-antiquark (tt), single top quark (tW), diboson (WW, WZ, and ZZ), and DY τ + τ − processes. The contribution of these sources is reduced at high invariant dilepton mass. These backgrounds are estimated using simulated events. For tt, tW, diboson, and DY τ + τ − production, the yield of eµ final states produced should be approximately equal to the sum of ee and µµ final states. Therefore the simulation predictions in these channels are compared to data in the eµ final state. The predictions from data and simulation are in agreement within uncertainties and no further action is taken.
Multijet-enriched data control samples are used to evaluate the contribution of jets misidentified as electrons, as described in Ref. [13]. Similarly, multijet data control samples are used to evaluate the probability for jets and nonisolated muons to be identified as isolated muons. The contribution of misidentified jets to the total background is 1-3%; therefore even with large uncertainties up to 50%, it has a negligible effect on the statistical analysis of the data.
The contribution of cosmic ray background events is negligible in this analysis due to the event selection, and is <0.2% for muons with p T > 300 GeV.
The observed invariant mass spectra of the dielectron and dimuon events are presented in Fig. 1. The highest observed dilepton mass is 2.6 TeV and appears in the dielectron final state as a barrel-endcap event. The highest observed dimuon mass is 2.3 TeV and appears as a barrelendcap event. The structure observed just above the Z boson peak in the dimuon channel is due to the high threshold (p T > 53 GeV) applied to the transverse momentum of the muon candidates. The uncertainty bands in the ratio plots represent the systematic uncertainty in the background yields arising from the dilepton mass scale, trigger efficiency, acceptance times efficiency, Z boson normalization (1% for dielectron and 5% for dimuon channel), PDFs, non-DY background cross section determination (7%), and lepton misidentification rates, summed in quadrature. The selection criteria for the dielectron channel have a small sensitivity to pileup, so the contribution from pileup reweighting is included in the uncertainty band. The dimuon selection criteria are by construction insensitive to pileup.
The corresponding cumulative distributions are shown in Fig. 2. The SM expected yields in various mass bins are compared to the observed yields in Tables 2 and 3. Observations agree with expectations in the entire mass region for the dielectron events. A deficit of dimuon events is observed in the high-mass region compared to the expectations from SM processes (as shown on the right plot of Fig. 2). This deficit appears in the barrel (|η| < 1.2). In the barrel-barrel category we observe two dimuon candidates in the region m µµ > 1600 GeV, where we expect ten dimuon candidates from MC simulations. This leads to a local significance of the discrepancy equal to 2.9 standard deviations (s.d.). This significance is reduced to 1.8 s.d. when considering the entire pseudorapidity range and is considered to be compatible with a statistical fluctuation. No experimental sources for the deficit were identified. Table 2: The number of dielectron events in various invariant mass ranges. The total background is the sum of the events for the SM processes listed. The yields from simulation are normalized relative to the expected cross sections, and overall the simulation is normalized to the observed yield using the number of events in the mass window 60-120 GeV. Uncertainties include both statistical and systematic components, summed in quadrature.

Statistical analysis, results, and interpretation
The mass distributions are scrutinized for possible deviations from the SM background predictions. No significant deviations are observed. The limits are expressed as a function of R σ , which is the ratio of the cross section for dilepton production via a Z boson to the measured  cross section for dilepton production via the Z boson in the mass window 60-120 GeV: Expressing the limits as a ratio, reduces the dependency on the theoretical prediction of the Z boson cross section as well as the correlated experimental uncertainties.
For the electron and muon channel combination, the branching fractions of these two channels are assumed to be the same. The signal cross section corresponds to that obtained in the narrow width approximation; specifically, off-shell contributions from PDFs and interference effects are not included. The limits are set using a Bayesian method with an unbinned extended likelihood function [13] using the framework developed for statistically combining Higgs boson searches [53], which is based on the ROOSTATS package [54]. The signal probability density function (pdf) used is a convolution of a Breit-Wigner (BW) function and a Gaussian function with exponential tails to either side (Cruijff [55]). The BW function models the intrinsic width of the particle, while the Cruijff function models the detector response. A Crystal Ball (CB) function [56] better describes the signal shape at higher dielectron masses compared to the Cruijff function. Therefore for m ee > 2300 GeV, a CB function is used. The background pdf for the dielectron channel is provided by the following formula: while for the dimuon channel the following functional form is used: The different background pdfs for dielectrons and dimuons reflect the different background composition in the two channels. For each final state, the parameters of the background pdf are obtained by fitting the total background distribution produced using SM MC generators and the background arising from misidentified jets deduced from the data. The fits to the background distribution are set for masses above 120 GeV.
For the signal cross section we use a positive uniform prior. The systematic uncertainties in the dilepton mass originating from efficiencies, resolution and scale, as discussed in Sect. 4, are treated as nuisance parameters and assigned log normal priors. The relative mass scale of the different channels is the only uncertainty with a noticeable impact. The limits are calculated in a mass window of ±6 times the signal width, with this window being symmetrically enlarged until there is a minimum of 100 data events in it. This procedure sets the level of the statistical uncertainty in the local background amplitude; the level is chosen to dominate the expected systematic uncertainties in the background shape at high mass. The total background uncertainty ranges from approximately 3% at 200 GeV to 12% at 5 TeV for both electrons and muons. At low mass it is driven by normalization uncertainties; while at high mass, PDFs and higher-order corrections are the dominant sources of uncertainty.
The expected and observed limits for a resonance width equal to 0.6% of the resonance mass are shown in Fig. 3 for the dielectron channel, dimuon channels, and their combination. Table 4 presents the observed and expected 95% confidence level (CL) lower limits on the masses of spin-1 Z SSM and Z ψ bosons. Results for widths equal to 0.6, 3, 5 and 10% of the resonance mass are shown in Fig. 4 for the dielectron channel, dimuon channel, and their combination. For masses below 2 TeV the expected limits become less stringent with increasing resonance widths. At high masses, however, the experimental mass resolution dominates and the limits do not exhibit any dependence on the assumed resonance width. Compared to the default width of 0.6%, an increased width results in less stringent observed and expected limits, and a smoother variation of the observed limit as a function of mass.
The R σ curves are shown on the plots to obtain mass limits for Z signal models. The curves are constructed by dividing the LO cross section of a given model, calculated using the PYTHIA 8.2 program with the NNPDF2.3 [57] PDFs, by the NNLO Z boson cross section of 1928 ± 73 pb obtained with FEWZ 3.1 [40]. As the limits presented here are set on the on-shell cross section and the PYTHIA event generator includes off-shell effects, the cross section is calculated in a mass window of ±5% √ s centred on the resonance mass, following the prescription of Ref. [9]. The validity of this procedure for the Z SSM and Z ψ bosons was explicitly checked in Ref. [9] and is found to be accurate to approximately 5-7%. To account for NNLO QCD effects, the LO cross sections are multiplied by a mass independent K-factor. The value of the K-factor is estimated at a dilepton mass of 4.5 TeV and found to be consistent with unity. Applying a mass dependent K-factor, the Z ψ resonance mass limit differs by only 50 GeV, justifying the use of the simpler mass independent K-factor.    The upper limits at 95% CL on the product of production cross section and branching fraction for a spin-1 resonance with a width equal to 0.6% of the resonance mass, relative to the product of production cross section and branching fraction of a Z boson, for the dielectron channel (left), dimuon channel (right), and their combination (lower). The shaded bands correspond to the 68 and 95% quantiles for the expected limits. Theoretical predictions for the spin-1 Z SSM and Z ψ resonances are shown for comparison. the range 0.5-0.6 for the results shown here) and hence a single value of the cross section is represented by a straight line in the (c d , c u ) plane. In the log-log plot shown in Fig. 5, the thin lines labelled with a mass value correspond to the cross section limit at that mass. The closed contours representing the GSM, LR, and E 6 model classes are composed of thick line segments which correspond to ranges of the particular mixing angle for each considered model. A brief description of the models is given in Section 1 with further information provided in Table 1 and the exact definition of the models discussed in Ref. [6]. The mass limit on the relevant Z boson in any model, where c d and c u have been determined, can be read off this plot. For completeness we quantify a possible presence of an excess of events over what is expected for the background by computing the p-value. The p-value for different signal width hypotheses is shown for both the separate and combined channels in Fig. 6. The largest excess in the combined result is observed around M = 1300 GeV having a local significance of around 2.5 s.d. for a spin-1 resonance with widths 0.6 to 5.0%. This corresponds to a global significance of −0.92 s.d. after taking into consideration the look elsewhere effect [58] in the mass range 200 to 5500 GeV. The global significance is expressed as the corresponding number of standard deviations using the one-sided Gaussian tail convention. The methodology of the p-value computation is described in Ref. [59]. Obs. 95% CL limit  Obs. 95% CL limit   The upper limits at 95% CL on the product of production cross section and branching fraction for a spin-1 resonance, for widths equal to 0.6, 3, 5, and 10% of the resonance mass, relative to the product of production cross section and branching fraction for a Z boson, for the dielectron channel (left), dimuon channel (right), and their combination (lower). Theoretical predictions for the spin-1 Z SSM and Z ψ resonances are also shown. and 1.42 GeV corresponding to coupling parameters k/M Pl of 0.01, 0.05, and 0.10, are shown in Fig. 7 for the dielectron channel, dimuon channel, and their combination. Table 5 presents the values of the observed and expected 95% CL lower limits of the aforementioned models. The signal production cross sections, calculated using the PYTHIA 8.2 program with the NNPDF2.3 PDFs at LO, are multiplied by a K-factor of 1.6 to account for NLO effects [60]. The PI contribution to the production cross sections is small enough to be ignored. The results are also interpreted in the context of a simplified model with a DM particle that has sizeable interactions with SM fermions through an additional spin-1 high-mass particle mediating the SM-DM interaction. In the simplified model under consideration [26], only one DM particle exists, which is assumed to be a Dirac fermion. Limits are presented in Fig. 8 for two

GSM
Mixing angle: (-π/2,-π/4) (-π/4,0) (0,π/4) (π/4,π/2) cases with different sets of benchmark coupling values. The first case corresponds to a vector mediator with small couplings to leptons while the second one corresponds to an axial-vector mediator with equal couplings to quarks and leptons. The cross sections for lepton production are calculated at NLO in QCD, using the DMSIMP implementation [27] of the simplified model in MADGRAPH5 aMC@NLO version 2.5.2 [46]. Assuming the optimistic axial-vector coupling scenario and m DM > m Med /2, the signal cross section for the production of an electron or muon pair within the analysis acceptance ranges between approximately 100 pb at low values of the mediator mass (around 200 GeV), and 0.1 fb for higher values (around 4 TeV). The partial and total mediator decay widths, calculated at LO in QCD, are included via the MADWIDTH package [61].
While the DM particle is not probed directly, its mass indirectly modulates the sensitivity of the dilepton search. For low values of the DM particle mass m DM < m Med /2, the mediator boson will dominantly decay into DM particles, thus reducing the branching fraction to leptons, and making the mediator harder to probe in this search. At high values of the DM particle mass m DM > m Med /2, the mediator cannot decay to the DM particles and the leptonic branching fraction becomes sizeable. In the vector mediator model, the relatively small leptonic couplings mostly limit the sensitivity of this analysis to the regime of m DM > m Med /2. This regime is especially interesting to probe since it is almost inaccessible to typical searches based on missing transverse momentum [64]. In the axial-vector mediator model, the leptonic couplings of the mediator are sizeable and an exclusion is also possible for m DM < m Med /2. In the former case, the limit on the mediator mass reaches up to 1.8 TeV, depending on the mass of the DM particle; while in the latter case the limit reaches the 3.0-4.0 TeV, depending again on the mass of the DM particle. In the vector mediator model, the observed exclusion reaches up to values of the mediator mass equal to approximately 1.8 (0.6) TeV above (below) the diagonal, m DM = m Med /2. In the region where the mediator mass is equal to approximately 1.3 TeV, an upward fluctuation in the data (Fig. 3) results in a small above-diagonal region that is not excluded. Assuming that there is no new physics other than the mediator and DM particle, the relic density of DM in the universe, Ωh 2 , can be calculated. Regions of parameter space that reproduce the observed value Ωh 2 ≈ 0.12 [4], are indicated in Fig. 8. Numerical values are obtained from Ref. [26], where they were calculated using MADDM in version 2.06 [62,63]. In the hatched area, the DM would be overabundant in the universe. However, it is unlikely that new physics is fully described by the simplified model considered here, and the relic density constraint is not a stringent constraint, as additional new phenomena may modify its calculated value.    Figure 8: Limits at 95% confidence level for the masses of the DM particle, which is assumed to be Dirac fermion, and its associated mediator, in a simplified model of DM production via a vector (left) or axial vector (right) mediator. The parameter exclusion is obtained by comparing the limits on product of the production cross section and the branching fraction for decay to a Z boson, with the values obtained from calculations in the simplified model. For each combination of the DM particle and mediator mass values, the width of the mediator is taken into account in the limit calculation. The lines with the hatching represents the excluded regions. The solid grey lines, marked as "Ωh 2 ≥ 0.12", correspond to parameter regions that reproduce the observed DM relic density in the universe [4,26,62,63], with the hatched area indicating the region where the DM relic abundance exceeds the observed value.

Summary
A search for narrow resonances in dielectron and dimuon invariant mass spectra has been performed using data recorded in 2016 from proton-proton collisions at √ s = 13 TeV. The integrated luminosity for the dielectron sample is 35.9 fb −1 and for the dimuon sample is 36.3 fb −1 .
Observations are in agreement with standard model expectations. Upper limits at 95% confidence level on the product of a narrow-resonance production cross section and branching fraction to dileptons have been calculated in a model-independent manner to enable interpretation in the framework of models predicting a narrow dielectron or dimuon resonance. A scan of different intrinsic width hypotheses is performed.
Limits are set on the masses of various hypothetical particles. For the Z SSM particle, which arises in the sequential standard model, and for the superstring-inspired Z ψ particle, 95% confidence level lower mass limits for the combined channels are found to be 4. Finally, limits at 95% confidence level are obtained for the masses of the dark matter particle and its associated mediator, in a simplified model of dark matter production via a vector or axial vector mediator.

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 [19] L. Randall