Search for resonant and nonresonant new phenomena in high-mass dilepton final states at $\sqrt{s} = $ 13 TeV

A search is presented for physics beyond the standard model (SM) using electron or muon pairs with high invariant mass. A data set of proton-proton collisions collected by the CMS experiment at the LHC at $\sqrt{s} =$ 13 TeV from 2016 to 2018 corresponding to a total integrated luminosity of up to 140 fb$^{-1}$ is analyzed. No significant deviation is observed with respect to the SM background expectations. Upper limits are presented on the ratio of the product of the production cross section and the branching fraction to dileptons of a new narrow resonance to that of the Z boson. These provide the most stringent lower limits to date on the masses for various spin-1 particles, spin-2 gravitons in the Randall--Sundrum model, as well as spin-1 mediators between the SM and dark matter particles. Lower limits on the ultraviolet cutoff parameter are set both for four-fermion contact interactions and for the Arkani-Hamed, Dimopoulos, and Dvali model with large extra dimensions. Lepton flavor universality is tested at the TeV scale for the first time by comparing the dimuon and dielectron mass spectra. No significant deviation from the SM expectation of unity is observed.


Introduction
The presence of new phenomena, both resonant and nonresonant, in the high-mass dilepton final state is predicted by various theoretical models aiming to extend the standard model (SM) of particle physics. In this paper, the production of new spin-1 or spin-2 resonances, as well as the nonresonant production of high-mass lepton pairs, is considered.
The existence of new neutral gauge bosons is a possible signature of grand unified theories, such as superstring and left-right-symmetric (LRS) models, that include a unification of the three forces at a high energy scale [1,2].
One persistent puzzle in modern particle physics is the large difference in the energy scale of electroweak symmetry breaking and the energy scale of gravitation. This could be explained in theories including spatial extra dimensions, where the gravitational force can propagate into additional dimensions. In models by Arkani-Hamed, Dimopoulos, and Dvali (ADD) [3,4], the SM particles are confined to the traditional four dimensions of space and time, while in models proposed by Randall and Sundrum (RS) [5] the SM particles could also propagate into the additional dimensions. Possible signatures of these theories at LHC energies are graviton excitations, either as distinct spin-2 high-mass resonances in the RS model, or as a series of nearly mass-degenerate excitations that result in an overall nonresonant excess of events at high mass in the ADD model.
Based on many astrophysical and cosmological observations, it is assumed that dark matter (DM) accounts for the majority of matter in the universe [6]. Models have been proposed in which the DM consists of particles that can interact with those of the SM via high-mass, weakly coupled mediator particles [7]. These mediators could then be observed via their decay into SM particles, including the dilepton final state.
It has long been speculated that the presence of three generations of quarks and leptons is a sign that these particles are not fundamental but rather composed of constituent particles commonly called "preons" [8]. At energies observable at collider experiments, the preons would be confined into bound states by a new interaction analogous to quantum chromodynamics (QCD). This new interaction is characterized by an energy scale Λ, at which its effects would be directly observable. At center-of-mass energies far below Λ, the presence of the preon bound states would manifest itself as a flavor-diagonal "contact interaction" (CI) [9], resulting in a nonresonant excess of events at high mass.
Hints for lepton flavor universality violation in several measurements recently reported by the LHCb Collaboration [10][11][12], together with other flavor anomalies in B-meson decays summarized in Ref. [13], have sparked interest in models for physics beyond the SM that could explain these effects. These include models with heavy neutral gauge bosons [2] or leptoquarks [14]. Some of these models would result in a significant deviation from unity of the ratio of the dimuon to dielectron differential cross section as a function of dilepton mass m , been performed by the ATLAS Collaboration with data collected at 7  [34][35][36]. The additional data in the dimuon channel are recorded with the detector in conditions unsuitable for electron reconstruction. Deviations from the SM predictions are searched for in the dilepton invariant mass spectrum for both dielectron and dimuon events. The contribution from dileptonic decays of tau lepton pairs to the signal is small and neglected. The shapes of contributions from SM processes are estimated from simulation, except for backgrounds containing leptons produced inside jets or jets misidentified as leptons, which are estimated from control regions in data. The simulated events describing the dominant Drell-Yan (DY) background are corrected to the highest order calculations available. The resulting background shape is normalized to the observed data yields in a mass window of 60-120 GeV around the Z boson peak, separately for the dielectron and dimuon channels. Two analysis strategies are followed, targeting resonant and nonresonant signatures. For resonant signatures, the search is performed in a mass window around the assumed resonance mass, whose size depends on the assumed intrinsic decay width of the resonance and the massdependent detector resolution. A range of masses and widths is scanned to provide results covering a wide selection of signal models. Unbinned maximum likelihood fits are performed inside the mass windows, allowing the background normalization to be determined from data. The background parametrizations are obtained from fits to the m distribution of the background estimates.
Upper limits are set on the ratio of the product of the production cross section and the branching fraction of a new narrow dilepton resonance to that of the SM Z boson. Thus, many experimental and theoretical uncertainties common to both measurements cancel out or are reduced, leaving only uncertainties in the ratio that vary with the dilepton mass to be considered.
In the case of nonresonant signatures, the event sample is divided into several bins in invariant mass and, improving upon previous CMS analyses, the scattering angle cos θ * in the Collins-Soper frame [37]. The result is then obtained from a combination of the individual counting experiments within these bins. While most of the sensitivity to new physics stems from the highest mass bins where the signal-to-background ratio is most favorable, the bins at lower masses still contain valuable information, for example on possible interference between the signal and SM backgrounds. In this approach, the signal is not normalized to the Z boson cross section, and the background estimate cannot be obtained from the data. Therefore, this part of the analysis is more affected by uncertainties in the signal and background modeling.
In addition to these two search strategies, lepton flavor universality is tested for the first time in these final states by measuring the dimuon to dielectron ratio R µ + µ − /e + e − and comparing it to the SM expectation of unity, including corrections for detector effects, lepton acceptances, and lepton efficiencies. This paper is structured as follows. The CMS detector is briefly described in Section 2. A description of the signal models is given in Section 3, followed by a description of simulated event samples used in the analysis in Section 4. The event reconstruction and selection is given in Section 5 and the estimation of SM backgrounds is described in Section 6. Systematic uncertainties are described in Section 7, followed by the results and their statistical interpretation in Sections 8 and 9, respectively. The paper is summarized in Section 10. Tabulated results are provided in HEPDATA [38].

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. Before the data taking period in 2017, an upgraded pixel detector was installed, adding an additional barrel layer closer to the interaction point and additional disks in the two forward parts of the detector [39]. 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. 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].
Events of interest are selected using a two-tiered trigger system [41]. 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 fixed 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.

Signal models
A large number of models predict new phenomena resulting in high-mass dilepton signatures. The available models may be classified depending on whether they predict resonant or nonresonant production.

