Simplified Models for Exotic BSM Searches

Simplified models are a successful way of interpreting current LHC searches for models beyond the standard model (BSM). So far simplified models have focused on topologies featuring a missing transverse energy (MET) signature. However, in some BSM theories other, more exotic, signatures occur. If a charged particle becomes long-lived on collider time scales - as it is the case in parts of the SUSY parameter space - it leads to a very distinct signature. We present an extension of the computer package SModelS which includes simplified models for heavy stable charged particles (HSCP). As a physical application we investigate the CMSSM stau co-annihilation strip containing long-lived staus, which presents a potential solution to the Lithium problem. Applying both MET and HSCP constraints we show that, for low values of $\tan\beta$, all this region of parameter space either violates Dark Matter constraints or is excluded by LHC searches.


Introduction
In theories beyond the standard model (BSM) with an unbroken Z 2 -symmetry the lightest Z 2 -odd particle (LOP) is stable and hence usually required to be neutral as there exist strong bounds on the presence of stable charged particles in the universe [1][2][3]. The typical collider signature of such a BSM scenario is missing transverse energy (MET) caused by the invisible LOP escaping the collider. However, there are scenarios where a heavy Z 2odd particle can become sufficiently long-lived to appear as stable in a collider experiment. This particle can be charged and thus produce a very distinct signature.
Heavy stable charged particles (or HSCP) can appear when the next-to-lightest Z 2odd particle is nearly mass degenerate to the LOP and hence its decay is kinematically suppressed. A prominent example for such a situation is a supersymmetric scenario where a wino-or higgsino-like neutralino is the lightest supersymmetric particle (LSP). The typically small mass splitting between the lightest neutralino and the lightest chargino in such a scenario can render the chargino long-lived, see e.g. [4]. A similar situation might occur in models of extra dimensions [5]. Other supersymmetric scenarios with mass degenerate sparticles have been proposed in order to address the discrepancy between the 6 Li, 7 Li abundances predicted in standard big bang nucleosynthesis (BBN) and those inferred from astrophysical observations [6]. In this scenario a bino-like neutralino and the lightest charged slepton are close in mass, resulting in a long-lived slepton [7]. Another possible scenario leading to HSCPs corresponds to the case where the LOP does not share the SM gauge couplings and only interacts super weakly. An important example are supersymmetric scenarios where the gravitino or axino is the LSP [8,9]. In this case the couplings are suppressed by powers of the Planck scale or the Peccei-Quinn scale and the life-time of the next-to-LSP (NLSP) can easily exceed the typical time for passing the detectors by many orders of magnitude.
Collider constraints on HSCPs have been mostly presented in specific BSM models [10][11][12] and in most cases cannot be directly applied to other BSM scenarios. Currently, Simplified Model Spectra (SMS) [13][14][15][16] have become a popular alternative for presenting less model-dependent constraints, which can then be systematically applied to specific BSM scenarios. The CMS and ATLAS collaborations typically interpret their results in terms of simplified models and methods to use these interpretations in a systematical way have been made publicly available for missing energy (MET) topologies by tools such as SMod-elS [17] and Fastlim [18]. As a result it is now possible (under some approximations) to test general BSM scenarios with MET signatures against LHC data. So far a similar approach has not been considered for HSCP or mixed MET-HSCP scenarios, although it has been argued [19,20] that SMS are particularly suitable for parametrizing the LHC sensitivity to HSCP signatures, since these are rather inclusive and depend almost only on the kinematics of the HSCP itself.
In this study we examine limits on the particle spectrum of BSM theories containing HSCPs making use of simplified models. For this purpose we introduce a set of eight simplified model topologies containing either one or two HSCPs in the final states and compute efficiencies for several values of the BSM masses appearing in each topology. The resulting efficiency grids or efficiency maps allows us to re-interpret previous results in a wide range of BSM models. In order to do this in a more general framework we incorporate these efficiency maps to the program package SModelS [17,21]. Since the public version of SModelS already contains a large number of LHC constraints on simplified models containing MET signatures, our modified version allows us to simultaneously apply both MET and HSCP constraints to full BSM models. Hence we are able to test scenarios where HSCP and MET constraints compete. As we will show, current LHC searches provide a high sensitivity to HSCPs [11,12] and the HSCP signature can support an exclusion or discovery even if its contribution to the total BSM signal is subdominant.
As a HSCP lead to a non-standard signature that is not supported by common fast detector simulations particular attention has to be drawn to a reliable implementation of such an analysis. In this study we use a novel approach for the computation of signal efficiencies presented in Ref. [11]. This method uses a parametrization of the CMS detector response as a function of the kinematic properties of the HSCP and allows us to accurately compute the CMS signal efficiencies for arbitrary models.
As an application of the SMS framework to HSCPs, we consider the stau co-annihilation strip in the Constrained Minimal Supersymmetric Standard Model (CMSSM). We focus on the nearly mass degenerate neutralino LSP and stau NLSP, which has been proposed in order to solve the 7 Li problem [22]. Using our extension of SModelS, we study the implications of both MET and HSCP searches and show that all of the parameter space (with tan β = 10) consistent with a potential solution to the 7 Li problem and the observed Dark Matter relic abundance is excluded by either HSCP or MET constraints (or both). This region of the CMSSM parameter space has also been studied in [23,24], although not focussing on the region of interest for the solution of the 7 Li problem. However, in [24] no particular attention was drawn to the derivation of the efficiencies for the HSCP signal and as an approximation the cross section limits from the inclusive stau production presented in [10] were used. By applying the SMS framework to the considered slice of the CMSSM our derivation of the HSCP constraints provide a significant improvement to the previous work.
The remainder of this paper is organized as follows. In Sec. 2 we will define the simplified models used in our results and present the computation of efficiency maps as well as their validation against the results presented in the CMS search [11]. In Sec. 3 we explain our implementation of the decomposition of a full BSM model into simplified model topologies and how these are used to constrain the full model. An application of our results to a CMSSM scenario with a nearly mass degenerate neutralino and stau will be presented in Sec. 4. We conclude in Sec. 5.

