Search for black holes and other new phenomena in high-multiplicity ﬁnal states in proton–proton collisions at √ s = 13 TeV

A search for new physics in energetic, high-multiplicity ﬁnal states has been performed using proton– proton collision data collected with the CMS detector at a center-of-mass energy of 13TeV and corresponding to an integrated luminosity of 2.3fb − 1 . The standard model background, dominated by multijet production, is determined exclusively from control regions in data. No statistically signiﬁcant excess of events is observed. Model-independent limits on the product of the cross section and the acceptance of a new physics signal in these ﬁnal states are set and further interpreted in terms of limits on the production of black holes. Semiclassical black holes and string balls with masses as high as 9.5TeV, and quantum black holes with masses as high as 9.0TeV are excluded by this search in the context of models with extra dimensions, thus signiﬁcantly extending limits set at a center-of-mass energy of 8TeV with the LHC Run 1 data. © 2017 The Author(s). Published B.V. is an open access article the CC license Funded by SCOAP 3 .


Introduction
The standard model (SM) [1][2][3] of particle physics is a remarkably successful theory. However, several outstanding problems remain. One of these is the "hierarchy problem" [4], i.e., the vast separation between the electroweak energy scale and the scale at which gravity becomes strong. The latter, referred to as the "Planck scale" (M Pl ), is some 17 orders of magnitude greater than the former. There are various theoretical extensions of the SM that address the hierarchy problem, such as supersymmetry (SUSY) and models with extra dimensions.
In many of these models, high-multiplicity, energetic final states naturally occur. Strong single or pair production of various new physics signals result in multijet final states, often accompanied by energetic leptons and/or invisible particles resulting in transverse momentum (p T ) imbalance in the event. Examples include a large variety of SUSY signals, both with R-parity [5] conservation [6] and violation [7], and signals associated with technicolor models [8], axigluons [9], colorons [10][11][12][13], and various models with low-scale gravity.
In this Letter, we describe a model-independent search for new physics in high-multiplicity final states, and explicitly test the predictions of two possible solutions to the hierarchy problem. One of these solutions invokes a model with n large extra dimensions, col-E-mail address: cms-publication-committee-chair@cern.ch. loquially known as the "ADD model", named after its proponents, Arkani-Hamed, Dimopoulos, and Dvali [14][15][16]. The other solution is based on the Randall-Sundrum model [17,18], called the "RS1 model". In this model a single, compact extra spatial dimension is warped, and the SM particles are localized on a TeV-scale brane, while gravity originates on the second, Planck brane, separated from the TeV brane in the extra dimension.
In the ADD model, the fundamental multidimensional Planck scale (M D ) is related to the "apparent" 3-dimensional Planck scale M Pl as: where r is the compactification radius or the characteristic size of extra dimensions.
In the RS1 model, the analog of the ADD scale M D is defined as a function of the exponential warp factor k and the compactification radius r: In both models, M D can be of order 1 TeV, thus eliminating the hierarchy of scales and alleviating the hierarchy problem.
At high-energy hadron colliders, such as the CERN LHC, if the collision energy exceeds M D , both the ADD and RS1 models allow for the formation of microscopic black holes (BHs) [19][20][21][22][23]. In the simplest scenario, microscopic BHs are produced when the distance of closest approach between two colliding particles is less than the Schwarzschild radius R S , which for a BH in a (3 + n)-dimensional space is given by [24]: where M BH is the mass of the BH. The parton-level production cross section of such processes is expected to be simply π R 2 S [19,20]. In more complicated scenarios with energy loss during the formation of the BH horizon, the production cross section could significantly depart from this simple geometrical formula. Since the production of BHs is a threshold phenomenon, we assume a minimum mass threshold M min In the semiclassical case, corresponding to M BH M D , BHs evaporate rapidly via Hawking radiation [25], with a lifetime of order 10 −27 seconds. As gravity couples universally to the energymomentum tensor and does not distinguish between various particle species, microscopic BHs decay democratically into all SM degrees of freedom, i.e., all SM particle species with all possible values of quantum numbers, such as spin, color, and charge. The final state is therefore populated by a variety of energetic particles, such as hadrons (jets), leptons, photons, and neutrinos. Due to the large number of color degrees of freedom, about 75% of the particles produced are expected to be quarks and gluons. The final state may contain significant transverse momentum imbalance from the presence of neutrinos, which constitute about 5% of the decay products. Other processes, such as the decay of W and Z bosons, or of heavy-flavor quarks, and the possible emission of gravitons or the formation of a noninteracting stable BH remnant, contribute to the transverse momentum imbalance as well.
The semiclassical approximation breaks down when the mass of the BH approaches M D and the BH becomes a quantum object or a quantum black hole (QBH). These objects do not obey the usual BH thermodynamics and hence decay much more rapidly than their semiclassical counterparts. Their decays are characterized by the presence of only a few particles, e.g., a pair of jets [26][27][28]. These QBHs could also decay into lepton flavor violating final states, as preserving baryon number or lepton numbers separately is typically not a requirement of the decay process [26,27].
In addition to semiclassical BHs and QBHs, one could also explore stringy precursors of BHs, called "string balls" (SBs) [29]. Such objects, which arise in string theory, are highly excited, long, folded strings that form below the BH production threshold. Like semiclassical BHs, SBs evaporate thermally, but at a constant Hagedorn temperature [30] independent of the SB mass, and also produce a large number of energetic particles in the final state, with the composition similar to that for a semiclassical BH. String balls undergo a phase transition into ordinary semiclassical BHs when their mass reaches M S /g 2 S [29], where M S and g S are the string scale and the string coupling constant, respectively. For an SB mass between M S /g S and the BH transition, M S /g 2 S , the parton-level cross section saturates at σ ∼ 1/M 2 S , while for lighter SBs, it grows [29]. While our choice of final states is inspired by the production of microscopic BHs, in this Letter we focus on a generic search (Section 8) that can be used to probe a large class of new-physics models. Consequently, our emphasis is on the exploration of a multiparticle final state with a model-independent search.
During Run 1 of the LHC, a number of searches for semiclassical and quantum BHs were performed at a center-of-mass energy of 8 TeV. A review of these results can be found in Ref. [31]. The limits on the minimum BH mass set by these searches lie in the 6 TeV range. With the increased LHC center-of-mass energy of 13 TeV, the BH phase space can be probed much more extensively, as was demonstrated in the recent ATLAS publications [32][33][34], which set BH mass limits reaching 9 TeV.