Models with resonant signatures
We consider three types of models resulting in resonances at high dilepton masses: Z particles in a class of U (1) gauge group models, spin-2 gravitons in the RS model of extra dimensions, and spin-1 dark matter mediators. The most prevalent are the theories that introduce new U (1) gauge groups. The symmetry properties of these groups, or a combination of them, are broken at some energy scale so that the SM group structure emerges as the low-energy limit of the new theory. If this breaking happens at the TeV scale, new gauge bosons Z with masses accessible at the LHC could exist. The properties of the new Z bosons depend on the mixing of the U (1) generators, which can be described by a continuously varying angle within a class of models. Commonly considered classes are the generalized sequential model (GSM) [42], containing the sequential SM boson Z SSM that has SM-like couplings to SM fermions [43]; LRS 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 the baryon and lepton numbers [42], which predict high-mass neutral bosons; and grand unified theories based on the E 6 gauge group, containing the Z ψ boson [1,44].
In the narrow width approximation [45], the Z production cross section can be expressed as c u w u + c d w d , where c u (c d ) are its couplings to up-type (down-type) quarks and w u (w d ) are the parton distribution functions (PDFs) for these quarks [42,46]. As the values of the modeldependent couplings depend on the aforementioned mixing of the U (1) generators, each class of models can be represented by a unique contour in the (c d , c u ) plane as a function of the mixing angle. The properties of a variety of benchmark models in the three model classes mentioned above are shown in Table 1. Ref. [42] shows that the finite width of the resonance, which can reach up to 12% for the Z Q , has negligible effects on the translation of experimental limits obtained in the narrow width approximation into the (c d , c u ) plane. Similarly, the effect of interference between the signal and SM backgrounds can be neglected as long as the signal cross section is considered in a sufficiently narrow window around the Z mass [45]. Table 1: Various benchmark models in the GSM [42], LRS [42], and E 6 [1,44] model classes, with their corresponding mixing angles, their branching fraction (B) to dileptons, the c u , c d parameters and their ratio, and the width-to-mass ratio of the associated Z boson. The RS model of extra dimensions [5,47] introduces an additional warped spatial dimension where gravity can propagate through all of the five-dimensional bulk. This leads to the prediction of spin-2 Kaluza-Klein (KK) excitations of the graviton (G KK ) that could be resonantly produced at the LHC and observed in their decay into lepton pairs. These new resonances are characterized by their masses and the coupling k/M Pl , where k is the warp factor of the fivedimensional anti-de Sitter space and M Pl is the reduced Planck mass. The ratios of the intrinsic widths of the first excitation of the graviton to its mass for the coupling parameters k/M Pl of 0.01, 0.05, and 0.10 are 0.01%, 0.36%, and 1.42%, respectively.
To ensure a consistent search strategy and comparability between the results of different searches, the LHC Dark Matter Working Group has defined a variety of simplified models to be used as benchmarks for DM searches at the LHC [48,49]. In this paper, we consider a model that assumes the existence of a single DM particle that interacts with the SM particles through a spin-1 mediator, which can be either a vector or axial-vector boson. The model is fully described by the following parameters: 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. Two signal scenarios have been defined that feature nonzero values of g , so that the mediator could be directly observed in its decay into the dilepton final state. The first choice represents a relative strength of the couplings to quarks and charged leptons typical of some models with a pure vector mediator, while the second choice is a representative case found in the simplest complete models with axial-vector Z bosons [48]: • Vector mediator with small couplings to leptons: g DM = 1.0, g q = 0.1, g = 0.01; • Axial-vector mediator with equal couplings to quark and leptons: g DM = 1.0, Interference between the mediator exchange and the SM processes resulting in the dilepton final state has been studied and found to be negligible for experimental searches [48].

Models with nonresonant signatures
We consider two types of nonresonant excess at high m : a series of virtual spin-2 graviton excitations in the ADD model of large extra dimensions [3], and a four-fermion CI caused by fermion substructure [8]. In these models, the differential cross section for dilepton pair production can be expressed as: where dσ DY /dm is the SM DY differential cross section, η X is a model-specific parameter, and the signal contribution terms are separated into an interference term (I) and a pure signal term (S). Interference between the signal process and the SM DY background is possible when the new process acts on the same initial state and yields the same final state, and can be either constructive or destructive, depending on the sign of η X .
In the ADD model, which is an attempt to describe quantum gravity with an effective field theory, spacetime is extended by n additional compactified spatial dimensions of size L. SM particles are confined to the four-dimensional subspace (the brane), while gravity can propagate to all D = n + 4 dimensions (the bulk). If L is sufficiently large, the D-dimensional fundamental Planck mass M D , which is related to the effective Planck mass M Pl in three dimensions by can then be probed at the TeV scale. In contrast to the RS model, the ADD model predicts the presence of a quasi-continuous spectrum of KK graviton modes as a result of the compactification of the extra dimensions with a large compactification radius. The number of excited modes increases with the interaction scale, resulting in a nonresonant excess at high dilepton masses from the decay of the virtual gravitons, which cannot be individually resolved with the LHC detectors.
These processes can be characterized by the single energy cutoff scale Λ T in the Giudice-Rattazzi-Wells (GRW) convention [50], the string scale M S in the Hewett convention [51], or the number of additional dimensions n in conjunction with M S in the Han-Lykken-Zhang (HLZ) convention [52]. The generic form factor η X in Eq. (2) is replaced by η G , which depends on the chosen convention: Hewett: HLZ: Hereŝ is the square of the parton center-of-mass energy. In the ADD model, interference with DY is limited, as the production of virtual gravitons is dominated by gluon-induced processes. Positive interference is assumed in this analysis, but the choice has negligible effects on the results. The effective field theory will only yield reliable results up to a certain energy scale, above which a more exact theory of quantum gravity is necessary. Both Λ T and M S act as this ultraviolet cutoff parameter in their respective conventions.
In case of the CI model, assuming that quarks and leptons share common constituents, the Lagrangian for the CI process qq → , where is a charged lepton, can be expressed as where q L = (q u , q d ) L is a left-handed quark doublet; q R represents a sum over the righthanded quark singlets (q u -and q d -type); the η ij are real numbers between −1 and 1; and L and R are the left-and right-handed lepton fields, respectively. By convention, g 2 contact /4π = 1 and the helicity parameters η ij are taken to have unit magnitude. The compositeness scale, represented by Λ, is potentially different for each of the individual terms in the Lagrangian. Therefore, the individual helicity currents for "left-left" (LL), "right-right" (RR), and the combination of "left-right" (LR) and "right-left" (RL) in Eq. (7), together with their scales (Λ LL , Λ RR , Λ LR , and Λ RL ), are considered separately in this search, and in each case all other currents are assumed to be zero. While the LL and RR models favor positive values of cos θ * , as does the DY background, the distribution is inverted for the LR and RL models, which favor negative values.
A given η ij can be related to the form factor in the differential cross section in Eq. (2) by