Simplified Models for HSCPs
In this section we will briefly review how CMS searched for HSCPs as it has a direct impact on some of the choices that we make in this paper. We will then introduce a number of simplified models that contain one or two HSCPs in the final state. These models correspond to the simplest topologies and appear in several BSM theories. The models are summarized in Tab. 1 (two HSCPs in the final state) and Tab. 2 (one HSCP and one neutral BSM particle in the final state). For computing the efficiencies as a function of the topologies and the BSM masses appearing in the cascade decays, we must choose a specific BSM model for production and decay of the Z 2 -odd particles. Here we use the supersymmetric (SUSY) simplified models, which are listed in the last column of Tabs. 1 and 2. In these, the HSCP is either the lightest chargino, χ ± 1 , or the lighter stau, τ 1 . As part of the Simplified Model approximations, we assume that the efficiencies are weakly dependent on the spin of the HSCP.

Overview of the CMS Search for HSCP
Heavy stable charged particles are highly penetrating particles that are expected to cross the entire CMS detector and reach the muon system 1 , and are therefore experimentally reconstructed and identified as muon particles. However, because of their large mass and the limited energy available in LHC collisions, they will be travelling through the detector with a velocity (β) significantly slower than the speed-of-light. Consequently, they will have an anomalously high ionization energy loss (dE/dx) and a longer time-of-flight (TOF) than relativistic standard model particles. In the CMS search for HSCP [10], the CMS silicon tracker is used to measure the particle dE/dx, while the CMS muon system is used to measure the particle's TOF. The events are mostly selected online by a muon name diagram parameters SUSY topology (p T > 45 GeV) trigger. However, the trigger becomes inefficient when the particle velocity is too low (β < 0.45) due to the too long delay (> 25 ns) for the particle to reach the muon system causing a mismatch between the muon system information and the inner tracker information.
While there is no real standard model background to this search, instrumental back- grounds due to the mis-measurement of either dE/dx or TOF is not negligible. To predict the amount of backgrounds in the signal region CMS exploits the fact that the dE/dx and TOF measurements are uncorrelated for backgrounds. The track dE/dx and momentum variables are used to reconstruct the particle mass and further discriminate the HSCP signal from mis-reconstruction background peaking at low values of the reconstructed mass. Although the mass threshold used in [10] is continuous, the required inputs for the reinterpretation of these results [11] are only provided in 100 GeV steps. Below we use the results presented in Ref. [11] to compute the signal efficiencies for the simplified models introduced in Tables 1 and 2.