Analysis strategy
The most challenging aspect of the analysis presented in this Letter is accurately describing the QCD multijet background, since the BH signal leads to a broad excess in the S T spectrum, rather than a narrow peak. Here, S T is defined as the scalar sum of the transverse energies of jets, leptons, photons, as well as the missing transverse energy (E miss T , defined as the magnitude of the transverse momentum imbalance in an event, as detailed in Section 4): (4) where N is the total number of final-state objects (excluding the E miss T ), or the object multiplicity. For the QCD background, the final-state objects are almost exclusively jets and the E miss T is expected to be small, so the S T variable is reduced to a scalar sum of the transverse momenta of the jets. The signal region for this search typically lies in the high-multiplicity regime where the QCD multijet background is dominated by higher-order effects. These effects have not been calculated for high-multiplicity final states, and therefore an accurate simulation of the QCD multijet background, pertinent to our signal region, does not yet exist. This significant hurdle is mitigated by predicting the QCD multijet background directly from data using a new technique developed in Run 1 of the LHC [35][36][37]. Studies performed with simulated QCD multijet events and with data at low object multiplicities show that the shape of the S T distribution above its turn-on threshold is approximately independent of the multiplicity of the final state. This observation is consistent with the development of the parton shower via nearly collinear emission, which approximately conserves the S T value, up to the effects of additional jets falling below the kinematic threshold. For this reason one can predict the S T spectrum of a multijet final state using samples of dijet or trijet events. This feature provides a powerful tool to predict the shape of the S T spectrum at higher multiplicity using a low-multiplicity control region. The method has found wide applicability in various CMS searches, such as a search for stealth SUSY [38] and a search for multijet resonances [39]. An earlier CMS analysis [35] also considered other kinematic variables, such as the invariant mass or transverse invariant mass of the event. However, the multiplicity invariance is not exhibited by these variable to the degree shown by the S T variable.
In this Letter we follow closely the methodology of Refs. [35][36][37] geared toward a multiparticle final state, dominated by QCD multijets in the case of semiclassical BHs, and toward a dijet final state for the QBHs. The variable S T is the single discriminating variable used in the analysis, chosen for its robustness against variations in the BH evaporation model and its lack of sensitivity to the relative abundance of various particles produced. This variable encompasses the total transverse energy in an event and is therefore useful in discriminating between the signal and the background. There is a minimum transverse energy (E T ) threshold of 50 GeV [35] that each of the objects (including the E miss T ) has to satisfy to be counted toward the definition of S T . The exact choice of the E T threshold is not particularly important; the 50 GeV threshold is chosen as it makes the analysis insensitive to additional interactions in the same or adjacent bunch crossings (pileup) and moderates the effect of initial-state radiation, which generally spoils the S T invariance.
The S T distributions are produced in a number of inclusive object multiplicity bins (N ≥ N min ). The background is estimated exclusively from collision data via the appropriately chosen control regions. This approach does not rely on the Monte Carlo (MC) description of the backgrounds. In addition, this method has the advantage of being more sensitive and less model dependent than exclusive searches in specific final states, e.g., the lepton + jets final state [40,41]. The ATLAS Collaboration has recently adopted a similar inclusive search strategy based on the 2012 and the 2015 √ s = 13 TeV data [33].

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3. ECAL and HCAL cells are summed to define the calorimeter tower energies, subsequently used to provide the energies and directions of hadronic jets. The first level of the CMS trigger system [44], composed of custom hardware processors, combines information from the calorimeters and muon detectors to select the most interesting events in a fixed time interval of 3.2 μs. The high-level trigger (HLT) processor farm further decreases the event rate from around 100 kHz to less than 1 kHz, before data storage.
A 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. [45].