Simulated event samples
Numerous simulated event samples are used to describe both background and signal processes. Dedicated samples are generated for each of the three years of data taking considered in this analysis, reflecting the changing beam and detector conditions from year to year. Different generators are used to generate the individual processes.
The dominant DY background is generated with POWHEG v2 [53][54][55][56][57][58] using next-to-leading order (NLO) matrix elements. In the sample generation, PDFs are evaluated using the LHAPDF library [59][60][61]. For the 2016 samples, the NNPDF3.0 [62] PDF set at NLO is used, while for the samples describing the 2017 and 2018 data, the NNPDF3.1 PDF set computed at the next-tonext-to-leading order (NNLO) [63]  The NLO cross sections obtained from POWHEG are corrected for NNLO effects in perturbative QCD, as well as for missing electroweak effects at NLO, using a correction factor that depends on the dilepton invariant mass, obtained using the FEWZ 3.1.b2 program [67]. In these calculations, the LUXqed plus PDF4LHC15 nnlo 100 [68] PDF set is used. It combines the QCD PDFs, based on the PDF4LHC recommendations [69], with the photon PDFs to account for pure quantum electrodynamics effects. This allows the inclusion of the background arising from a γγ initial state via tand u-channel processes [70,71] in the correction.
The tt, q t W, and WW backgrounds are simulated using POWHEG v2, with parton showering and hadronization described by PYTHIA, using the same PDF sets as for the DY samples. The tt background is normalized to the NNLO cross section calculated with TOP++ [72] assuming a top quark mass of 172.5 GeV, while the q t W simulation is normalized to the cross section computed up to next-to-next-to-leading logarithmic accuracy [73].
In addition to the WW samples, WZ and ZZ events are considered. As the analysis in the dimuon channel for 2016 is unchanged from Ref. [21], samples simulated at leading order (LO) using the PYTHIA program along with the NNPDF2.3 PDFs at LO are used. For the analysis of the 2016 data in the electron channel, as well as the 2017 and 2018 data in both channels, samples are generated at NLO with POWHEG and MADGRAPH5 aMC@NLO version 2.2.2 [74], using the NNPDF3.1 PDF set. The process involving γ * /Z → τ + τ − decays, which is not included in the POWHEG DY samples described above, as well as the W+jets process, are simulated at LO with the MADGRAPH5 aMC@NLO generator. Cross sections for these processes have been calculated up to NNLO with MCFM 6.6 [75][76][77][78].
The NNPDF3.1 NNLO PDF set was found to predict unphysical cross sections at very high dilepton invariant masses (>4 TeV). The DY and tt events generated with POWHEG are therefore reweighted so that the generated m distribution matches the prediction of the NNPDF3.0 NLO set that was used in Ref.
[21], which does not exhibit this undesired behavior. For the DY samples, these weights are applied before they are further corrected to the higher order cross section calculated with FEWZ described above. Thus, the DY background estimate is consistent among the three years of data taking. Other processes also generated with the NNPDF3.1 PDF set are less affected and have only small contributions to the overall background predictions. No reweighting is applied to them.
As the properties of high-mass events containing a spin-1 resonance can be studied by using the DY background samples, only one sample containing a Z SSM boson with a mass of 5000 GeV is generated for illustrative purposes. The RS G KK samples with the graviton generated at different mass values from 250 to 4000 GeV have been generated for all three years. In all samples, the high-mass resonances decay to electron and muon pairs. These signal samples are generated using the PYTHIA  The CI samples, as well as ADD samples for 2017 and 2018, are generated with a lower threshold of 300 GeV on the dilepton invariant mass, while for the ADD samples for 2016, a threshold of 1700 GeV is applied. In case of the CI samples, separate samples for the LL, LR+RL, and RR helicity states are generated, where the PYTHIA implementation of this signal model summed the LR and RL currents into a single process. The generated events in these samples are therefore reweighted to recover the individual LR and RL models. Dedicated PYTHIA DY samples are produced with the same generator settings and subtracted from the signal samples to obtain the respective pure signal yields. In the case of ADD samples, NNLO calculations in QCD show that the correction factor to the LO cross section can be as high as 1.6 [79], and that it always exceeds 1.3 in the considered dilepton mass range. Furthermore, NLO electroweak corrections are not taken into account. This motivates applying the minimal correction factor of 1.3 in the analysis, which also allows a direct comparison to previous results [19].
During the 2016 and 2017 data taking, the CMS L1 trigger was affected by a slowly increasing shift of the reconstructed cluster time in the ECAL, predominantly at high η. This led to a loss of a fraction of events in the trigger, as these clusters could be assigned to the wrong LHC bunch crossing ("trigger prefiring") [80]. As this effect is not present in simulation, simulated events in the dielectron channel are reweighted to account for this inefficiency.
The detector response is simulated using the GEANT4 [81] package. The presence of additional pp interactions in the same or adjacent bunch crossing observed in data (pileup) is incorporated into simulated events by including extra pp interactions that have been generated with PYTHIA. In the dielectron channel, simulated samples are reweighted to reproduce the same pileup distribution as measured in data (pileup reweighting). In the dimuon channel, the trigger, reconstruction, and identification efficiencies exhibit no significant pileup dependence, and no such reweighting is performed.

Lepton reconstruction and event selection
The electron and muon reconstruction algorithms remain close to those used in Ref.
[21], with a few updates implemented. In the electron channel, the 2016 data are reprocessed using an updated electron calibration that corrected for a drift in the energy response at very high energy. For the 2017 and 2018 data sets, there are some modifications to the criteria for both electron and muon identification and electron isolation, as well as for the overall event selection. These improve the treatment of the higher pileup and the selection efficiency for leptons with very high (several TeV) transverse momentum p T .

Electron reconstruction and selection
To select dielectron events in the L1 trigger, events are required to contain at least two electron candidates. The p T thresholds evolved with time, but never exceeded 25 GeV for the electron with the higher p T and 17 GeV for the electron with the lower p T . As a safeguard against cases where one of the electrons fails the L1 trigger, events containing at least one electron, jet, or tau lepton candidate passing the higher p T thresholds at L1 are also considered at HLT level. In the HLT, dielectron events were collected with a trigger requiring two electrons with p T > 33 GeV and |η| < 2.5 in 2016 and 2017. These electrons were also required to pass loose identification criteria using the shower shape in the calorimeter and the matching of a pixel detector track with the energy deposit. The p T threshold was lowered to 25 GeV in 2018, since the electron trigger was improved by making use of the upgraded pixel detector. Secondary trigger paths with progressively higher p T thresholds and fewer selection requirements are employed to monitor the performance of the main trigger.
Electron candidates are reconstructed from energy clusters in the ECAL and tracks from the inner tracker. Typically, electrons lose some of their energy due to bremsstrahlung emission in the tracking detector, which is recovered by adding compatible energy deposits to the ECAL cluster. The final clusters are then associated geometrically with tracks to form electron candidates. As the measurement of the track parameters can be unreliable for high-energy electrons because of bremsstrahlung, the energy of the candidate is taken directly from the ECAL cluster without combination with track information. However, the direction of the electron used in the analysis is still obtained from the combined track and ECAL cluster information. This minor mis-reconstruction has a negligible effect on the dielectron mass measurement.
We select electron candidates with p T > 35 GeV that 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 making up the electron with respect to the nominal center of the CMS detector. The transition region 1.44 < |η C | < 1.57 is excluded since it leads to lower quality reconstructed clusters, owing mainly to inactive material consisting of services and cables exiting between the barrel and endcap calorimeters.
Electrons are required to pass a set of selection criteria optimized for the identification of highenergy electrons [82]. The lateral spread of energy deposits in the ECAL associated with the electron's cluster is required to be consistent with that of a single electron. The electron's track is required to be matched geometrically to the ECAL cluster and to be compatible with originating from the nominal interaction point. Energy deposits in the HCAL in the direction of the electron, corrected for noise and pileup, are required to be less than 5% of the electron's energy. Misreconstructed electrons, and electrons in jets, are suppressed by requiring that the electron be isolated in a cone of radius ∆R = √ (∆η) 2 + (∆φ) 2 = 0.3 in both the calorimeter and tracker [20]. 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, measured in data as discussed below, is at least 95% for barrel (endcap) electrons satisfying p T > 36(37) GeV in 2016, p T > 38(40) GeV in 2017, and p T > 27(29) GeV in 2018.
In the endcap region of the detector, there is a large background of events with hadronic jets misreconstructed as electrons, which comes from W+jets events and from events composed uniquely of jets produced through the strong interaction, referred to as QCD multijet events. To reduce this background, we require that at least one of the electrons be located in the barrel region. The subset of events with both electrons in the endcap region is expected to be overwhelmingly dominated by this background, and is used as a control region for the jet background estimate.
No requirement is made on the sign of electrons because of the increasing probability of charge misidentification for high-energy electrons, which would result in a significant efficiency loss from an opposite-sign selection. The fraction of same-sign events in DY simulation was found to increase from about 1% at the Z boson peak to ≈10% for masses of several TeV. In data, it was found to be consistent with this expectation.
The efficiency to trigger, reconstruct, and select an electron pair with invariant mass m ee = 1 TeV within the detector acceptance is 72 (68)% for 2016 (2017 and 2018) for barrel-barrel events and 67, 60, and 67% for the three years, for barrel-endcap events. It remains stable within a few percent for larger masses. The lower efficiency for barrel-barrel events in 2017 and 2018 is caused by changes to the trigger making use of the upgraded pixel detector to control the trigger rate. The reduction for barrel-endcap events in 2017 is mostly due to trigger prefiring. The total efficiency is estimated using simulated DY events, and the efficiencies of the various selection steps are compared in data and simulation using a method similar to that described in Ref. [83]. Using high-p T Z bosons, it is possible to probe the efficiencies up to p T = 1 (0.7) TeV for electrons in the barrel (endcap) region, with a precision of less than 10%. Based on those measurements, a correction factor is applied to the trigger efficiency in simulation to model the lower p T turn-on observed in data. No correction is found to be necessary for the offline selection.