Computation of signal efficiencies
In order to compute the efficiencies for the simplified models, we perform a Monte Carlo simulation of the signal at the 8 TeV LHC. For each topology listed in Tabs. 1 and 2 we scan over the respective BSM masses (listed in the third column) and generate 30 k events for each set of masses. For the event generation we use MadGraph 5 [25] to generate parton level events and then Pythia 6 [26] to perform the decays, as well as showering and hadronization. No detector simulation is performed, since we follow the fast simulation procedure defined in Ref. [11], where signal acceptances for HSCP candidates are provided as a function of the HSCP's kinematics. In order to identify HSCP candidates in each event we must first apply the following isolation criteria: where the first (second) sum includes all the charged (visible) particles in a cone of ∆R = ∆η 2 + ∆φ 2 < 0.3 around the direction of the long-lived particle, p T j denotes their transverse momenta, E j their energy and |p| is the magnitude of the long-lived particle's three-momentum. In both sums the long-lived particle candidate itself is not included. As muons release very little energy in the calorimeters they are not considered as visible particles. The purpose of these isolation requirements is to mimic the event selection used in the CMS analysis [10]. Long-lived particles failing any of these isolation requirements are not considered as HSCP candidates.
Once the HSCP candidates are identified, we can compute the signal efficiencies using the acceptances provided in Ref. [11]. These acceptances are given as probabilities for the candidate to pass the on-and off-line selection criteria (P on and P off ) and depend on the candidate's pseudo-rapidity η, transverse momentum p T and velocity β. The final signal efficiency (ǫ) is then given by: where P on (P off ) is the on-line (off-line) probability for each event, the sum runs over all generated events, N , and k i = (η i , p Ti , β i ) contains the kinematic properties for the HSCP candidate in the ith event. For events containing two HSCP candidates, the above probabilities must be replaced by [11] where k 1,2 i are the kinematical vectors of the HSCPs. There are two main effects governing P on/off . On the one hand, the velocity β should considerably deviate from 1 in order to allow for a discrimination against muons. Hence, for β → 1 the acceptance goes to zero. On the other hand for too small β (β 0.45) the particle may not be assigned to the right bunch crossing anymore. In this case, the trigger efficiencies (online selection) go down very drastically.
The CMS analysis also requires a minimum reconstructed mass (m rec ) for the candidate. For the fast simulation method used here, the collaboration provides the P on/off probabilities for four distinct mass cuts, which we consider as four different signal regions: Due to detector resolution effects, the reconstructed mass is typically m rec ≃ 0.6m HSCP [11] and the above requirements must be translated to the real HSCP mass. Therefore, when computing the efficiencies for each signal region, we take ǫ = 0, if m HSCP < 166 GeV, 334 GeV and 500 GeV for the signal regions SR 100 , SR 200 and SR 300 , respectively.

Validation
In order to validate the procedure described in Sec. 2.2, we compute the efficiencies for the gauge-mediated supersymmetry breaking (GMSB) models considered by the CMS collaboration in Ref. [11]. These models have a gravitino LSP and a long-lived stau as the NLSP. As a result, for collider purposes, all the sparticles cascade decay to the lightest stau, which is the HSCP candidate. We simulated the signal with Pythia 6 and analyzed the generated events as described in Sec. 2.2. The results obtained for the inclusive production of staus are shown in Fig. 1, where we also show the corresponding efficiencies obtained by the CMS collaboration. As in Ref. [11], we choose SR 0 for m HSCP < 166 GeV, SR 100 for 166 GeV < m HSCP < 334 GeV and SR 200 for higher masses. Our efficiencies agree within 3% with the ones obtained by CMS, where the differences are likely due to Monte Carlo statistical uncertainties. Therefore our procedure for computing the signal efficiencies reproduce very well the experimental results and can be used to produce the efficiency maps for the simplified models listed in Tabs. 1 and 2. In Fig. 1 we also reproduce the 95% CL limits on the inclusive production cross sections, which again agree very well (within ∼ 3%) with the ones obtained by CMS from Ref. [11]. Note that this limit is based on the discrete mass cuts on m rec mentioned above, the full CMS analysis allows for a event-based mass cut, resulting in somewhat stronger constraints for some values of the HSCP mass.   for the GMSB model as the function of the stau mass. We compare the CMS analysis (CMS-EXO-13-006 [11]) from the full detector simulation (red solid lines) with our implementation of the analysis described in Sec. 2.2 (blue dashed lines). In the lower frames we show the respective ratios ǫ CMS /ǫ Our , σ CMS limit /σ Our limit . Signal efficiency ǫ as a function of the mass of the HSCP for model M1 (direct production of two HSCPs, blue solid curve) and model M2 (direct production of one HSCP and one invisible particle, red dashed curve).