Event reconstruction
The analysis is based on proton-proton collision data corresponding to an integrated luminosity of 2.3 fb −1 , collected with the CMS detector in 2015. The trigger chosen for the analysis is based on the total transverse energy of jets in an event. At the HLT, events are selected if they passed an 800 GeV threshold on the scalar sum of the transverse momenta of all hadronic jets, which are reconstructed with the particle-flow (PF) algorithm [46], as described below. The trigger is nearly 100% efficient for S T above 1 TeV.
In the subsequent analysis, we select events with at least one reconstructed vertex [43] within 24 (2) cm of the nominal interaction point measured parallel (perpendicular) to the LHC beamline. The vertex with the highest sum of the transverse momenta squared of the associated tracks is chosen as the hard-scattering vertex.
The reconstruction of physics objects in the event is based on the PF algorithm that identifies each single particle in an event (photon, electron, muon, charged hadron, neutral hadron) using an optimized combination of all subdetector information. The energy of photons is directly obtained from the ECAL measurement, corrected for zero-suppression effects [47]. The energy of electrons is determined from a combination of the track momentum at the main interaction vertex, the corresponding ECAL cluster energy, and the energy sum of all bremsstrahlung photons attached to the track. The energy of muons is obtained from the corresponding track momentum. The energy of charged hadrons is determined from a combination of the track momentum and the corresponding ECAL and HCAL energy, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energy. Finally, the E miss T is defined as the absolute value of the vectorial p T sum of all the PF candidates reconstructed in an event [48].
For each event, hadronic jets are clustered from the PF candidates using the anti-k T algorithm [49] with a distance parameter of 0.4, as implemented in the FastJet package [50]. The jet momentum is determined as the vectorial sum of the momenta of all its constituents, and is found to be within 5 to 10% of the true (particle-level) momentum over the whole p T spectrum and the detector acceptance. Jet energy corrections are derived from simulation and in situ measurements of the energy balance in dijet, multijet, γ + jet, and leptonic Z + jet events [51,52]. The contribution of charged hadrons that do not originate from the hard-scattering vertex ("pileup"), but are clustered during the reconstruction of a jet, is subtracted. Corrections based on the jet area [53] are applied to remove the energy contribution of neutral hadrons from pileup interactions. The jet energy resolution (JER) is approximately 8% at 100 GeV and 4% at 1 TeV [51,52]. The minimum threshold on the corrected E T of the jets used in the analysis is 50 GeV, and jets are accepted in the full pseudorapidity range of the CMS calorimeters (|η| < 5). To reduce contamination from poorly reconstructed muons, the fraction of the jet momentum carried by a muon is required to be less than 80%. Also, to suppress jets due to rare, spurious anomalous calorimeter signals, jets must contain at least two particles, one of which is a charged hadron, and the jet energy fraction carried by neutral PF candidates (neutral hadrons and photons) should be less than 90%. These criteria have an efficiency greater than 99% per jet.
Details of muon reconstruction can be found in Ref. [54]. Muon candidates are required to satisfy a minimum p T threshold of 50 GeV and to be within |η| < 2.4. The transverse impact parameter and the longitudinal distance of the track associated with the muon with respect to the primary vertex is required to be less than 2 and 5 mm, respectively, to reduce contamination from cosmic muons. The muon candidate is required to have at least one energy deposit in the pixel tracker and at least six deposits in the silicon strip tracker. The global track fit to the tracker trajectory and to at least two muon detector segments must have a χ 2 per degree of freedom of less than 10.
Details of electron reconstruction can be found in Ref. [55].
Electron candidates are required to have E T > 50 GeV, |η| < 2.5, excluding the 1.44 < |η| < 1.57 transition region between the ECAL barrel and endcap detectors where the reconstruction is suboptimal, and to pass standard identification criteria, corresponding to a working point with an average efficiency of 80% per electron. Both muons and electrons are required to be isolated from other energy deposits in the tracker and the calorimeters. The relative isolation I is defined as the ratio of the sum of transverse Table 1 Generator settings for various semiclassical BH model points probed in this analysis. These parameters are defined in Refs. [56] and [57] for the BlackMax and charybdis 2 generators, respectively. The generator settings not specified are kept at their default values.