Muon reconstruction and selection
For the selection of dimuon events, in the L1 trigger, events with at least one single muon with p T > 22 GeV are accepted. At HLT, single-muon triggers requiring p T > 50 GeV and |η| < 2.4 and without an isolation requirement were available during the whole data taking period. During 2017 and 2018, an upgraded version of the muon reconstruction in the HLT was deployed, and as a backup the algorithm used in 2016 was used in an HLT path requiring p T > 100 GeV [84]. The backup trigger increased the overall trigger efficiency for dimuon events with invariant masses above 1 TeV by ≈1%. To collect an event sample unbiased by the trigger turn-on for a data normalization region around the Z boson mass (60 < m < 120 GeV), a prescaled HLT path with a reduced pT threshold of 27 GeV was used. The average prescale value ranged from about 140 in 2016 to about 460 in 2018.
Muon candidates are reconstructed by combining hits in the inner tracking detector with those in the muon system. For muons with p T > 200 GeV, dedicated reconstruction algorithms take into account radiative energy losses in the detector material [85,86]. As no calorimeter information is used in the muon reconstruction, data corresponding to an integrated luminosity of about 3 fb −1 recorded with ECAL or HCAL in a degraded condition are included in the analysis of the dimuon final state. Muons are required to have p T > 53 GeV, |η| < 2.4, and to pass a set of selection criteria optimized for high-p T muons [19]. To reject muons produced inside jets, we require that the scalar p T sum of all tracks identified to come from the same pp collision and within a cone of ∆R = 0.3 around the muon candidate, excluding the muon candidate itself, does not exceed 10% of the p T of the muon.
Dimuon candidates are formed from two muons with invariant mass m µµ > 150 GeV. While the searched-for neutral dilepton final states would decay to truly opposite-sign dileptons, in principle leptons with misreconstructed charge are usable (and are in fact used in the case of electrons, as described above). For muons, misreconstruction of the charge indicates poor resolution of the track fit and thus an unreliable p T measurement. Therefore, the two muons are required to have opposite charge. As the charge misidentification probability is of the order of 10 −4 for muon momenta of 2 TeV [86], this requirement has a negligible effect on the signal efficiency. The muon tracks are fit to a common vertex, which improves mass resolution at high mass. A loose requirement on the quality of the vertex fit rejects some residual background and potentially badly measured signal events (less than 1% in simulation). Backgrounds from cosmic ray muons are reduced to a negligible level by requiring that the three-dimensional angle between the two muons be less than π − 0.02. The remaining contribution of this background after the full selection for events containing muons with p T > 300 GeV is less than 0.2%. The trigger efficiency for dimuon events that pass these offline muon selection criteria has been measured in data using Z boson events. For all three years of data taking, it is above 99.3% for events where both muons are in the barrel region (|η| < 1.2) of the detector and above 99.0% if one or more muons are located in the endcaps (|η| > 1.2). Measured in data using Z events, the overall efficiency to trigger on, reconstruct, and identify muon pairs within the detector acceptance is 93%, with little dependence on the dimuon invariant mass.
Detailed studies of the muon reconstruction performance revealed a slight inefficiency for muons with |η| > 1.6, which increases with muon momentum, caused by electromagnetic showers in the detector material [86]. This effect is more pronounced in data than in simulation. This difference is included in a mass-dependent systematic uncertainty.
The dimuon mass resolution for events with high-p T muons has been studied using highly Lorentz-boosted Z boson events. It was found that for events where at least one muon is in the endcap, the simulation predicted better mass resolution than observed in the data. To correct for this effect, the reconstructed dimuon mass for simulated events in this category is smeared by 15%, bringing data and simulation into good agreement.
The combined product of the acceptance and the efficiency as a function of mass is obtained from simulation for both dielectron and dimuon pairs. This product is shown in Fig. 1, separately for spin-1 and spin-2 particles. The values for the different years of data taking are averaged, weighted by the integrated luminosity of the data sets for each year. Mixtures of polynomials and exponential functions are used to parameterize the mass dependence of the product of the acceptance and the efficiency. The observed turn-on in the low-mass region is a consequence of the lepton p T threshold and the p T -dependent selection efficiency. For m > 200 GeV, the reconstruction and selection efficiencies do not significantly depend on mass, and changes in the product of acceptance and efficiency are due to the changing acceptance as the angular distribution of the leptons changes. For the spin-2 RS graviton, the slight drop at high mass is related to the different production modes (gg and qq) that lead to different angular distribution for the leptons. The fraction of the latter increases with mass and results in a slightly larger fraction of leptons emitted outside the detector acceptance.

Event disambiguation
In a small fraction of the events (<1%), originating mostly from WZ and ZZ production, more than two leptons are present. If several dilepton pairs can be formed in the same event, and if at least one pair has a mass within 20 GeV of the Z boson mass, then the pair closest in mass to the Z boson mass is chosen. Otherwise, the pair with the two highest-p T leptons is selected. For dimuon events, the analysis has not been changed from the published result based on the 2016 data, and the pair with the two highest p T muons is always selected [21]. Events with both a dimuon and dielectron pair that pass the event selection are very rare and both pairs are considered in the analysis.
For both the CI and ADD models, the sensitivity of the analysis can be improved by studying the angle of the outgoing negatively charged lepton with respect to the z axis in the Collins-Soper frame, given by Here p ± represents 1 √ 2 (E ± p z ), and the indices 1 and 2 correspond to the negatively and positively charged leptons, respectively. Events with same-sign electrons are included in the selected event sample, which introduces an ambiguity in the calculation of cos θ * . To resolve the ambiguity, same-sign events in the dielectron channel are assumed to originate from charge misidentification. The measured charge of the electron with the lower value of η C is assumed to be correct and the charge of the other electron is taken to be mismeasured and assigned to be the opposite. This results in the correct lepton being chosen in 76% of same-sign events, according to the simulation.