Results for the efficiency maps
Our procedure to compute the efficiency maps using the method outlined in Sec. 2.2 is as follows. For the eight models listed in Tabs. 1 and 2 we generate events and compute the efficiencies in a wide range of sparticle masses between 50 GeV and 3 TeV, varying the HSCP mass and all other masses listed in the tables in steps of 10 GeV and 50 GeV, respectively. In order to allow for a fast processing within SModelS we then reduce the number of mass points in regions of parameter space where the efficiencies do not vary considerably. As an example, we show in Figs. 2 and 3 the efficiencies for the SR 100 signal region (m rec > 100 GeV or m HSCP > 166 GeV) for the simplified models M1, M2 and M3, M8, respectively.
The signal efficiencies for the models with two HSCPs can be as high as 70% and strongly depend on the HSCP mass. For direct production of two HSCPs (simplified model M1), the signal efficiency stays above 20% for masses between 166 GeV and 1.4 TeV. As discussed in Sec. 2.2, we take ǫ = 0 for m HSCP < 166 GeV in the SR 100 signal region. This is the reason for the sharp drop in the efficiencies in the light HSCP region, seen in Fig. 2. For very large masses the particles are produced very close to threshold, with extremely small β ( 0.5) and the signal efficiency drops again. This is a result of the low detection probabilities (in particular the on-line probability) due to small trigger efficiencies for β < 0.5. In models M3 and M8 the HSCPs are produced in the decay of heavier particles and the efficiencies are largest in the intermediate mass range 500 GeV to 1.2 TeV. For large mass gaps between the produced particle and the HSCP, the latter becomes extremely boosted making the discrimination against the muon background extremely hard. Hence, the signal efficiency decreases rapidly in this region. This decrease is less pronounced in model M8, as the three-body decay leaves less energy to the HSCP.
As mentioned above, in the SMS framework used here we neglect sub-leading effects due to the spin of the HSCP or any other BSM particle. Nonetheless, we checked that, for direct production of the HSCP, the differences in the efficiencies are smaller than 20% when comparing the s-channel production of spin 0, 1/2 and 1 particles as well as the t-channel production of spin 1/2 particles. For the production of BSM particles through longer decay chains we expect the differences to be even smaller, as the kinematics of the HSCP are more strongly influenced by the mass spectrum of the model. We hence expect that neglecting spin effects as well as the effects of the production channels (s-channel versus t-channel) generates uncertainties of the order of 20% or below. Since these uncertainties are of the order of other theoretical uncertainties (such as NLO corrections to the sparticles production cross sections), we consider them acceptable.