Model
Choose_a_case Mass_loss_factor Momentum_loss_factor turn_on_graviton bhspin  mjlost  yrcsc  nbodyaverage  nbodyphase  nbodyvar  rmstab   C1  TRUE  FALSE  FALSE  FALSE  TRUE  TRUE  FALSE  C2  FALSE  FALSE  FALSE  FALSE  TRUE  TRUE  FALSE  C3  TRUE  FALSE  FALSE  TRUE  FALSE  FALSE  FALSE  C4  TRUE  TRUE  TRUE  FALSE  TRUE  TRUE  FALSE  C5  TRUE  TRUE  TRUE  FALSE  FALSE  FALSE  TRUE energies of photons, and charged and neutral hadrons in a cone of radius of 0.4 (0.3) in the η-φ space centered on the muon (electron) candidate to the p T of the lepton: where the sum runs over charged hadrons originating from the hard-scattering vertex, neutral hadrons, and photons. The numerator of the ratio is corrected for contributions due to pileup (E PU T ), using the fraction of energy carried by the charged hadrons originating from other vertices to estimate the contribution of neutral particles from pileup for muons, and an average area method [53], as estimated with FastJet for electrons. Muons (electrons) are required to have relative isolation values less than 0.15 (0.10).
Photons are required to have p T > 50 GeV and to be in the same fiducial region as the electrons, and are reconstructed with a standard algorithm [47] using a medium working point. They are also required to be isolated, with the definition of isolation tuned to yield constant efficiency as a function of photon p T . The pileup-corrected charged-hadron transverse energy sum within the isolation cone of radius of 0.3 is required to be less than 1.37 (1.10) GeV in the barrel (endcaps). A similarly defined neutral-hadron isolation sum is required to be less than 1.06 GeV + 0.014p T + 0.000019 GeV −1 p 2 T in the barrel and 2.69 GeV + 0.0139p T + 0.000025 GeV −1 p 2 T in the endcaps. Finally, the pileup-corrected isolation sum of any additional photon candidates in the isolation cone must be less than 0.28 GeV + 0.0053p T (0.39 GeV + 0.0034p T ) in the barrel (endcaps).
In the semiclassical case, several models are explored. We simulate the following scenarios: nonrotating BHs (models B1, C2), rotating BHs without energy loss (models B2, C1) and with an alternative evaporation model (C3), rotating BHs with 10% loss of mass and angular momentum (B3), rotating BHs with Yoshino-Rychkov bounds [60] (model C4), and rotating BHs with a stable remnant with the mass equal to the multidimensional Planck scale M D (model C5, for which additionally the Charybdis2 parameter NBODY was changed from its default value of 2 to 0). The generator parameters used for semiclassical BH signal models are given in Table 1.
The parameters associated with the BH signal generation, as can be inferred from Eq. String balls are generated using the charybdis 2 event generator for the nonrotating scenario. The number of extra dimensions is fixed at n = 6, as it was in earlier CMS publications. Note also that as the dependence of the SB dynamics on n is only implicit, via the relationship between the Planck and string scales. The mass of the string ball (M SB ) is varied between 5 and 10 TeV. Four different benchmark scenarios are considered: These benchmarks are chosen such that the various regimes of the SB dynamics are probed, mainly between the first transition at M SB = M S /g S and the BH transition at M SB = M S /g 2 S . Above the second transition, SBs turn into semiclassical BHs, so the standard BH analysis fully applies. (We note that a significant fraction of the SB parameter space probed by the ATLAS Collaboration [33] in fact falls into the BH, and not the SB regime.) In all cases the fundamental Planck scale M D is chosen to satisfy the matching condition at the SB/BH transition point: where σ BH (M BH ) and σ SB (M SB ) are the BH and SB parton-level production cross sections, as functions of their masses, respectively.
The choice of parton distribution functions (PDFs) used in this analysis was made as a result of detailed studies of the PDF dependence of signal cross sections. The leading order (LO) MSTW2008LO [61,62] PDFs are chosen for all the signal samples. This PDF set provides a conservative estimate of the signal cross section at high BH masses with respect to the modern NNPDF3.0 [63] PDFs. The cross section at large BH masses obtained with different PDF sets can vary up to an order of magnitude. The use of the MSTW2008LO PDF corresponds to the cross section at the lower edge of the uncertainty range for the NNPDF3.0 PDFs, therefore indicating a reasonable and conservative choice for signal simulation. The MSTW2008LO PDF has been recently superseded by the MMHT2014LO [64] set. However, a numerical comparison of the signal cross sections computed using both of these PDF sets reveals no differences. This is expected, as no new data constrain-ing PDFs for processes with such high momentum transfer have been utilized in the global fit used to extract the MMHT2014LO PDF set. Given that MSTW2008LO was the PDF of choice for the Run 1 version of this analysis, the use of this PDF also allows one to directly compare Run 1 and Run 2 results.
The CMS detector simulation is performed using both detailed simulation via Geant4 [65] (semiclassical BlackMax and charybdis 2 C1-C3 BH samples) and fast parametric simulation via the FastSim [66] package (QBH and SB samples, as well as charybdis 2 C4-C5 samples). The fast simulation samples are validated against the full simulation, and the small observed differences between the two approaches were found to have a negligible impact on the signal acceptance.
In addition we use simulated samples of W + jets, Z + jets, γ + jets, tt, and QCD multijet events for auxiliary studies. These events are generated with the MadGraph5_amc@nlo [67] event generator, followed by pythia for description of hadronization and fragmentation. The NNPDF3.0 PDF set with the tune CUETP8M1 is used for the background generation, and the CMS detector response is simulated via Geant4.