Background estimation
The dominant (and irreducible) source of background in the signal selection is the DY process. Further sources of high-mass dilepton events are tt and q t W production and diboson production (especially WW events). In addition, both QCD multijet events and W+jets events can be reconstructed as dilepton events if jets or nonisolated leptons are identified as isolated leptons; this category is called "jet misidentification" below.
The DY background shape is estimated directly from the POWHEG samples described above, including the mass-dependent higher-order corrections to the cross section. For the DY → ττ process, the MADGRAPH5 aMC@NLO sample is used. Given the small size of this background, no mass-dependent corrections are applied. The other backgrounds, with the exception of the jet misidentification backgrounds, are also estimated from simulation. The second-largest background contribution originates from tt events and is estimated from the POWHEG samples. The simulation of this process is validated in a data control region containing an eµ pair passing the same lepton requirements as discussed above and is found to describe both the shape of the mass distribution and the event yield observed in data, within uncertainties.
Backgrounds arising from jets misidentified as electrons are estimated from data control regions enriched in QCD multijet events with the methods described in Ref. [19]. Similar methods are used to evaluate the background contributions from jets or nonisolated muons to the sample of selected muons. Even though the uncertainty in these estimates is as high as 50% for the entire mass range, this background has a negligible effect on the results of the analysis, as its contribution to the overall event yield is about 1-3%. The jet background estimate is found to be independent of cos θ * and is therefore scaled by 0.5 when the event sample is split into cos θ * < 0 and ≥0 in the nonresonant interpretation.
Finally, the combined background shape is normalized to the data in the control region around the Z boson mass (60 < m < 120 GeV). For the 2017 and 2018 data, the DY simulation gener-ated with the NNPDF3.1 NNLO PDF set does not accurately describe the Z boson p T spectrum at low p T . The simulated events are therefore reweighted to match the prediction obtained with the NNPDF3.0 NLO PDF, as used for the 2016 samples. For the dimuon channel, a prescaled trigger with a muon p T threshold of 27 GeV was used to collect events in this region. The corresponding offline threshold on the muon p T is 30 GeV. To calculate the normalization factors, the numbers of events in that region are multiplied by the corresponding trigger prescale factor. In the electron channel the normalization factor is 0.95, independent of the year of data taking. In the dimuon channel, these factors are 0.97 for 2016, 1.03 for 2017, and 1.00 for 2018 data.

Systematic uncertainties
Various effects can impact the shape or the normalization of the dilepton invariant mass distribution and lead to systematic uncertainties in the signal and background estimates.
One source of uncertainty describes possible residual systematic effects leading to mismodeling in the simulation of the lepton selection efficiency and momentum measurements at high mass, which are not absorbed in the normalization at the Z boson peak. These uncertainties are estimated from measurements performed with on-shell Z bosons with large transverse boost and measured up to 1 TeV in electron p T and above 2 TeV in muon momentum. High-energy cosmic ray muons are also used to study the muon energy scale [86].
Based on these studies, the overall event selection efficiency including trigger, reconstruction, and identification is assigned a relative uncertainty of 6% for barrel-barrel and 8% for barrelendcap events in the dielectron channel, while an uncertainty of 1-2% is assigned everywhere for the dimuon channel. As discussed in Section 5, a loss of efficiency for very energetic muons caused by electromagnetic showers in the detector material was observed. As this effect is significantly stronger in data compared to simulations, an additional mass-dependent, onesided uncertainty is assigned. Its magnitude is the largest in 2016 data where it reaches −1.0 (−6.5)% at m µµ = 4 TeV for barrel-barrel (barrel-endcap plus endcap-endcap) events. It is −2% or smaller in other years [86]. The difference in the uncertainty between the years is due to the limited size of the data control samples.
The uncertainty in the mass scale is 2 (1)% for dielectron pairs in the barrel-barrel (barrelendcap) category. For dimuon events, potential biases in the mass scale are studied using the generalized endpoint method [86]. The associated uncertainty is 1 (3)% in the barrel-barrel (barrel-endcap plus endcap-endcap) category for 2016, independent of mass. In 2017 and 2018, a mass-dependent uncertainty is assigned that reaches 1-2% at m µµ = 5 TeV.
For dimuon events, an uncertainty in the mass resolution of 15% is assigned for 2016 and 8.5% for 2017 and 2018, based on the studies discussed in Section 5. The difference in the uncertainty between the years is the result of different parameterizations of the mass resolution being chosen. The uncertainty in the dielectron mass resolution is negligible.
The uncertainty in the pileup reweighting procedure in the dielectron channel is evaluated by recalculating the weights, shifting the assumed total inelastic cross section by 4.6% around the nominal value [87]. The resulting changes in the event yield estimates are below 0.6%.
The uncertainty in the trigger prefiring rate in the forward region of the ECAL endcap in 2016 and 2017 is 20%, resulting in a 2-3% uncertainty in the acceptance for dielectron events in which one electron is in the endcap for 2017. Averaged over the full data set, the uncertainty is well below 1%.
Theoretical uncertainties in the DY cross section at high mass related to higher-order corrections and PDF uncertainties grow with mass and reach 20% at m = 6 TeV. Here, the PDF uncertainty is determined with the PDF4LHC procedure [69] using replicas of the NNPDF3.0 PDF set. The contribution from other backgrounds estimated from simulations is mostly relevant for masses below 1-2 TeV, and a 7% uncertainty is applied based on the uncertainties in the cross section calculations of the individual processes.
The estimate for the background from jet misidentification (W+jets, QCD), obtained from data, is assigned a 50% uncertainty in both dimuon and dielectron channels for all years. These estimates are based on the validation of the method in a control region enriched in events with one genuine electron and one misidentified electron in the dielectron channel, and on uncertainties in the measurement of the misidentification rate at high mass in the dimuon channel. The stability of the analysis results against potential further uncertainties in the extrapolation of this background estimate to the highest mass values has been tested and the effect on the results was found to be negligible.
The uncertainty in the Z peak normalization factor in the dielectron channel is 1% in 2016 and 2 (4)% for barrel-barrel (barrel-endcap) events for 2017 and 2018, driven by the typical electron p T in this region being close to the trigger threshold. It is 5% in the dimuon channel for all three years of data taking, arising from the use of a trigger, whose prescale changed with instantaneous luminosity, to collect the Z boson sample.
The uncertainties arising from the limited sizes of the simulated background and signal samples affect the results differently for the resonant and nonresonant signals. In the resonant case they are accounted for in the uncertainty in the parameters in the background parametrization discussed below. For nonresonant signals they are included as a systematic uncertainty in the background and signal yield estimates.
For nonresonant signals, an additional uncertainty is assigned to the reweighting of the 2017 and 2018 samples to match the 2016 PDF set. The size of this uncertainty is mass dependent and reaches 30% at a mass of 5 TeV. The impact of this uncertainty on the limits on the model parameters Λ and Λ T in these signal models is as large as 5%. As no signal simulation is used in the resonant case, it is unaffected by this uncertainty.
In the statistical interpretation for resonant signals, the signal is normalized to the data at the Z boson peak and the backgrounds are normalized to the data in the mass windows around the resonance mass. Uncertainties in the background normalization are taken into account with an overall statistical uncertainty and any remaining uncertainties in the background shape are covered by it. Therefore, only uncertainties in the signal modeling whose impact on the analysis depend on the dilepton mass have to be considered in this case. PDF uncertainties in the signal cross section are considered to be theoretical uncertainties and are not included here. These uncertainties are summarized in Table 2.

Source
Uncertainty Electron selection efficiency 6-8% Muon selection efficiency 1-2% (two-sided), 0-6.5% (one-sided) Mass scale uncertainty 0-3% Dimuon mass resolution uncertainty 8.5-15% For the limits on nonresonant signals and the measurement of R µ + µ − /e + e − , the full set of uncertainties is taken into account. The full set of uncertainties is also taken into account in the uncertainty bands for data to simulation comparisons. The impact of these uncertainties on the background estimate for different mass thresholds is shown in Table 3.