Using Simplified Models to Constrain Full Models
The efficiency maps described in Sec. 2.4 allows us to compute the predicted signal cross section (σ th ) for a given simplified model M i in one of the four signal regions (SR j ): where σ M i is the cross section for the simplified model and ǫ M i SR j the respective efficiency for the signal region SR j . Comparing σ th with the experimental upper limit for the signal cross section in the respective signal region (σ UL ), it is possible to determine if the simplified model is excluded or not by the experimental searches. However, simplified models rarely match any model of interest and their usefulness relies on the fact that, under some approximations, it is possible to decompose a full model in terms of a coherent sum of simplified models (see Ref. [17] for details). In this case, the full model signal cross section (σ Full th ) to be confronted with σ UL is approximately given by: In the above expressionσ M i is the corresponding weight for the simplified model M i in the full model. These weights are computed by the decomposition procedure, which maps the full model into a sum of simplified models topologies. Since the decomposition method used here follows closely the one used by SModelS [17], we only outline the main steps, focusing on the differences required to treat long-lived particles. First, using as input an SLHA card, all the widths (Γ ), branching ratios (BRs), masses and production cross sections (σ) of the BSM states are defined for the input model. 2 Second, for each of the particles appearing in the production cross sections, we generate all possible cascade decays 3 , using the branching ratios and widths read from the SLHA card. However, since now the input model may contain quasi-stable states, we must determine if the particle at the end of the cascade decay is long-lived or not. More specifically, we must estimate what is the probability for a BSM state to have a prompt decay or to decay outside the detector. The fraction of particles which survive after traveling a distance l in the detector is given by: where Γ is the particle's width, γ = 1/ 1 − β 2 and β is the particle's boost. Therefore, the probability for a particle to decay promptly is: where l inner corresponds to the inner size of the detector, for which all decays are seen as prompt. For our subsequent results we take l inner = 10 mm. On the other hand, the probability for a particle to decay outside the detector is given by: where l outer corresponds to the detector size, which we take to be 10 m (for CMS). Clearly the above probabilities are event-dependent, since they depend on the boost of the unstable particle, through the factor γβ. Nonetheless we can still conservatively estimate these.
Since the efficiencies for a long-lived particle to be identified as a charged track fall sharply below β ≃ 0.45, here we take (γβ) outer = 0.6 (or β ≃ 0.5), which gives a mostly conservative estimate of F prompt . On the other hand, for prompt decays we take (γβ) inner = 10, which corresponds to β ≃ 0.995. Notice that for most models, Γ is such that the particle can be considered as decaying promptly or long-lived for a wide range of γβ values. Therefore, for most cases, our results are only mildly dependent on our choice of γβ.
Once F prompt and F long are known for each state, during the decomposition each unstable particle (with a non-zero width) will generate two possible topologies 4 : • one where the particle does not decay (it is considered as long-lived). In this case the topology weight will be proportional to the probability for the particle to decay outside the detector (F long ); 2 In order to read the SLHA file, SModelS uses the tools provided by the PySLHA [27] code. 3 Since the total number of all possible cascade decays is typically of the order of hundreds of thousands, we neglect all topologies which haveσ M i <σmin. In the results presented below we takeσmin = 5 × 10 −3 fb 5 × 10 −4 fb for the points with mHSCP ≤ 400 GeV (mHSCP > 400 GeV). 4 The case of displaced vertices is not considered in the present work.
... • one where the particle decays promptly. In this case the topology weight will be proportional to the probability for the particle to decay inside the inner detector (F prompt ).
Notice that in most cases we have (F long , F prompt ) ≃ (1, 0) for quasi-stable (or stable) particles or (F long , F prompt ) ≃ (0, 1) for unstable particles. Using the information in the SLHA input and a modified version of the package SMod-elS, which includes the computation of the F prompt/long probabilities, it is possible to decompose full models into distinct simplified model topologies and compute their weights. As mentioned above, the decomposition procedure starts with the pair production of BSM states (determined by the cross sections read from the SLHA card) and generates all possible cascade decays for each particle produced. For each step in the cascade decay, the topology weight is given by the product of the production cross section (σ prod ), the BRs for the decays appearing in the decay chain and the prompt decay/long-lived fractions, F prompt/long :σ where F X long F Y long is the product of the (non-decay) probabilities for the final states (X, Y ) appearing in the cascade decay chain. This procedure is outlined in Fig. 4. Since we only keep the topologies with a final weight above a minimum value (σ min ), the topologies containing very small probabilities (F prompt ≪ 1 or F long ≪ 1) are automatically discarded. This procedure has the advantage of allowing us to probe scenarios where more than one particle is (meta-)stable and to automatically determine which states can be considered as long-lived or decaying promptly. Furthermore, for models containing both neutral and charged (meta-)stable particles, the above procedure will produce topologies with both HSCP and missing energy (MET) signatures (or mixed MET-HSCP). Therefore it allows us to simultaneously confront the corresponding model with both MET and HSCP searches.
Once the topology weights (σ M i ) are known, through Eq. 3.2 it is possible to compute the full model signal cross section for each of the signal regions considered by the experimental searches. It is important to notice that, since we only computed efficiencies for the simplified models appearing in Tabs. 1 and 2, the signal from any other topolo-  gies appearing during the decomposition procedure are not included in the final signal.
Although this leads to conservative predictions, we will show below that, for the models studied here we only underestimate the signal cross section by 20% or less. At this point it is also relevant to stress that the current public version of SModelS does not include efficiency maps for the MET searches. For these the topology weights cannot be summed up and the experimental upper limits for the individualσ M i weights are used instead (for more details see Ref. [17]). Since the decomposition procedure can produce topologies with one invisible and one HSCP in the final states (such as the ones shown in Tab. 2), these can be constrained by both MET and HSCP searches. However, since HSCP constraints are typically stronger, the mixed topologies are considered as containing a single HSCP and we only apply the constraints from HSCP searches. Furthermore, since the HSCP efficiencies for these mixed topologies are almost independent of the cascade decay ending in the invisible state, we can neglect the kinematics of the MET branch and compress it to a single step decay, as shown in Fig. 5. This allows us to use the efficiencies computed for the simplified models M 4 and M 6 for a wide range of HSCP-MET topologies and improve the coverage of our efficiency maps. Using the efficiency maps for the simplified models M1-M8 computed in Sec. 2.4 and the decomposition procedure described above, we proceed to apply our modified version of SModelS to a physical model of interest. As we will show, both MET and HSCP searches can be relevant (although the former are typically stronger) and allows us to impose strong constraints on the stau co-annihilation region of the CMSSM.