Backgrounds
The main SM backgrounds to the multiparticle final states from new-physics processes are QCD multijet and γ + jets production, vector bosons produced in association with jets (V + jets, where V stands for a W or Z boson), and tt process. The QCD multijet background is by far the dominant one. The additional backgrounds are not explicitly considered for the remainder of the analysis, as they together contribute only a few percent of the total background. They also have a very similar shape to the QCD multijet background, as evident from Fig. 1, which shows the comparison in the S T distribution of the total contribution of all the simulated backgrounds and the data for both the low-multiplicity N = 2 and high-multiplicity N ≥ 6 selections. While simulated samples are not used to estimate the background in this analysis, we note that the simulation nevertheless describes data reasonably well.
The background estimation strategy hinges on the shape of the S T distribution in the tail region of the spectrum being independent of the object multiplicity. As a consequence of this invariance, the S T distribution in higher multiplicity regions can be obtained by appropriately rescaling the low-multiplicity S T spectrum, where signal contamination is expected to be minimal. Given the similarity of the shapes of the QCD multijet and other backgrounds, this technique also implicitly accounts for most of the contribution of the non-QCD backgrounds. The technique has been extensively validated in previous CMS publications [35][36][37], using both lowmultiplicity data samples and simulated QCD multijet events.
The shape of the background is obtained by performing a binned likelihood fit of the S T spectrum to a functional form given , where the P i are free parameters of the fit, in the range between 1.4 and 2.4 TeV. The functional form used to fit the background is taken from earlier searches at the LHC and at the Fermilab Tevatron [68], and represents a wellestablished characterization of rapidly falling E T spectra. Since the goal is to operate in the background-only control region for the extraction of these shapes or templates, the multiplicities chosen for the background extraction are N = 2 and 3. This choice is justified by the dedicated Run 2 analyses (e.g., [69]) as well as by the earlier iterations of this analysis (e.g., [70]) where no signal of new physics was observed in these low-multiplicity regions.
The background template extracted by fitting the S T spectrum corresponding to the exclusive object multiplicity of N = 2 is normalized appropriately to obtain an estimate of the background for the S T spectra at higher multiplicities. The normalization regions in S T depend on the object multiplicity. Their definition is based on the studies of the S T invariance in simulated QCD multijet samples. The lower S T bound of the normalization region is chosen such that it is above the S T turn-on region, while the choice of the upper S T bound is guided by the need to operate in a regime of low signal contamination, namely where one does not expect a significant event yield for signals that have cross sections below those already excluded by the previous CMS search [37]. The Table 2 The normalization regions and the corresponding N = 2 background template normalization factors s and their uncertainties for inclusive multiplicities, N ≥ 2 . . . 10 Table 2.
To ascertain the uncertainties associated with this method of background extraction, two additional fitting functions are considered, given by P 0 (P 1 + x) −P 2 and P 0 (P 1 + P 2 x + x 2 ) −P 3 . With an aim to compute a conservative estimate of the uncertainty, both the N = 2 and 3 S T spectra are used to estimate the background shape. This leads to six different functions used in the fit: one nominal and five additional ones that form an uncertainty envelope that is symmetrized around the main fit. The shape systematic uncertainty is indicated with the gray shaded bands in Figs. 2 and 3, which show the S T spectrum observed in data for various inclusive multiplicities and the background predictions with their uncertainties. It ranges between 1 and 200% and rapidly increases in the high-S T range because of the limited number of events in the N = 2 and 3 S T spectra used to derive the background templates. We assign an additional 5% systematic uncertainty to address possible deviation from the S T invariance by estimating the difference between the background predictions based on the N = 2 (default) and N = 3 templates. This is a subdominant uncertainty in the high-S T range relevant for the analysis. The systematic uncertainties in the background estimate are detailed in Table 3 of Section 7.