Results
The invariant mass distributions of electron and muon pairs are shown in Fig. 2, combining the 2016-2018 data sets. For illustration, simulated G KK and Z SSM signals with masses close to the exclusion sensitivity of the analysis are shown.
Agreement is observed between the data yields and the expected background, which is quantified in Section 9, although a slight excess of events is seen in the dielectron channel for masses above 1.8 TeV. Careful inspection of those events did not reveal any pathologies. Four events are observed with mass above 3 TeV, two each in the dielectron and dimuon channels. The measured masses are 3.35 and 3.47 TeV in the dielectron channel and 3.07 and 3.34 TeV in the dimuon channel. Table 4 presents the observed and expected number of events for various mass ranges for the dielectron and dimuon pairs. The uncertainty in the total background is calculated taking into account correlations in the uncertainties among the different background sources. No jet background estimate is available in the dimuon channel for 60 < m < 120 GeV, since the data in this control region are collected with the prescaled trigger. However, this background is negligible in this mass range given the large cross section for the Z boson resonance.
The mass distributions are shown again in Fig. 3, split into the two bins in cos θ * used in the search for nonresonant signals, and using the same mass bins as those used in the statistical interpretation for the CI signal, as described in Section 9.2. A CI signal from the LR model is shown to illustrate the improved signal-to-background ratio in the cos θ * < 0 bin for these types of signals. The invariant mass distribution of pairs of (left) electrons and (right) muons observed in data (black dots with statistical uncertainties) and expected from the SM processes (stacked histograms). For the dimuon channel, a prescaled trigger with a p T threshold of 27 GeV was used to collect events in the normalization region (NR) with m µµ < 120 GeV. The corresponding offline threshold is 30 GeV. Events in the signal region (SR) corresponding to masses above 120 GeV are collected using an unprescaled single-muon trigger. The bin width gradually increases with mass. The ratios of the data yields after background subtraction to the expected background yields are shown in the lower plots. The blue shaded band represents the combined statistical and systematic uncertainties in the background. Signal contributions expected from simulated G KK and Z SSM resonances with masses of 3.5 and 5 TeV, respectively, are shown.

Statistical interpretation
Limits are calculated at 95% confidence level (CL) with Bayesian techniques known to have good frequentist coverage properties [88], using the framework developed for statistically combining Higgs boson searches [89], which is based on the ROOSTATS package [90]. For the signal cross section, we use a positive uniform prior, while the nuisance parameters for the uncertainties in dilepton efficiencies, resolution, and scale are modeled with log-normal priors. Two different approaches are followed for the statistical interpretation of the results. To search for resonances in the dilepton invariant mass spectrum, an unbinned maximum likelihood fit is performed to the m spectrum in data, while for nonresonant signals a binned likelihood in m is constructed. In both cases the likelihood fit is done simultaneously for dielectron and dimuon events, years of data taking, and the |η| categories, within the two channels. To further increase the sensitivity to some of the models, the event sample in each of these categories is split into two bins, cos θ * < 0 and cos θ * ≥ 0, in the nonresonant case. The likelihoods for all subcategories are then combined to obtain the results. When the electron and muon channels are combined, we assume that branching fractions of these two signals are the same.

Search for resonant signals
The signal is modeled with the convolution of a Breit-Wigner function to model the intrinsic decay width of the resonance and a double-sided Crystal Ball function [91] to model the mass-dependent mass resolution, except for 2016 data in the muon channel, where a Gaussian function with exponential tails to either side is used [92]. Selection efficiency, mass resolution, and mass scale are taken into account as nuisance parameters in the signal description. The parameters are treated as uncorrelated between the different channels. The relative mass scale between them is the only uncertainty with a noticeable impact. The background is modeled with two different functions for the two dilepton channels. In the dielectron channel it is given by: while for the dimuon channel the following functional form is used: Here m is the dilepton invariant mass; the other parameters do not represent physical quantities. In the dimuon channel, these functions are required to be continuously differentiable at the transition point m threshold between the low-and high-mass parameterizations, which is left free in the fits and ranges from 350 to 750 GeV. In the electron channel, m threshold is set to 600 GeV. They are first fit to the sum of the background estimates, the vast majority of which comes from simulation, as described in Section 6, for masses above 150 GeV. In the limit calculation, the resulting background shapes are then normalized to data within mass windows. The function parameters are treated as nuisance parameters and left floating in the limit calculation. Their initial values are set to those obtained from the fit to the background estimates and they are constrained by log-normal priors whose widths are set to the uncertainties obtained from the same fits. Correlations between the parameters are not considered by these constraints.
The limits are calculated in a mass window of ±4 times the signal width, defined as the sum of intrinsic width and mass resolution, with this window being symmetrically enlarged until there is a minimum of 100 data events in it. This sets an upper limit of 10% on the statistical uncertainty in the local background estimate in the mass window. It is chosen to be large enough to dominate the expected systematic uncertainties in the background shape at high mass, which do not have to be considered explicitly for this reason. To allow the background yield to be constrained by its statistical uncertainties, a log-normal prior with a width of three times these uncertainties is used.
The results have been found to be largely robust against the impact of any remaining uncertainty in the background shape, as well as the choice of background shape parameterization for the full range of considered resonance masses, in the case of a narrow resonance. However, for the largest signal width hypothesis considered, 10%, a significant bias in the limit values for low resonance masses from the choice of the parameterization can occur. Therefore, no results below a resonance mass of 700 GeV are shown for this width.
The parameter of interest is chosen to be the ratio of the cross section for dilepton production via a Z boson to the observed cross section for Z → in the data control region of 60 < m < 120 GeV, R σ : This variable has reduced dependence on the theoretical predictions for the Z cross section and on systematic uncertainties that are correlated between high and low m , such as trigger and selection efficiencies, for which only uncertainties in their stability at high mass remain. To ease comparison with the predicted cross section for various signal models, the limits on R σ are presented multiplied by the theoretical prediction for σ(pp → Z + X → + X) of 1928 pb, calculated at NNLO in QCD and NLO in EWK theory with the FEWZ 3.1 program [67].
The resulting limits for a narrow resonance with an intrinsic width of 0.6% are shown in Fig. 4. At high mass, where the background estimate approaches zero, the limits are independent of the resonance width. The limits are therefore applicable for both Z SSM and Z ψ . The observed limits are consistent with the expectation from SM backgrounds. The expected and observed mass limits in the Z SSM and Z ψ models are shown in   Figure 4: The upper limits at 95% CL on the product of the production cross section and the branching fraction for a spin-1 resonance with a width equal to 0.6% of the resonance mass, relative to the product of the production cross section and the branching fraction of a Z boson, multiplied by the theoretical value of σ(pp → Z + X → + X) of 1928 pb, for (top left) the dielectron channel, (top right) the dimuon channel, and (bottom) their combination. The shaded bands correspond to the 68 and 95% quantiles for the expected limits. Simulated predictions for the spin-1 Z SSM and Z ψ resonances are shown for comparison.
Limits are computed for several hypotheses for the resonance width. In Fig. 5, expected and observed limits are shown for 0.6% (identical to curves in Fig. 4), 3, 5, and 10%. As discussed above, lower thresholds on the resonance mass are chosen to avoid biasing the limits due to the choice of background shape for the wider resonances. The thresholds are 200 GeV for 0.6, 3, and 5%, and, as mentioned earlier, 700 GeV for 10% width. The limit curves for the differ- The limits in the narrow width approximation are interpreted as limits on the generalized couplings in the (c d , c u ) plane as shown in Fig. 6 for the combination of the dielectron and dimuon channels. The closed colored curves show the observed lower limits in the three classes of models (GSM, LRS, and E 6 ), as functions of the mixing angle within each class. For a given point in the plane, a lower mass limit can be obtained with the help of the thin black lines, which are obtained from the observed limit and the ratio of parton luminosities [42]. The observed limit exceeds 4.5 TeV for all models considered and reaches close to 7 TeV in the GSM class of models.
To quantify any discrepancies between the observed data and the background estimates, the p-value for the background-only hypothesis is calculated, as a function of resonance mass hypothesis, following the methodology described in Ref. [93]. The same mass windows as for the limit calculation are used and the resulting p-values therefore depend on the assumed resonance width. Similar to the exclusion limits computed for different width hypotheses, lower mass thresholds are determined based on the potential bias in the fitted number of signal events introduced by the choice of the background parameterization. As this extraction of signal yields is found to be more susceptible to these biases, the resulting thresholds are higher than for the limit setting. The threshold is 200 GeV for a width of 0.6%, 600 GeV for a width of 3%, and 1 TeV for widths of 5 and 10%.
The results are shown in Fig. 7. In the narrow width case, the most significant deviations in the dimuon channel are observed at 520 GeV with a local p-value of 0.0074, corresponding to a local significance of 2.4 standard deviations. The local p-value corresponds to the hypothesis in which a peak at the observed location (520 GeV in this case) is predicted. In both the dielectron channel and the combination, the largest deviation is located at 710 GeV, with local p-values of 0.0010 and 0.0014, respectively, corresponding to local significances of 3.1 and 3.0 standard deviations. By repeating the p-value scan on an ensemble of toy data sets generated from the background parameterizations, a global p-value can also be computed, giving the probability of finding a local p-value as great or greater than that observed, anywhere in a specified mass range. The global p-value for the combination of the two channels is 0.76 in the mass range 500 to 5500 GeV. For such one-tailed tests (for which a p-value of 0.5 corresponds to 0 standard deviations) this formally corresponds to −1.4 standard deviations. Focusing on the vicinity of the observed fluctuation, the global p-value in the mass range 700 to 800 GeV is 0.19, corresponding to 0.9 standard deviations.
In addition to the spin-1 resonances considered thus far, limits are also calculated for spin-  Figure 5: The upper limits at 95% CL on the product of the production cross section and the 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 the production cross section and the branching fraction for a Z boson, multiplied by the theoretical value of σ(pp → Z + X → + X) of 1928 pb, for (upper left) the dielectron channel, (upper right) the dimuon channel, and (lower) their combination. Theoretical predictions for the spin-1 Z SSM and Z ψ resonances are also shown.
2 resonances and the results are shown in Fig. 8. The calculation is performed assuming a narrow width and differs from that for spin-1 resonances only in the detector acceptance for the different spin configurations. The limits in the spin-2 case are more stringent at low mass as the leptons are produced more centrally in the detector, which increases the acceptance, especially in the electron channel where events with two electrons in the endcaps are rejected. Lower mass limits for three example values of the coupling parameter k/M Pl of 0.01, 0.05, and 0.10, corresponding to intrinsic widths of 0.01, 0.36, and 1.42%, are shown in Table 6 and range from slightly over 2 to almost 5 TeV. This improves the previous CMS limits by about 370 to 530 GeV.