Application to the Lithium-7 Problem
In this section we apply the procedure outlined in Sec. 3 to a full model containing long-lived particles. One interesting motivation for the existence of long-lived particles (in cosmological scales) is the Lithium-7 problem [28,29] (for a recent review, see Ref. [30] is highly inconsistent with the experimentally measured Lithium abundance [31]: Although some of the proposed solutions to the above discrepancy do not involve new physics (e.g. stellar depletion or inclusion of new nuclear reactions), these are usually highly tuned or require modification of nuclear rates well outside the expected uncertainties [32]. A popular alternative is to invoke new physics during BBN, which could explain a smaller Lithium production rate (or an annihilation of the original Lithium abundance). A well studied scenario is supersymmetry with long-lived staus ( τ ). If τ s are still present during BBN, they can form bound states with nuclei (such as 7 Li) and deplete the Lithium abundance, thus providing a viable solution to the Lithium problem. Since such solutions have been discussed at length in the literature [6,22,[33][34][35][36][37][38][39][40][41][42][43], here we concentrate on their features relevant for the LHC searches.
In this work we focus on the case of a neutralino LSP and consider the CMSSM, closely following the discussion presented in Ref. [22]. In order to cover the CMSSM parameter space we perform a Monte Carlo scan in the input parameters: where m 0 is the universal soft scalar mass, M 1/2 is the universal soft gaugino mass and A 0 the trilinear soft term, all defined at the unification scale, M GUT ≃ 2 × 10 16 GeV. We take the supersymmetric mass term µ to be positive (µ > 0), while we fix the ratio of the Higgs vacuum expectation values to be tan β = 10. We also limit our results to negative values of A 0 , since these enhance the radiative corrections to the Higgs mass. Similar results would be obtained for µ < 0 and A 0 > 0 except in this case we would have a larger discrepancy between the predicted and measured values for the anomalous magnetic moment of the muon (g − 2) µ . The supersymmetric spectrum is generated with SPheno 3.2.1 [44], the sparticle production cross sections are computed using Pythia 6 and NLLfast [45][46][47][48][49][50][51] and the neutralino relic abundance is computed with micrOMEGAs 3.0.24 [52].
As computed in Refs. [6,41], in order to solve the Lithium problem, the stau yield and life-time must satisfy: Y 0 τ 10 −13 and τ τ 1 − 100s. (4.4) The latter condition requires the stau to be nearly degenerate with the LSP, which we assume to be the lightest neutralino. In particular, the mass difference δm = m τ 1 − m χ 0 1 must be significantly smaller than the τ mass (δm < 0.1 GeV). In this quasi-degenerate scenario, the stau abundance before its decay is related to the neutralino relic abundance by [41] Y 0 τ ≃ where Y 0 τ is the stau yield after freeze-out (and before decay) and Y χ 0 is the final neutralino yield (after staus have decayed), which can be obtained from its final relic abundance: (4.6) The neutralino freeze-out temperature, T f , can be well approximated by T f ≃ m χ 0 1 /25 for the parameter space considered below.
Before discussing the LHC constraints from MET and HSCP searches, we first impose the following set of minimal constraints to the CMSSM: • a neutralino LSP; • 120 GeV < m h < 130 GeV. 5 Although we keep points with Ω χ 0 h 2 > 0.1289 (which violate the Planck's 3σ upper bound on the Dark Matter relic abundance [53]), we will explicitly identify in our results the region consistent with Planck.
Since the left-right mixing of staus is proportional to A τ − µ tan β, the stau mass not only depends on the scalar mass parameter m 0 but also strongly on A 0 and tan β. On the other hand, the neutralino mass is mainly dependent on the gaugino mass parameter M 1/2 . Therefore, the requirement m τ 1 ≃ m χ 0 1 introduces a correlation between M 1/2 and m 0 , A 0 , tan β. In particular, as m 0 increases (for a fixed M 1/2 value), A 0 must increase (in absolute value) in order to enhance the stau mixing and reduce its mass, which must satisfy m τ 1 ≃ m χ 0 1 ∝ M 1/2 . Therefore, whilst scanning over A 0 and M 1/2 with flat priors, we limit the scan over m 0 to a gaussian distribution around the value predicted by the linear relations found in [22]. This dramatically increases the efficiency of obtaining points fulfilling the δm < 0.1 GeV requirement. We also checked that points outside the 2σ-band of the gaussian distribution never satisfy the mass splitting condition.

Scan Results
For the results presented below we scan (with 14 k points) over the ranges: −42000 GeV < A 0 < −1000 GeV, 630 GeV < M 1/2 < 1100 GeV and 144 GeV < m 0 < 463 GeV with tan β = 10 and µ > 0. The requirements on the stau-neutralino mass splitting drastically restricts the allowed CMSSM parameter space. In Fig. 6 we show the m 0 -M 1/2 plane along with the corresponding values for the Higgs mass and Ω χ 0 h 2 . All the points shown satisfy δm < 0.1 GeV and contain a neutralino LSP. As we can see, at the right edge of the points shown, the Higgs mass falls below 120 GeV, while for M 1/2 1 TeV, the neutralino relic density violates Planck's upper bound. Furthermore, to the left of the points shown (low m 0 ), the stau becomes the LSP. We also show values for the stau relic abundance, computed according to Eqs. 4.5 and 4.6. As we can see all points allow for the minimum stau yield required to provide a solution to the 7 Li-problem. Finally, we show in the bottom-right plot values for the lightest stau mass after imposing all the minimal conditions listed above. As shown, the constraints on the Higgs mass and the relic density (as well as δm) imply 200 GeV < m τ 1 < 460 GeV.

LHC Constraints
After identifying the region of the CMSSM parameter space consistent with Planck (Ω χ 0 h 2 < 0.1283), the Higgs mass and the solution to the 7 Li-problem we proceed to discuss the constraints from MET and HSCP searches at the LHC. As described in Sec. 3, we have modified the public version of SModelS to include the CMS search for HSCPs. Within this framework we can simultaneously apply the MET and HSCP constraints to the CMSSM parameter space. While the MET constraints directly make use of the upper limits on the production cross sections (for a given simplified model) provided by ATLAS and CMS, the HSCP constraints use the efficiency maps for the CMS exotic search [11], as described in Sec. 2.2. Since the cascade decays of the SUSY particles in the scenario considered here may end up either on the lightest stau or on the lightest neutralino, we expect both the MET and the HSCP searches to be relevant for the parameter space shown in Fig. 6. Therefore we simultaneously apply both constraints and we consider a point in the parameter space excluded if, for at least one of the MET or the HSCP searches, the signal cross section for a given topology (or sum of topologies in the case of HSCP searches), σ th , is higher than its corresponding experimental upper limit (σ UL ).
In Fig. 7 we show the σ th /σ UL ratio 6 in the m 0 -M 1/2 and m 0 -m τ 1 planes, as well as the lines for the upper limit on Ω χ 0 h 2 (dashed gray) and σ th /σ UL = 1 (solid black). Values up to M 1/2 ≃ 1 TeV can be excluded, which corresponds to m τ 1 ≃ 450 GeV. We also notice that the excluded region (where σ th /σ UL > 1) extends to higher M 1/2 values when m 0 reaches its highest values (right edge). In this region the stop mass is suppressed due to large |A 0 | values and the pair production of stops is considerably enhanced, thus resulting in higher signal cross sections. The CMS paper [11] provides a conservative limit on the stau mass, m τ 1 > 260 GeV, which is based on exclusive stau pair production and is fairly model independent. Here, however, we can derive a much more stringent bound, since we are able to consider the inclusive production of staus from production of heavier sparticles and their decays. 7 We can also compare the bound derived here with the one obtained by CMS for the inclusive stau production in the GMSB scenario (see Ref. [11] for details). In the GMSB case the inclusive production (for a given m τ 1 ) is smaller than in the CMSSM scenario discussed above, due to the presence of heavier squarks and gluinos. Nonetheless, CMS quotes m τ 1 > 500 GeV, which is higher than the one found here for the CMSSM scenario. This is mainly due to the fact that, while 100% of the GMSB signal considered by CMS goes to HSCP final states, a considerable fraction of our (CMSSM) signal goes into final states with missing energy (neutralino final states), thus reducing the reach of HSCP searches. Furthermore, the signal efficiencies for the events containing one or two HSCPs are smaller for the CMSSM than for the GMSB scenario, as a result of the different spectra as well as the fact that most of the events in the CMSSM signal contain only one HSCP. Nonetheless, we are still able to exclude all the region of parameter space (for tan β = 10) consistent with Planck's upper bound on the Dark Matter relic abundance. Therefore we conclude that, for low values of tan β, the solution to the 7 Li within the CMSSM is no longer compatible with the LHC and Dark Matter constraints.
Finally, we comment on the feature appearing around M 1/2 ≃ 750 GeV or m τ 1 ≃ 320 GeV in Fig. 7. As discussed in Sec. 2.2, for the HSCP searches we consider four signal regions (SR 0 , SR 100 , SR 200 and SR 300 ). For m τ 1 < 334 GeV the efficiencies for SR 200 and SR 300 are taken as zero, so the parameter space is constrained by the upper limits for SR 100 . Once m τ 1 > 334 GeV (M 1/2 > 800 GeV), the efficiencies for SR 200 are no longer zero and this signal region becomes the most constraining one, as shown by the sharp transition σ th /σUL Figure 7. Values for the signal cross section (σ th ) over the experimental 95% CL upper limit (σ UL ) in the m 0 -M 1/2 (left) and m 0 -m τ1 (right) planes. Points with σ th /σ UL > 1 are excluded by either MET searches or HSCP searches at the LHC.
m0 [GeV] Figure 8. Points excluded at 95% CL by HSCP (left) and MET (right) searches in the m 0 -M 1/2 plane. The distinct signal regions for the HSCP search from CMS-EXO-13-006 [11] are shown as light orange (SR 100 ) and dark orange (SR 200 ). Signal regions SR 0 and SR 300 were also considered but are less constraining than SR 100 and SR 200 for this model. For the MET searches we show by distinct colors the constraints from CMS [54] (dark blue) and ATLAS [55] (light blue) analyses. seen in Fig. 7. This transition, however, does not affect our results, since all the points in this region are excluded. Since in the model considered here the signal cross section splits into a HSCP signal and a MET signal, we expect both the MET and the HSCP searches to have a smaller reach than in a scenario where the signature is pure MET or pure HSCP. In order to compare the reach of MET searches against the one of HSCP searches, we show in Fig. 8 the most constraining (the one with the highest σ th /σ UL ratio) HSCP analysis (left) and MET analysis (right). As we can see, the constraints from MET searches exclude points up to M 1/2 ≃ 650 GeV, which corresponds to mg ≃ 1500 GeV or mq ≃ 1350 GeV. The most constraining topologies in this case are simply squark pair production followed by direct decay to the neutralino LSP. On the other hand, the HSCP searches can exclude up to M 1/2 = 1 TeV or mg = 2250 GeV and mq = 2050 GeV. As expected, the HSCP constraints allows us to exclude a much larger region of the parameter space.

Comparison with Full Simulation
As already mentioned, the method outlined in Sec. 3 used to obtain the results above relies on a few approximations. First, several details of the full model are neglected, such as the spin of the intermediate particles, the production channel (t-channel or s-channel), kinematical effects due to off-shell decays and others. Second, when computing the signal cross sections (σ th ) for the full model we can only include the simplified models M1-M8 for which we have computed efficiencies. Finally other small effects such as the interpolation of the efficiency map grid and the effect of neglecting topologies with weights belowσ min can also affect the final result obtained through the simplified models approach. Therefore it is relevant to compare the results obtained in Sec. 4.1 with a full Monte Carlo simulation. In order to make this comparison we select O(100) representative points from the CMSSM scan discussed in Sec. 4.1 and compute their signal cross sections for each of the HSCP signal regions performing the full simulation via Pythia 6 followed by the analysis detailed in Sec. 2.2.
In Fig. 9 we show the ratio σ th /σ UL for the best signal region as a function of the stau mass obtained from the modified SModelS version (red points) and from the full simulation (blue points). The lower frame shows the ratio between σ th /σ UL for the full simulation and SModelS. As shown in the figure, the agreement is within ∼ 20% for all of the mass range shown. As already mentioned, since we do not compute efficiencies for σ th /σUL Ratio m τ 1 [GeV] Figure 9. Comparison between the results obtained with the full simulation (blue) and the modified SModelS version (red) for the HSCP search. In the upper frame we show the ratios between the signal cross section (σ th ) and the 95% CL upper limit (σ UL ) for the two methods, whilst the lower frame shows their ratio, i.e. (σ th /σ UL ) Full /(σ th /σ UL ) SModelS .
all possible simplified models, we expect the excluded region obtained by SModelS to be conservative. This is indeed what is seen in Fig. 9, where the SModelS value for σ th is always below the one obtained with the full simulation. In order to guide the eye we also show σ th /σ UL = 1 as a dashed line, hence all points above this line are excluded. Even though SModelS is conservative, both SModelS and the full simulation exclude stau masses up to ≃ 450 GeV, as already found in Sec. 4.1.

Conclusions
Heavy stable charged particles (HSCPs) provide a prominent signature at the LHC and are present in several well motivated BSM scenarios. Most of the current experimental searches for HSCPs present their results for specific BSM models, thus making it difficult to apply them to other scenarios of interest. Here we have re-interpreted these results in terms of simplified models, what allowed us to derive constraints to a wide range of arbitrary BSM models not considered previously. To this end we have presented a new method to systematically decompose full BSM models as a coherent sum of simplified models containing both stable and quasi-stable particles. To this end we have defined a set of eight simplified models containing one or two HSCPs in the final states and computed the corresponding signal efficiencies as a function of the model parameters. With the inclusion of both the new decomposition method and the efficiency maps to the program package SModelS we are able to apply both MET and HSCP constraints to arbitrary BSM models containing a Z 2 symmetry.
We showed that HSCP constraints on full BSM models can be reliably applied through the simplified model framework presented here. The constraints obtained by SModelS on the signal cross-sections for the scenario studied here agreed within ∼ 20% with the full Monte Carlo simulation. These differences are similar to other theoretical uncertainties and, when translated to constraints on the sparticles masses, do not lead to any significant difference between the full simulation and SModelS results. Therefore we conclude that the simplified models introduced here along with the SModelS tools are well suited for confronting full models with experimental HSCP searches.
We then applied our modified SModelS program to the CMSSM stau co-annihilation strip, particularly considering the case of a nearly mass degenerate stau and neutralino. In this part of the parameter space the stau becomes long-lived providing the HSCP signature and presenting a potential solution to the Lithium problem. As the decay chains following the production of heavier SUSY particles can terminate in either the neutralino or the stau, we encounter in this scenario both MET and HSCP signatures. We have shown that the MET constraints for single simplified model topologies allow us to exclude points up to m τ 1 = 275 GeV (or mg ≃ 1500 GeV and mq ≃ 1350 GeV). On the other hand, using the efficiency maps computed here, the HSCP searches can exclude up to m τ 1 = 450 GeV (or mg = 2250 GeV and mq = 2050 GeV). For small tan β values, the HSCP searches exclude the whole parameter space consistent with a potential solution to the Lithium problem and Planck's bound on the neutralino relic abundance.