±5%
As seen from Figs. 2 and 3, the data agree well with the predicted background and no evidence for new physics production is observed in any of the multiparticle final states studied.

Systematic uncertainties
A detailed study of the systematic uncertainties associated with this analysis was carried out, and these uncertainties were taken into account while computing exclusion limits, as detailed in Section 8. The dominant sources of systematic uncertainties arise from the effects of the jet energy scale (JES) and JER uncertainties, from the choice of the PDFs, and from the migration between event categories as a result of final-state radiation (FSR). The uncertainties vary for different signal samples and different inclusive multiplicities. The following sections give details on their estimation and the typical range for each uncertainty.

Parton distribution functions
There are two different sources of the PDF uncertainty that could affect the signal acceptance. The uncertainty could simply come from the choice of the PDF set. Alternatively, the systematic uncertainty could come from the variations induced by the uncertainties associated with the fit parameters in a chosen PDF set. The latter uncertainty is estimated by calculating 2k + 1 weights per event, where k is the number of alternative PDF sets with parameters varied within their associated uncertainties either up or down. These alternative sets are provided by the authors of each PDF set, along with the main PDFs used in the signal generation. The signal acceptance is then calculated for each of these weights, resulting in 2k + 1 values of the acceptance. For a particular choice of PDF set, the systematic uncertainty in the acceptance is quoted as the sum in quadrature of the deviations from the central value for each of the k PDF eigenvectors. This analysis is further extended by using various PDFs, such as MSTW2008LO, CTEQ6.1L [71], and CT10 [72]. The uncertainty is computed for a particular benchmark scenario, a nonrotating BH with M D = 3 TeV, M BH = 5.5 TeV, and n = 2, representative of the uncertainties for other benchmark points in the range probed by this analysis. The uncertainty in the acceptance using the variation of the chosen PDF, MSTW2008LO, does not exceed 0.5%. The uncertainty due to the choice of the PDF is significantly larger and can be as high as 6%. Following the recommendations of the PDF4LHC group [73,74], we assign a total uncertainty of 6% associated with this source of the systematic uncertainty.