GSM
Mixing angle: (-π/2,-π/4) (-π/4,0) (0,π/4) (π/4,π/2) GSM : Q,T 3L ,SM LRS : Y,R,LR,B-L E 6 : η,ψ,N,S,χ The results are also interpreted in the context of the simplified DM models described in Section 3. For this purpose, dilepton production cross sections in these models are calculated at NLO in QCD using the DMSIMP implementation [49] of the simplified model in MAD-GRAPH5 aMC@NLO version 2.6.5 [74]. The partial and total mediator decay widths are included using the MADWIDTH package [94] and are calculated at LO. In case of the axial-vector scenario with coupling g DM = 1.0 and g q = g = 0.1, and assuming m DM > m med /2, the production cross section for electron or muon pairs within detector acceptance ranges from approximately 100 pb for m med around 200 GeV to 0.1 fb for m med around 4 TeV.
In addition to the size of the mediator's coupling to leptons, the sensitivity of this search also depends on m DM relative to m med . For small m DM , the mediator will dominantly decay to DM particles. This decay becomes suppressed for m DM > m med /2, enhancing the decays into leptons and increasing the sensitivity of the dilepton channel. Additionally, the decay width of the mediator is larger if the decay into DM particles is available. Therefore, the width of the mediator changes over the m med -m DM plane. Expected and observed limits are therefore calculated for widths between 0.5 and 3.5% in steps of 0.25%. To determine if a point in the mass plane is excluded, the limit obtained with the width closest to the theoretical value for that point is used. The resulting exclusion contours are shown in Fig. 9

Search for nonresonant signals
To set limits on nonresonant signal models, the parameter of interest in the statistical analysis is the ratio µ of observed to predicted signal cross section. The dilepton mass spectra are divided into multiple exclusive bins. The bin width is optimized separately for the CI and ADD models based on the expected limits. For CI, bins with lower edges of 400, 500, 700, 1100, 1900, and 3500 GeV are used. For the ADD model, the most sensitive part of the invariant mass spectrum, m > 1.8 TeV, is subdivided into 400 GeV-wide search regions, with the final region covering all events above 3 TeV. For the 2016 data, the lowest bin ranges from 1900 to 2200 GeV to ensure that the reconstructed mass distribution is not affected by the lower cutoff on the generated mass of 1700 GeV in the signal samples for that year.
In each bin the background estimate as described in Section 6 and the signal prediction obtained from the samples described in Section 4 are compared. For both signal and backgrounds, all systematic uncertainties described in Section 7 are considered as nuisance parameters modeled with log-normal distributions, taking into account correlations between the different categories where applicable.
The limits for ADD models are obtained using the signal samples produced in the GRW convention. The theoretical description of the model is only valid up to the ultraviolet cutoff pa- rameter Λ T . This is not taken into account in the signal simulation with PYTHIA. However, the signal cross section above the cutoff is small in the Λ T range where limits are set, and has negligible impact on the results of this analysis.
In Fig. 10, the limits on the model parameters in the different ADD conventions are shown in the dielectron, dimuon, and combined channels. In the GRW convention, taking into account the NNLO correction factor of 1.3, the observed (expected) limits are 6.9 (7.2) TeV in the dielectron channel, 7.2 (7.4) TeV in the dimuon channel, and 7.5 (7.8) TeV for the combination. The limits are also translated into the Hewett convention and the HLZ convention for different numbers of extra dimensions using Eqs. (5) and (6). All of the limits are shown in Table 7. The case of n = 2 in the HLZ convention is not considered as there are very strong bounds on M S from gravitational and astrophysical experiments far in excess of the reach of this analysis [98,99]. For the considered models, the limits range from 5.9 to 8.9 TeV, depending on the model. These results are the best to date and improve on the previous most stringent limits by 0.5-1.0 TeV [28]. For the CI model, there can be significant negative interference between the signal and the SM DY process. This can lead to a negative effective signal prediction in some of the mass bins h Ω Figure 9: Summary of upper limits at 95% CL on the masses of the DM particle, which is assumed to be a Dirac fermion, and its associated mediator, in a simplified model of DM production via a (left) vector or (right) axial-vector mediator. The parameter exclusion regions are obtained by comparing the limits on the 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 curves with the hatching represent the excluded regions. The solid gray curves, marked as "Ωh 2 ≥ 0.12", correspond to parameter regions that reproduce the observed DM relic density in the universe [6,48,96,97], with the hatched area indicating the region where the DM relic abundance exceeds the observed value.
or even an overall negative signal contribution for high values of Λ. The interference term is therefore taken into account explicitly in the limit calculation. For the combination of the electron and muon channels, a universal CI is assumed.
The resulting limits for the dielectron, dimuon, and combined channels are shown in Fig. 11. The limits are given separately for the eight distinct CI models. In the electron channel, the observed limit is weaker than the expected limit by up to two σ. This is caused by the excess of electron events discussed in Section 8. The effect is more pronounced here than in the limits in the ADD model as the shape of the mass distribution for the CI model makes it more sensitive to the mass range in which the excess is present. The lower limit on Λ ranges from 23.9 to 36.4 TeV, depending on the model, an improvement over previous CMS results by . Using the same framework of contact interactions, ATLAS has set limits of up to 35. 8 TeV [33]. Signal yields in the LR and RL models are reduced (enhanced) compared to LL and RR in the constructive (destructive) case. However, due to the improved signal-tobackground ratio in the negative cos θ * bin for these models, the reduced sensitivity is mostly recovered in the constructive case, while in the destructive case the increased sensitivity compared to LL and RR is increased even further. The splitting of the event sample into two cos θ * ranges improves the limits by ≈1 TeV in the case of destructive interference and ≈3 TeV in the case of constructive interference.