Jet energy scale and resolution
The CMS experiment has adopted a factorized approach [51,52] to the application of corrections associated with the energy of a jet. After the subtraction of the additional energy due to pileup, the energies of the reconstructed jets are corrected for the nonlinear response of the calorimeter. These corrections are parametrized in p T and η, and are derived from simulation and in situ measurements of the energy balance in dijet, multijet, γ + jet, and leptonic Z + jet events [51,52]. Given the predominance of jets in the final state of interest, the uncertainty related to the JES can significantly affect the signal acceptance. The JES is modified by one standard deviation in both upward and downward directions, and the corresponding changes in the jet energies are propagated into the E miss T and S T computation. For various model-specific scenarios, we find the variation to lie between 1 and 5% and conservatively assign a uniform uncertainty in signal acceptance of 5% due to the JES uncertainty.
The uncertainty due to the JER is estimated by oversmearing jet energies in simulated signal samples with η-dependent factors to match the resolution observed in data. The effect of this change to the jet energy is propagated to the E miss T and S T calculations. The JER uncertainty in the signal acceptance varies between 1.2 and 4.0%; we therefore assign a conservative uniform uncertainty of 4% to account for this effect.
In estimating both the JES and JER uncertainties, we find that there are migrations of events in different multiplicity categories due to the jet p T values moving either above or below the 50 GeV threshold. The respective values of uncertainties associated with both of these sources sufficiently cover the effect of event migration. These two uncertainties affect only the simulated signal yields, as the background determination does not rely on simulation.

Final-state radiation
The presence of imperfectly modeled FSR may result in event migration between various multiplicity bins due to the change in the number of objects counted toward the S T . The uncertainty in the signal acceptance due to the imperfect modeling of FSR is estimated by varying the parameters in the pythia 8 generator that govern FSR modeling, and is found to be 1.2%.
The sources of systematic uncertainties considered in the analysis and their magnitude are listed in Table 3.
In addition, an uncertainty in the integrated luminosity of 2.7% [75] is propagated to both the model-independent and modelspecific limits. All other uncertainties are negligible [35][36][37].

Model-independent limits
No statistically significant deviations of data relative to the background predictions are observed in any of the spectra. Exclusion limits on various signals are set using the modified frequentist CL s approach [76,77] with the profile likelihood as a test statistic and nuisance parameters implemented via log-normal priors. The likelihood of the observed number of events is modeled with a Poisson distribution. The limit calculation is performed using the methodology developed by the ATLAS and CMS Collaborations in the search for the Higgs boson [78]. The systematic uncertainties taken into account are detailed in Section 7.
As mentioned in Section 1, the main emphasis of the analysis is the computation of model-independent limits on the product of the hypothetical signal cross section and acceptance (σ A) in inclusive N ≥ N min multiplicity bins, as a function of the minimum S T requirement (S min T ). These limits represent the minimum σ A of a potential signal excluded by the analysis and do not assume any particular signal model. This makes the limits applicable to a large variety of BH models and other theoretical models with processes that result in high-multiplicity final states. To obtain the model-independent limits, a counting experiment is performed for Model-independent limits can be used straightforwardly to test any model of new physics that predicts signal corresponding to high-multiplicity, energetic final states. Given that the object selection efficiency is close to 100%, all is needed to test a particular model against the results of this search is to estimate the signal acceptance as a function of the minimum S T requirement, corresponding to the selections described in Section 4, which can be done with sufficient precision even at the generator level. A comparison of the product of the signal cross section and the acceptance as a function of the minimum S T requirement with the model-independent limit at the relevant object multiplicity would then indicate whether a particular model has been excluded by this analysis.

Model-specific limits
In the case of the model-specific limits, the selection criteria (N min , S min T ) were chosen for optimum sensitivity to each particular model. This was done by using the lowest value of the expected limit, as well as the maximum value of the Z bi [79] statistic, and converting the corresponding model-independent limit into a limit on the BH or SB mass, using the known signal cross section and acceptance. The Z bi statistic is a measure of equivalent Gaussian signal significance obtained by considering the binomial probability of the events in data being distributed at least as signal-like as observed, under the assumption of the background-only hypothesis. In the majority of cases, the limit-and significance-based optimizations are in agreement and lead to the same set of optimum selection requirements. Since the signal MC benchmark points were produced on a grid of discrete M BH values (with a typical spacing of 1 TeV), in order to obtain the mass limit, we smoothly interpolate the signal cross section times acceptance and the exclusion limit between the adjacent points on the grid. Typical precision of this procedure expressed as a limit on the minimum BH mass is ∼0.1 TeV, although for some of the points where the grid spacing was not as fine, the precision can be somewhat worse.
Model-specific limits span the entire set of models discussed in Section 5. The limits on semiclassical BHs in all cases, except for the model with stable remnant (C5), come from the optimized values of N min equal to 9 or 10 for low BH masses and from lower   values of N min for high BH masses, thanks to smaller backgrounds at high values of S T . In Fig. 6, we show the observed lower limits on the semiclassical BH mass at 95% CL for BlackMax signal samples. In these models, we exclude minimum BHs masses below 7.0-9.5 TeV, for a large set of model parameters. Similar limits for models generated with the charybdis 2 generator are shown in Fig. 7. Considering that an independent optimization was used for each model, and for each M D and n point, these limits are in general agreement with those obtained using BlackMax for analogous models. The limits significantly extend those from LHC Run 1 [37], which only reach 6.5 TeV.
In the case of QBHs and the case of semiclassical BHs with a stable remnant, where the number of produced particles is small as the massive stable remnant escapes detection, the best sensitivity is achieved for N min = 2. This sample partially overlaps with the N = 2 data set used for deriving the background template. Nevertheless, even in this case we still can use the background estimate based on the N = 2 spectrum, because the background is determined using only the low-S T range of the spectrum  (1.4-2.4 TeV), where the signal has already been excluded up to ≈5 TeV by the 8 TeV analysis [37]. Thus the potential signal contamination of the S T range used in deriving the background shape and normalization is negligible.
The lower limits on the minimum QBH mass are shown in Fig. 8 and span the 7.3-9.0 TeV range for the ADD (n > 2) and 5.1-6.2 TeV range for the RS1 (n = 1) QBHs. Again, these limits extend significantly those obtained in the LHC Run 1 (4.6-6.2 TeV). Finally, for the case of the SBs, the mass exclusion limits shown in Fig. 9 reach 8.0-8.5 TeV. Most of these limits correspond to the saturated string ball regime, M S /g S < M SB < M S /g 2 S . The transitions between the saturated SB and BH regimes are not clearly visible in theoretical cross section curves, as the parton-level cross section that exhibits this transition as a change in slope is significantly modified by the rapidly falling PDFs.

Summary
We have conducted a search for new physics in multiparticle final states in a data sample of proton-proton collisions at √ s = 13 TeV collected with the CMS detector, corresponding to an integrated luminosity of 2.3 fb −1 . The discriminating variable between signal and the dominant QCD multijet background is the scalar sum of the transverse energies of all reconstructed objects in the event, S T . The shape of the S T distribution in low-multiplicity data is used to predict the QCD multijet background in high-multiplicity signal regions. No significant excess of events over the standard model expectation is observed in any of the analyzed final-state multiplicities. Comparing the S T distribution in data with that from the background prediction, we set model-independent upper limits at 95% confidence level on the product of the cross section and the acceptance for hypothetical signals. In addition, we set limits on various theoretical black hole and string ball models, including models of rotating and nonrotating black holes and quantum black holes. In all cases the exclusions represent significant improvements over the limits achieved in Run 1 of the LHC.

Acknowledgements
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 for delivering so effectively the computing infrastructure essential to our analyses.