Search for lepton flavor universality violation
To search for signs of lepton flavor universality violation at high dilepton mass, we consider the ratio R µ + µ − /e + e − of the differential dilepton production cross sections in the muon and electron channels. The reconstructed distributions are distorted compared to particle level by bin-to-bin migration caused by mass scale and resolution effects, and by the detector acceptance and lepton efficiencies.
The mass distributions in data are unfolded after subtracting all backgrounds except for DY. For this purpose, the bin-to-bin migration is quantified as response matrices, i.e. 2D-histograms of the generated versus reconstructed dilepton mass in an event. They are obtained from the DY simulation for events passing the full analysis selection. The response matrices are then inverted and applied to the reconstructed mass spectra to unfold them to particle level using the ROOUNFOLD [100] utility. Since detailed studies found the off-diagonal elements of the response matrices to be small, no regularization is necessary. No significant dependence of the unfolded result on the characteristics of the DY simulation was observed. In the dimuon channel, the size of the unfolding correction rises from around 2% at low mass to around 10% at 3 TeV, caused by the worsening m resolution. The size of the unfolding corrections is 3-5% in the dielectron channel, where the resolution is largely independent of m .
As the interest is in departures from the SM at high mass, we assume, as in the model of Table 7: Exclusion limits at 95% CL for the electron and muon channels, and their combination, for various parameter conventions of the ADD model. Signal model cross sections are calculated up to LO and the NNLO correction factor of 1.3 is applied. For each of the model parameters, the observed limit is shown first, followed by the expected limit in parentheses.
[15], that departures are negligible in the mass region 200-400 GeV. Thus, to correct for the differences in acceptance and efficiency between the dielectron and dimuon channels in this mass region, R µ + µ − /e + e − is normalized to unity here. To account for the remaining mass dependence of the difference between the channels, R µ + µ − /e + e − is measured in simulated DY events after application of the full event selection, the unfolding procedure, and the normalization of the ratio at low mass. The inverse of this ratio is then used to correct R µ + µ − /e + e − in data, which is obtained from the unfolded mass distributions. The size of this correction is up to 5% for events with central leptons and as large as 20% for events with forward leptons.
Several uncertainties, with the most significant originating from the PDFs, cancel in the flavor ratio. Uncertainties in the modeling of the mass resolution and scale are propagated through the unfolding procedure by obtaining response matrices for which the reconstructed mass values are smeared or shifted by the uncertainty. The uncertainties in the unfolded mass spectrum are then calculated as the differences between the mass spectra obtained with the different response matrices. The uncertainties in the detector acceptance and lepton efficiencies are propagated to the final flavor ratio via the acceptance and efficiency corrections.
The resulting ratio R µ + µ − /e + e − as a function of dilepton mass is shown in Fig. 12, for events with both leptons in the barrel region, events with at least one lepton in the endcaps, and their combination. Lepton flavor universality implies that this ratio is unity. Good agreement with this expectation is observed up to 1.5 TeV. At very high masses, the statistical uncertainties are large. Here, some deviations from unity are observed, caused by the slight excess in the dielectron channel discussed above. A χ 2 test for the mass range above 400 GeV is performed. The resulting χ 2 /dof values are 11.2/7 for the events with two barrel leptons, 9.4/7 for those with at least one lepton in the endcaps, and 17.9/7 for the combined distribution. These correspond to one-sided p-values of 0.130 and 0.225, and 0.012, respectively.
As the flavor ratio unfolded to the particle level has been measured, these results can serve as a basis to test models that predict deviations from lepton flavor universality.

Summary
A search for resonant and nonresonant new phenomena in the dilepton invariant mass spectrum in proton-proton collisions at √ s = 13 TeV corresponding to an integrated luminosity of up to 140 fb −1 has been presented. High-mass dielectron and dimuon events were reconstructed and selected with algorithms optimized for electrons and muons with high transverse momenta. Standard model (SM) backgrounds were primarily estimated from simulation, with the dominant Drell-Yan background corrected to the highest order calculations available, including the contribution from photon-induced processes. When searching for resonant signals, the background normalizations were obtained from sidebands in the data, while for nonresonant signals, the background was normalized to the data in a control region around the Z boson peak. No significant deviation from SM expectation is observed.
Upper limits are set on the ratio of the product of the production cross section and the branching fraction in a dilepton channel of a new resonance with an intrinsic width of up to 10% to that of the SM Z boson at 95% confidence level. The limits are interpreted in the context of a sequential SM (SSM) and a superstring-inspired model that predict spin-1 resonances. Lower  Figure 12: Ratio of the differential dilepton production cross section in the dimuon and dielectron channels R µ + µ − /e + e − , as a function of m for (upper left) events with two barrel leptons, (upper right) at least one lepton in the endcaps, and (lower) their combination. The ratio is obtained after correcting the reconstructed mass spectra to particle level. The error bars include both statistical and systematic uncertainties. Two models of nonresonant signatures have been considered. In case of a four-fermion contact interaction, lower limits on the ultraviolet cutoff parameter Λ range from 23.8 to 36.4 TeV depending on the helicity structure of the interaction and the sign of its interference with the SM Drell-Yan background. In the Arkani-Hamed, Dimopoulos, and Dvali model of large extra dimensions, lower limits on the ultraviolet cutoff ranging from 5.9 to 8.9 TeV are set, depending on the parameter convention.
The dimuon and dielectron invariant mass spectra are corrected for the detector effects and, for the first time in this kind of analysis, compared at the TeV scale. No significant deviation from lepton flavor universality is observed.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid and other centers 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, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies:     [35] CMS Collaboration, "CMS luminosity measurement for the 2017 data-taking period at √ s = 13 TeV", CMS Physics Analysis Summary CMS-PAS-LUM-17-004, 2018.