Long-Lived stop at the LHC with or without R-parity

We consider scenarios of gravitino LSP and DM with stop NLSP both within R-parity conserving and R-parity violating supersymmetry (RPC and RPV SUSY, respectively). We discuss cosmological bounds from Big Bang Nucleosynthesis (BBN) and the gravitino abundance and then concentrate on the signals of long-lived stops at the LHC as displaced vertices or metastable particles. Finally we discuss how to distinguish R-parity conserving and R-parity breaking stop decays if they happen within the detector and suppress SM backgrounds.


Introduction
at least one charged lepton is produced in the decay.
The paper is organized as follows: we discuss the general setting and the stop-gravitino couplings in section 2 and then proceed to the cosmological constraints from BBN and the DM density in section 3. In section 4 we give our analysis of the LHC reach for the case of displaced vertices in the pixel or tracker detectors and of metastable tracks. Such analysis is independent on the decay channel, even if we will mostly assume that the decay product contain at least one lepton, e.g. a muon, that can be well-measured also away from the central part of the detector. We will give the reach of LHC depending on the stop mass and lifetime for both the different signals and see that there are interesting parameter regions where two types of signals can be observed. Finally in section 5 we will discuss how to distinguish the RPC and RPV scenarios, assuming the decay happens in the detector. We will then conclude in section 6.

Supersymmetric stop in the MSSM with and without R-parity
We consider a supersymmetric model of the MSSM type given by the most general renormalizable superpotential. The R-parity, baryon and lepton conserving part includes the Higgs µ term and Yukawa interactions as [50] where the indices i, j, k run over the three generations of fermions. The labels L i , E i , Q i , D i , U i , H u and H d are for the chiral superfields containing the lepton doublet, lepton singlet, quark doublet, down-type quark singlet, up-type quark singlet, up-type Higgs doublet and down-type Higgs doublet, respectively. We denote with y u ij , y d ij , y e ij the dimensionless Yukawa coupling parameters, while µ is the supersymmetric higgsino mass parameter.
The baryon violating part of the superpotential is given by while the lepton violating part is where λ ijk , λ ijk , λ ijk are R-parity violating Yukawa couplings and µ i is a parameter with dimension of mass mixing the Higgs with the lepton multiplet. If they all appear together, these parameters are strongly constrained by low-energy observable, i.e. proton decay, therefore in the following we will consider only the case of either baryon violating or lepton violating RPV. Apart for the superpotential, we introduce in the lagrangian also the soft supersymmetry breaking terms, in particular gaugino and scalar masses and soft trilinear terms A i . Then the squared-mass matrix for the top squarks in the gauge-eigenstate basis (t L ,t R ) is given by where is a non-diagonal hermitian matrix where A t , m t , and tan β denote, respectively, the trilinear coupling of the Higgs with top sfermions, the top quark mass, and the ratio of the two Higgs vacuum expectation values tan β = v u /v d . The masses m 2 t R and m 2 t L arise from the soft breaking,the D term contribution, and the top Yukawa coupling as follows: where θ W denotes the weak mixing angle and m Z is the Z 0 boson mass. The soft breaking masses mŨ 3 and mQ 3 are model-dependent. We see that in general the stop mass matrix can have a large off-diagonal entry, in particular if A t is chosen large to explain the Higgs mass [51][52][53]. In such a case the two mass eigenstate repel each other, so that the lightest one can become much lighter than the average mass scale. The stop mass matrix m 2 t can be diagonalized by a unitary matrix to give mass eigenstates: (2.7) We define heret 1 as the lightest state. Note that different SUSY breaking scenarios can account for a stop NLSP with a gravitino LSP. As an example we mention that in gauge-mediated supersymmetry breaking (GMSB) scenarios [54] the supersymmetry-breaking scale is typically much smaller than in the gravity-mediated case, so that the gravitino is almost always the LSP. Moreover in the recently proposed model-independent framework of general gauge mediation (GGM) [55,56] essentially any MSSM superpartner can be the NLSP. Here we will assume that the NLSP is the lightest stopt 1 . At the present time many studies of light stops decaying promptly have been performed, see for instance [44,45] in the context of gravitino LSP, since the stop is the superpartner mostly connected to the Higgs and required to be light to solve the hierarchy problem. In this paper, we consider only non-prompt decays instead. We take the stop masses as strongly split and we concentrate on the lightest stop mass eigenstatet 1 , assuming thatt 2 and the rest of the colored supersymmetric particles are outside the reach of LHC. Note in any case that most of our results are very weakly dependent on the stop mixing angle and therefore valid also if the second stop is not too heavy, as long as its production is suppressed. In case the other colored states like the gluino and first two generation squarks are within the reach of the LHC, we have additional particle production channels and the search becomes more promising.

Stop NLSP couplings, production and decay channels
In this section we consider the interactions of the stop NLSP in supersymmetric models with gravitino LSP and DM candidate.
The main interactions of the stop NLSP are the R-parity conserving QCD couplings, which in general dominate the stop pair production. In fact the R-parity violating couplings considered here are many orders of magnitude smaller than the QCD gauge coupling and too suppressed to give a measurable single-stop production. In the limiting case when the rest of the colored states are too heavy to be produced efficiently, the stop production cross-section is dominated by the direct production via the quark-antiquark annihilation and the gluon fusion channels. Then the stop mass is the only supersymmetric parameter influencing the production cross-section at tree-level and the dependence on the stop mixing arises only at NLO [57]. In this paper we simulate the stop pair production at the LHC with MADGRAPH 5 [58], which includes only the LO cross-section and therefore we neglect any mixing angle dependence of the production. Note that NLO corrections can change the cross-section by a of factor 50-70% within the mass range investigated here [57]. We will take such correction into account in our results by multiplying our cross-section by a constant NLO k-factor.
Let us now consider the decay channels of the stop and anti-stop pairs in the RPC and RPV models. In order to do that, let us recall that the gravitino comes into play when the supergravity is taken into account, i.e. supersymmetry is promoted to a local symmetry. In a sense, the gravitino plays the role of "gauge fermion" for local supersymmetry and its couplings are fixed by supergravity. Therefore the coupling of stop with the gravitino is described by the SUGRA R-Parity Conserving lagrangian where D νtR = (∂ ν + ig i A i ν )t R and P L (P R ) is the projection operator projecting onto left-handed (right-handed) spinor. Here, A i ν denotes any SM gauge boson whereas M P = (8πG N ) −1/2 is the reduced Planck mass. The interaction lagrangian of the left-handedt L has an analogous form. We can easily obtain from this expression the couplings for the lightest stopt 1 as In the rest frame of the decayingt i , taking the matrix element at the order 1/m 2 3/2 , the decay rate is independent of the mixing angle, which appears only in the interference at order m t /m 3/2 , and is given by the formula where m t = 173 GeV is the top mass, mt is the stop mass and m 3/2 is the gravitino mass. Now, if we neglect the top mass and the gravitino mass in the phase-space, and normalize the stop mass to the value 500 GeV and the gravitino mass to 1 GeV, we obtain a a stop lifetime given by We see that the stop lifetime can cover a very large range of values since it strongly depends on the stop and gravitino masses. For instance, choosing mt 1 = 800 GeV, a lifetime of τt 1 0.018 s or τt 1 1.9 × 10 4 s for a gravitino mass of m 3/2 = 0.1 GeV or m 3/2 = 100 GeV respectively, can be achieved. Finally, it is worth noting that this decay rate depends only on the particle masses.
Let us now turn to the RPV couplings, in particular in the case of bilinear RPV given in [9]. Other models with unstable gravitino DM are given in [59][60][61]. In that case we assume that the only RPV coupling is given by which can be rotated away from the superpotential if we redefine L i and H d as However, in this new basis, R-parity breaking is reintroduced in the form of trilinear R-parity violation. So much so that, one obtains the new trilinear R-parity violating terms and i are the Bilinear R-Parity breaking parameters. We see therefore that the lightest stop can decay via this λ -type coupling h through its LH component. The corresponding R-Parity violating lagrangian is in fact given by where the last expression is given in term of the stop mass eigenstate, the SM vacuum expectation value (vev) υ and the bottom mass.
In the rest frame of the decaying stopt 1 , if the antilepton masses are neglected, the total decay rate for equal i = reads then 17) and so, using the explicit value for vev υ and bottom mass m b = 4 GeV, the corresponding lifetime τt 1 is Note that this stop decay is also present in the case of trilinear lepton number RPV, while it is absent in the case of baryon number violating RPV.
We consider here this RPV stop decay as particularly promising because it contains leptons in the final state, that can be more easily detected at collider experiments also away from the central collision region and are therefore a very favorable signal.
In the case of baryonic violating RPV or even in MFV models like [62,63], the relevant superpotential coupling is given by the λ coupling and the lagrangian reads instead (2.19) with λ 3jk antisymmetric on the last two indices, giving the decay rate into two light-quark jets as For small λ this decay can also lead to displaced vertices, with two jets originating far away from the stop pair production vertex.
In any of the scenarios we discussed here, we see that the stop lifetime is always much longer than the hadronization time and therefore we expect the stop and anti-stop to hadronize into an R-hadron before decay [64]. Such an R-hadron can in principle be both electromagnetically charged or not and even change its charge while it travels in the detector [65].
For the sake of clarity and transparency of exposition, from now on we adopt the convention that the lightest stopt 1 is simply calledt.

Stop NLSP and gravitino LSP in cosmology
In this section we discuss very shortly the cosmological bounds on the scenario with stop NLSP and the gravitino DM and LSP, in order to single out the cosmologically preferred parameter space in both RPC and RPV models.
Let us consider first the effect of a stop NLSP during BBN. The stop is a colored and EMcharged particle and therefore it can disrupt BBN not only through the energy release in the decay, but also because of bound state effects [66]. In the first case the light element abundances are more strongly affected by hadro-dissociation and therefore the limits are more stringent for hadronically decaying particles like the stop [67,68]. In the latter case the constraints are independent of the decay channel and just depend on the stop lifetime and density at the time of decay [66,69]. They therefore apply equally to any of the scenarios we discussed.
The limits on the abundance Y χ (τ χ ) = n χ /s from bound state effects for a hypothetical longlived strongly interacting massive particle χ, have been computed by M. Kusakabe et al. in [69]. Requiring that the primordial light element abundances remain within the observed ranges, they obtained the following constraints depending on the particle lifetime τ χ : In the window (30 s < τ χ < 200 − 300 s) the most stringent constraint comes from the upper limit on the 7 Li, while for (200 − 300 s < τ χ < 2 × 10 3 s) the strongest constraint is due to the upper limit on the B abundance. Finally the upper limit on the 9 Be abundance determines the bound for longer lifetimes. These constraints are very steep and become quickly dominant over the hadro-dissociation bounds [67,68]. For lifetimes shorter that approximately 30 s the constraints from bound states disappear and those on hadronic decays are very weak [70] and this corresponds in our two scenarios to a stop mass range We see therefore that BBN does not provide practically any bound on the RPV scenario, apart if the RPV coupling is very small, below 10 −12 . We will therefore in the following only consider in detail the constraints for the RPC case.

RPC decay of stop NLSP in cosmology
Assuming that the stop NLSP is in equilibrium in the thermal plasma, the density of the stop at chemical decoupling is given by the solution of a Boltzmann equation as in the case of a WIMP: where x = mt/T , H is the Hubble parameter, Y eq is the equilibrium stop abundance and σ ann v x is the thermally averaged annihilation cross section. The density of a colored relic and the BBN bounds have been studied in a model independent way by C. Berger et al. in [42] . We follow here their analysis just updating the constraints to those given above. In [42] the authors first considered the simplified case of a single annihilation channeltt * → gg. Such a choice was motivated by the fact that this channel just depends on the stop mass and its QCD representation, without dependence on the rest of the supersymmetric spectrum, and, in addition, it is always the dominant channel, contributing at least 50% of the total annihilation cross-section. This channel gives the most conservative result since it cannot be suppressed by particular choices of spectrum or mixing and it provides a conservative upper limit on the stop abundance. In fact, if more annihilation channels are taken into account, the cross section increases and, therefore, the stop density decreases. In [42] also the computation of the Sommerfeld enhancement [71] for this particular channel was performed with two different prescriptions for the higher orders. The Sommerfeld enhancement increases the cross-section at low velocity and can be obtained by resumming over the exchange of a ladder of gauge bosons between the initial particles. It was found that the averaged Sommerfeld factor reduces the treelevel yield by roughly a factor of 2, while the summed Sommerfeld one by roughly a factor of 3. We will therefore take the stop abundance from the leading order computation and vary it by a factor 2-3 to see the effect of both the Sommerfeld enhancement and the additional annihilation channels.
The stop abundance is proportional to the stop mass and it can be rescaled as up to logarithmic corrections, since, in general, the mass always enters linearly in the coefficient of Eq. In order to set limits on the RPC model, we first compute the stop density as a function of the stop mass from Eq. (3.3) and we compare it with the limits in [69]. We determine then the maximal allowed value of the stop lifetime for the particular mass by the value of τt 1 such that Through the analytical formula for the stop lifetime, we can then recast the bounds in the plane mt vs m 3/2 . We give these results in Fig. 1. We see that the constraints for the LO and the Sommerfeld enhanced case overlap and are almost in perfect agreement. This is due to the fact that the bound-state BBN constraints are very steep and do not change appreciably even if the stop density changes by a factor of a few. So the BBN bounds practically do not depend on the details of the stop freeze-out, as long as no strong resonant annihilation is present, and are quite robust and independent of the masses of the heavier superpartner and therefore of the particular supersymmetric model with stop NLSP 2 .

CDM constraints
We are assuming in this paper that gravitinos are Cold Dark Matter and therefore they must have obtained the required abundance in the course of the cosmological evolution. The gravitino production by the decay of the stop NLSP [72,73], is in most of the parameter space negligible since either the stop abundance or the gravitino mass are too small. Moreover such contribution is substantial only in the case of the RPC model and it is instead very much suppressed if the stop RPV decay is dominant.
On the other hand, gravitinos can be produced in substantial numbers by scatterings and decays of supersymmetric particles in equilibrium in the hot plasma. Their abundance is then proportional to the bath reheating temperature T R and can exceed the critical density of the universe if no restrictions on the reheating temperature is imposed [74][75][76]. In this context the gravitino abundance Ω 3/2 is given by [28,75,76] where M i are the physical gaugino masses and the coefficients γ i account for RGE effects between the reheating temperature scale and the scale of the physical gaugino masses. We have for those constants the ranges γ 1 = 0.17−0.22, γ 2 = 0.54−0.57, γ 3 = 0.48−0.52 from the 1-loop RGE for the gaugino masses and gauge couplings from T R = 10 7 − 10 9 GeV [28]. We neglect here possible decay of the heavier superpartners in equilibrium, the FIMP contribution [77], which may even play a dominant role in the case of hierarchical spectra between gauginos and scalar superpartners [78].
Note here the dependence on other supersymmetric masses than the LSP mass, such that the exact abundance becomes a model dependent quantity. Since we are here interested mostly on a constraint on the gravitino and stop mass, we will here consider as most conservative the case when the gaugino masses are not much heavier that the stop, M i /mt = (1.1 − 2), in order to minimize the gravitino production. The current best-fit values for Dark Matter density in the universe from the data of the Planck satellite is given by [79] Ω CDM h 2 = 0.1199 (27), (3.5) and results slightly larger that the previously obtained combination of the seven-years WMAP data, observations of baryon acoustic oscillations and determinations of the present Hubble parameter in [80] , We impose here that the gravitino energy density in Eq.(3.4) is smaller or equal to the the Cold Dark Matter density. From the equality to two measurements Ω CDM = 0.1126 − 0.1199 we obtain the yellow and brown curves in Fig. 2. On such lines the gravitino CDM density is fully produced by thermal scatterings, while below the line the gravitinos are overabundant and therefore the parameter space is excluded by the CDM constraint. Note that the position of the line depends on the particular T R assumed and the curves move up and down in the value of the gravitino mass exactly by the change in T R , since indeed the dependence is on T R /m 3/2 . Finally, comparing the two plots for different values of M i /mt, we can see that for a bigger ratio of physical gaugino masses to stop mass, the white allowed region is reduced.

All constraints
We plot now both the cosmological constraints together in the plane mt vs m 3/2 in Fig. 2. For future convenience we also show the region of the parameter space where the RPC stop lifetime is smaller than 10 −7 and therefore the stop decays inside the detector. Looking at this figure, we see that the allowed region is limited from above by the BBN constraint and from below by the CDM constraint, so that only a narrow allowed strip remains for T R = 10 7 TeV. The breadth of such strip depends both on the supersymmetric spectra and in particular on M i /mt and on the reheat temperature assumed. In particular for reheat temperatures above a few 10 7 GeV no allowed parameter space remains for stop masses below 2 TeV.
We note that the region where the RPC decay is sufficiently fast to happen in the detector correspond to a very low reheat temperature of the order of 10 3 − 10 4 GeV, so that in case of high reheating temperature the stops appear as metastable particles at the LHC.
In order to conclude the discussion of the cosmology of a gravitino CDM with stop NLSP, it is worth drawing the reheating temperature corresponding to the right value of Dark Matter density for m 3/2 = 1 GeV and the maximal allowed reheat temperature T max R as a function of the stop mass in n Fig. 3. Such curves are shown for two different values of the ratio of physical gaugino masses and stop mass in Fig. 3, e.g. M i /mt = 1.1 in green (solid line) and M i /mt = 2 in blue (dot-dashed line). We see that for larger NLSP mass, a smaller reheating temperature is needed to match the observed DM density at fixed gravitino mass, while on the other hand the bound on T R becomes relaxed for larger stop masses as the BBN bounds are weaker.

Decay of stop NLSP at LHC
Many extensions of the SM include heavy, long-lived, charged particles (HSCPs). These particles can travel distances comparable to the size of modern detectors, where they might be produced. Thus, they might appear to be stable if they have a lifetime bigger that a few nanoseconds. Moreover, the HSCPs can be singly charged (|Q| = 1e), fractionally charged (|Q| < 1e), or multiply charged (|Q| > 1e). Since the particle identification algorithms at hadron collider experiments generally assume signatures appropriate for SM particles, e.g., v ≈ c and Q = 0 or ±1e, nowdays the HSCPs might be misidentified or even completely missed without dedicated searches. The LHC experiments have already performed specific analysis, especially for the case of metastable particles [16][17][18][19][20][21][22]. In this very exciting background the goal of this section is to study two different classes of signal coming from a long-lived stop NLSP, which is produced by the proton-proton collision at LHC, at a centre of mass energy of √ s = 14 TeV and an integrated luminosities of L = 25 fb −1 and 3000 fb −1 . The first signal is represented by a displaced vertex inside the detector due to the stop decay inside the Pixel or Tracker detector. We will here mostly consider the kinematics and geometry of the CMS detector to estimate the number of decaying events within two adjacent detector parts. We neglect the interactions of the R-hadron with the detector material that could cause the particle to stop in the detector before the decay and the presence of a magnetic field bending the trajectory for charged R-hadron. Such effects could be taken into account only by a full detector simulation, which is beyond the scope of our study. Note in any case that interactions in the detector can only lead to a larger number of events inside the detector and that the magnetic field does not affect neutral R-mesons, which are expected to be around 50% of the cases [64].
The second type of signal is instead the HSCP track of a metastable stop that leaves the detector before decaying. Such a signal is actively searched for by the LHC collaborations.
In both cases we will use MadGraph 5 to compute the LO stop production at the LHC and we will correct it with a constant NLO k-factor of 1.6 corresponding to the k-factor given by Prospino [81] for a stop mass of 800 GeV. We checked that this factor remains in the range 1.5 − 1.7 for stop masses up to 2 TeV.
Regarding the decay, we will either include it within the MadGraph 5 analysis with a reference decay rate and then rescale the distances to probe the whole accessible lifetime range or use an analytical estimate for the distribution of the decay lengths. As we will see both approaches give similar results, with the semi-analytical one allowing to explore more easily the parameter space. Before going into detail, we want to highlight here that this analysis is independent from the stop decay channel, as long as the decay gives measurable tracks and a clear displaced vertex. In fact, we will apply our results to the parameter space of both RPC and RPV models that we have already discussed in Sec. 2. We are returning to this crucial point later in Sec. 4.6.
To this day, both the CMS and ATLAS experiments have published searches on this topic for HSCPs produced in proton-proton collisions. Their latest results can be found, respectively, in the papers [21,22]. In particular, the CMS analysis investigates signatures in three different parts of the detector using data recorded at a centre of mass energy of √ s = 7 TeV and 8 TeV and an integrated luminosity of L = 5 fb −1 and 18.8 fb −1 . The studied parts are the inner tracker only, the inner tracker and muon detector, and the muon detector only. On the other hand, the ATLAS collaboration is also investigating signatures producing a displaced multi-track vertex with at least a high transverse momentum muon at a distance between millimeters and tens of centimeters from the proton-proton interaction point, by using a data sample at √ s = 8 TeV and for L = 20.3 fb −1 [23].
As a detector, in this paper we will consider as explicit example the CMS experiment even if we expect ATLAS to have comparable reach, perhaps even larger due to the bigger size. So first of all we start describing very briefly how the CMS detector is made and what it searches for.

The CMS detector
The CMS detector uses a right-handed coordinate system where the origin is at the nominal interaction point. The x-axis points towards the centre of the LHC ring, the y-axis points up with respect to the plane of the LHC ring and, at last, the z-axis along the counterclockwise beam direction. The polar angle θ is measured from the positive z-axis, the azimuthal angle φ in the x − y plane and the radial coordinate in this plane is denoted by r. The transverse quantities, such as the transverse momentum ( p T ), always refer to the components in the x−y plane. In this context, the magnitude of the three-vector p T is indicated by p T and the transverse energy E T is defined as E sin θ.
In order to make easier the description of the CMS detector, the layout of one quarter of it was sketched in Fig. 4. Now, if we start from the innermost part of the detector and going outwards,   [25] • Interaction Point (IP) is the point in the centre of the detector at which proton-proton collisions occur between the two counter-rotating beams. We will assume that the stop and anti-stop pair is produced at this point.
• Pixel (Pi) detector contains 65 million pixels, allowing it to track the paths of particles emerging from the collision with extreme accuracy. It is also the closest detector to the beam pipe and, therefore, is vital in reconstructing the tracks of very short-lived particles. We therefore expect that the decay would be very well measured if it happens in this part of the detector.
• Tracker (Tr) can reconstruct the paths of high-energy muons, electrons and hadrons, as well as see tracks coming from the decay of very short-lived particles. It is also the second inner most layer and, so, receives (along with the Pixel) the highest number of particles. Even if it is less densely equipped than the Pixel detector, it can still recognize tracks coming from a displaced vertex instead than the interaction point.
• Electromagnetic Calorimeter (EC) is designed to measure the energies of electrons and photons with high accuracy via electromagnetic calorimeters. In our case it can allow to measure the energy of the lepton arising in the decay.
• Hadron Calorimeter (HC) measures the energy of hadrons and can give an estimate of the b-jet energy in the decay.
• Magnet (M) is the central device around which the experiment is built. The job of this big magnet ( B = 4 T), which contains all the parts above, is to bend the paths of particles and allow for an accurate measurement of the momentum of even high-energy particles.
• Muon System (MS) is able to detect muons and possibly other charged particles able to cross the whole detector.
The best detector parts to single out the presence of a displaced vertex are the pixel and tracker detectors and therefore we will restrict our discussion to the case of stop/anti-stop decaying there or surviving through the whole detector.

Numerical analysis
The numerical analysis in this paper is realized by means of the open source software MadGraph 5 which can generate matrix elements at tree-level, given a lagrangian based model, for the simulation of parton-level events for decay and collision processes at high energy colliders [58]. In order to study the stop production at the LHC experiment we choose the MSSM model from the Mad-Graph 5 library of models which is built upon that of the package FeynRules 3 [83]. Furthermore, we choose that the proton-proton colliding beams are not polarized and the centre of mass energy is √ s = 14 TeV for all our simulations. We run MadGraph 5 for several stop masses but for only one reference value of the stop decay rate that we take to be Γ cm t = 2.02159 × 10 −10 GeV, requiring that it has to generate 10, 000 events per run. By doing so, the kinematics of all particles in the process can be obtained. From those quantity we compute numerically the stop and anti-stop decay lengths t in each single event assuming that the particles propagate undisturbed from the interaction point and decay randomly according to an exponential distribution. From the decay length and the stop or anti-stop momentum direction we can draw using the software Mathematica the distribution of all decay vertices in the detector plane (r, z). If we note that any change in the stop decay rate can be compensated by a change in the distance that the stop travels, we can circumvent the problem of launching MadGraph 5 for all of decay rates we need for our analysis, by simply rescaling the dimensions of all parts of the detector consistently. This point will be clarified in the next subsection when we are discussing the semi-analytic analysis. By using this expedient, the stop length distribution for a stop mass of mt = 800 GeV and a stop decay rate of Γt = 2.0159×10 −16 GeV can be easily found from the reference data by rescaling the size of the detector by a factor of 10 −6 . Such a distribution is shown in Fig. 5 by red dots along with the size of pixel, tracker and the whole CMS detector. Now equipped with the stop decay length distributions, we can count how many stops or anti-stops decay within the Pixel or Tracker in the CMS detector at the integrated luminosity L = 25 fb −1 and the maximum expected one L = 3000 fb −1 . This allows us to obtain an estimate of the LHC reach in the plane of stop mass versus lifetime with just two very simple working hypotheses. First we neglect both the backgrounds of the SM and other SUSY particles than the lightest stop, which are respectively expected to have much shorter decay lengths and assumed to be too heavy to be produced at LHC. Secondly we set the detector efficiencies to 100% and declare 10 decays inside one CMS detector part sufficient for the discovery of a displaced vertex and 10 decays outside for the discovery of a metastable stop. We require 10 decays instead of just 2 or 3 in order to reduce the numerical fluctuations and obtain a more stable numerical result. The analytical estimates in the next section will allow to draw conclusions also for a different number of decays. By doing so, the LHC reaches in Fig. 6 are obtained.
Particularly, the Fig. 6 shows the MadGraph data-points in the stop mass-stop lifetime plane corresponding to 10 decays inside Pixel (blue dots), Tracker (red dots) and outside the detector (green dots) at both integrated luminosities, with L = 25 fb −1 shown in the left panel while L = 3000 fb −1 in the right panel. Here, as in the next plots, the data-points for the stop decays outside the detector are always labeled by "Out". Comparing the Pixel and Tracker reach, we are not surprised to see that they are pretty close to each other since their sizes are similar. The tracker detector is larger than the pixel one and this gives a slightly better reach, on the other hand we expect that the pixel detector would offer a better precision and efficiency in the measurement of a displaced vertex, so that this may overcome the geometrical advantage in a full detector simulation. It is interesting that the stop kinematical distribution is such that a similar number of events is often obtained in the two detector parts, allowing on one side for a cross-check and on the other to better disentangle a long-lived stop from any SM background like B-mesons decaying mostly in the pixel detector.
We see clearly from the plot that the searches for displaced vertices and escaping particles are complementary: the first covers the lower plane corresponding to short lifetimes, while the latter is mostly sensitive to the long lifetimes. Combining the two searches it is possible to cover the parameter space for any lifetime up to a maximal mass where the production cross-section starts to become too small to produce a sufficient number of stops. Assuming that the LHC does not observe any signal, we could therefore obtain a lifetime-independent lower limit on the stop mass at around mt 1300 GeV for the integrated luminosity L = 25 fb −1 and at around mt 2100 GeV for the integrated luminosity L = 3000 fb −1 .

Semi-analytic approximate analysis
The semi-analytic analysis is realized via analytical estimates for the decay length of the long-lived stop particles. Even in this case we use MadGraph 5 to compute the production cross section σ and, therefore, the number of generated stop particles N 0 at LHC, from the product of the

In both panels, the Pixel reach is denoted by blue points and its interpolation is dashed while the Tracker reach is denoted by red points and its interpolation is dotdashed. At last, the reach of metastable particles is denoted by green points and its interpolation is a solid line.
cross-section times the integrated luminosity (σL). We complement the previous analysis with this semi-analytic approach in order to have a better control of the physical parameter space, faster results and, at the same time, a useful check of the results of MadGraph 5. The semi-analytic analysis is based on the well-known exponential decay formula for a particle travelling in a straight line, giving the probability P (d) that a particle decays at d. It reads where Γ is the decay rate in the centre of mass frame, c the speed of light and βγ is the relativistic βγ factor, which is defined in term of the energy E and the three-momentum p of the decaying particle as βγ = (| p| /E)/ 1 − (| p| /E) 2 . The factor in front of the exponential in Eq. (4.1) is determined by the proper normalization of the probability P (d), i.e. the condition +∞ 0 P (d) dd = 1.
By using this exponential decay formula, the corresponding formula for the probability as function of a dimensionless coordinate y can be very easily obtained. It reads where y = dΓ/c and the normalization was obtained as for P (d). Here, we can explicitly see that any change in the decay rate Γ can be always compensated by an appropriate change in the distance d that the particle travels. Therefore, from the analytical expression of P (d) we can directly justify the rescaling procedure used for the MadGraph events, that allowed us to cover the whole range of decay rates from a single run. Both the formula for P (d) and P (y) can be generalized to describe an exponential decay of a sample of particles, integrating over the particles distribution in momentum and therefore βγ. To have a simpler expression, we will instead approximate such integral by assuming a single effective value of βγ and by multiplying simply by the initial number of particles N 0 . In this way, an estimate for the number of decaying particles of the sample as a function of d and y, called N (d) and N (y) respectively, are achieved. They are where the coordinate d = √ r 2 + z 2 can be taken as a function of the coordinates in the CMS detector. We neglect here in first approximation the bending of the trajectory by the magnetic field, which affects only the case of stop hadronization into a charged hadron, or the energy loss due to the interaction with the detector material, possibly negligible for a stop decay within the inner part of the detector.
To obtain the best approximation in Eqs. (4.3), (4.4), we compute the stop βγ distribution with MadGraph and compare the analytical formula above with different effective βγ with the decay length's distribution obtained directly from the MadGraph run. The βγ distribution is given in Fig. 7 for a stop mass of 800 GeV. We see that even for relatively small mass, the stop and anti-stop are mostly produced as non-relativistic, with a peak in the βγ distribution clearly below 1. We consider the analytical vertex distance distributions with different effective βγ, i.e. taking the value at the maximum βγ max = 0.66 or (1/βγ) max = 0.8026 or the average values βγ = 0.9207 and 1/βγ = 1.24595 . For the reference mass mt = 800 GeV we compare these distributions in the detector range 0 r(m) 17 with the numerically computed distribution of displaced vertices and we select the value of βγ = βγ max = 0.66 as the one giving the best fit. We Let us work now out more details from the analytical formula of N (d). Assuming just a spherical geometry for the detector parts, we can integrate the expression with respect to d from r i to r f , which respectively stand for the initial and final radial distance from the IP of the part of interest of the detector, and from r f to +∞ and so we have an analytical approximate expressions for both number of stop decays that occur inside the detector parts and outside the detector. They are given by the equations Since the number of particles generated by proton-proton collisions is given by the product of cross-section times luminosity, N 0 = σL, and a power-law formula for σ(mt) can be obtained fitting the MadGraph data, we can solve Eq. (4.5) and Eq. These results are plotted in Fig. 8, where the dashed blue line, the dot-dashed red line and the green solid line denote, respectively, the Pixel, Tracker and Outside reach. Observing Fig. 8, we can note that the crossing analytical curves identify in the plane four regions, which we have labelled A, B, C and D. In the three regions A, B and D at least one signal is seen for any stop lifetime. In particular in the region D both types of signal are accessible, allowing to cross-check the measurement of the stop lifetime. Only in region C no clear signal would be measured at the LHC. Of course such a region may be reduced by combining the two types of signals to have 10 events in total or by loosening the requirement to fewer events. To show the dependence on the requested number of vertices or metastable tracks, we give in Fig. 9 the LHC reach for N = {1, 10, 100} again for both luminosities L = 25 fb −1 (left) and L = 3000 fb −1 (right). It is clear here that the change in N affects strongly the reach of the metastable track search only at large masses, since for long lifetimes those events practically coincide to the total number of produced stops and anti-stops, which decreases very fast as a function of the stop mass. For displaced vertices, the change in N just shifts the curves to a larger/smaller stop lifetime.
Instead in Fig. 10 we give additional plots of the approximate reach for different LHC luminosities: in the left panel we plot the LHC reach for L = 100 fb −1 whereas in the right one the LHC reach for L = 300 fb −1 . We see from the four different plots in Figs. 8 and 10 that the LHC has the chance to cover the whole parameter space in lifetime up to stop masses of order 1300, 1500, 1700, 2100 GeV for a luminosity of 25, 100, 300, 3000 fb −1 respectively.

Comparison and discussion
To compare directly the numerical and approximate LHC reach in the stop lifetime-stop mass plane, we plot all the curves together first without NLO correction in Fig. 11. Looking at this figure it is blatant to see a good agreement between the MadGraph data and the approximate curves at both integrated luminosities: L = 25 fb −1 , showed in the left panel, and L = 3000 fb −1 , showed in the right panel. The analytical curves can be easily extended to consider ±1σ statistical error bars in the Poisson distribution, corresponding to 8 and 12 events respectively. We see that these curves give a very good description of the numerical data-points as they are mostly included in the ±1σ band. Only at low masses, where the statistics is large and the band becomes narrow, the analytical estimate deviates slightly from the numerical results. This is surely partially due to the spherical shape of the detector assumed in the approximated expression, which causes a substantial error in the region with a large number of events. In fact the data-points computed in the numerical analysis with the exact detector shape contain more signal events and give a wider reach in the low mass region. The current CMS excluded region for metastable particles (MP) [21], obtained at a centre of mass energy of √ s = 8 TeV and an integrated luminosity L = 18.8 fb −1 , given by our analytical curve for zero decays outside the detector, is also shown as a yellow region in the upper left corner of each panel. We see that for long lifetimes the curve excludes stop masses below 800 GeV and this coincide with the results obtained in [21]. Such bounds weakens for lifetimes below 10 −7 s. Here it is worth highlighting that the current excluded region still leaves a considerable allowed region to be investigated at LHC in the future.
Finally we investigate as well the impact of NLO corrections to the stop production crosssections. It is well known that processes involving colored particles are strongly affected by QCD NLO corrections and that those tend to increase the cross-section. It is therefore clear that the  LHC reach obtained with the LO cross-section is very conservative. We redraw the same plots by multiplying the production cross-section with the NLO k-factor of 1.6 and we observe that such a correction increases the yield substantially and extends the reach to larger stop masses. Indeed the NLO LHC reach in stop mass moves up by approximately 200 GeV, i.e. from 1200 to 1400 (2000 to 2200) GeV for 25 (3000) f b −1 integrated luminosity.

RPC and RPV models
As we have already mentioned beforehand, the discussion so far has been centered on the presence of a displaced vertex or metastable tracks and is so independent from the particular stop decay channel, as long as it contains sufficiently many charged tracks to make the secondary vertex visible. Therefore, our analysis is valid for both the RPC and the RPV stop decays we have discussed in Section 2.1. We can therefore interpret our results in terms of the two different model parameters.

LHC reach for the RPC stop decay
For the RPC stop decay into gravitino and top, the LHC reach can be easily reformulated in the plane mt vs m 3/2 by using the analytical formula for the RPC stop lifetime, Eq. (2.11). In Fig. 13 we display these curves at NLO for Pixel, Tracker and the part outside the detector by means of We note that the cosmological allowed region (white region) appears above all of three detector curves, which means that stop NLSPs with a consistent cosmology and high reheat temperature can be detected only as metastable particles. The detection of a diplaced vertex would instead point to the case of small gravitino mass and low reheating temperature.

LHC reach for the RPV stop decay
For the bilinear RPV model, the LHC reach can be reformulated in the R-Parity breaking parameter-stop mass plane vs mt by using the analytic formula of the Bilinear RPV stop lifetime given by Eq. (2.18).
In Fig. 14  Here we also display the corresponding MadGraph data by points which are blue for Pixel, red for Tracker and green for the decays that occur outside CMS. At the lower left-hand corner of each plots the current excluded region for metastable particle (MP) is painted yellow. At the top of each panel, on the other hand, we see the indirect detection excluded region for the gravitino dark matter decay. The gravitino decay, indeed, leads to a diffuse γ-ray flux which can be compared to the γ-ray flux observed by Fermi-LAT collaboration [10,11] so as to get a severe lower bound on the gravitino lifetime, and therefore a upper bound on the R-Parity breaking parameter . Note that the latter bound depends on the gravitino mass, but not the stop mass. Particularly, we take here the value of the upper bound on the R-Parity breaking parameter to be 2 × 10 −8 for a gravitino mass of the few GeVs from [13]. We see in this case that the analysis of displaced vertices is absolutely needed to close the gap between the indirect detection bound and the possible metastable particle constraint. Indeed in this case both displaced vertices and metastable stops can be a signature in the cosmologically favorable region. Let us conclude this discussion pointing out that the majority of the parameter space is not excluded by the current indirect detection upper bound on the R-Parity breaking parameter for a decaying DM gravitino and neither by the current excluded region for metastable particles. The LHC experiment will be able in the near future to explore all the parameter space up to the point where the stop NLSP is too heavy to be produced in sufficient numbers.

RPC and RPV stop decays at the LHC
In the previous section we have performed a decay-channel independent analysis, relying only on the presence or absence of a displaced vertex in the CMS detector. In case one of such displaced vertices is observed, the question arises as to which model is realized and if the decay is due to the RPC or RPV couplings. In Fig. 15 and Fig. 4.8 we show the Feynman diagrams for these two decay channels of stop.
From the point of view of the visible particles, both decays look similar, since they correspond to a two-body RPV decaỹ and the four-body RPC stop decay: i.e. both decays end in a charged lepton and a b-quark as visible particles, but they can be distinguished by the missing energy in the decay and the particle kinematics.
The different phase space of these decay channels can be observed in different kinematical variables related to the visible particles. In particular, we focus our attention on the three following observables: the antilepton transverse momentum P T , the transverse mass of the pair antileptonbottom M T and finally, the angle between the bottom and the antilepton momentum ϑ b , which are defined by the formulas The variables E , p = (p x , p y , p z ), E b , p b = (p bx , p by , p bz ) stand for the energy and the threemomentum of antilepton and bottom, respectively. In order to compare these kinematical quantities for the two decay channels of the stop, we simulate both processes with MadGraph 5 for a stop mass of mt = 800 GeV and the same stop decay rate of Γt = 2.02159 × 10 −10 GeV.
The transverse momentum distribution of the final antilepton for the two-body stop decay and the four-body stop decay are displayed in the left and the right panel of Fig. 17 respectively. Looking at the figure we see that the peak of the two-body distribution is located at a much larger transverse momentum than the peak of the four-body distribution, because in the two body decay the lepton takes away around half the rest energy of the stop. Note that the presence of a lepton or an antilepton in the final state is very useful since it provides a clear and robust signature in the detector, which significantly suppresses the large Standard Model background (e.g. QCD). In addition, the measurement of light leptons at CMS is more precise than that of jets. The transverse mass distribution of the final pair antilepton-bottom for the two-body stop decay and the four-body stop decay are plotted in the left and the right panel of Fig. 18 respectively. As we can see from the figure, the two distributions are very different. The two-body distribution has a peak and endpoint at a transverse mass of M T = 800 GeV corresponding to the stop mass, while the four-body distribution shows no peak and tends to be concentrated to much smaller range. The difference is due to the missing energy in the four-body decay, while the two-body decay distribution allows to infer the stop mass from the position of the end-point. The distribution of the angle between the final bottom momentum and the final antilepton momentum for the two-body stop decay and the four-body stop decay are displayed in the left and the right panel of Fig. 19 respectively. Observing the figure, we see that the two-body and the four-body distributions have no real peak. They are, instead, centered at different ranges of angle θ b , which are π/2 θ b 3π/4 for the two-body stop decay and θ b π/3 for the four-body distribution. This again is consistent with the fact that in the two-body decay the lepton and bottom are back-to-back in the (slowly-moving) stop rest system, while in the four-body decay they recoil against a non-negligible missing momentum. In general these distributions can be used not only to distinguish between the two RPC and RPV decays discussed here, but also to disentangle the signal from the SM background from top decays and from the case of a different model or LSP, e.g. to compare to the case of a right-handed sneutrino LSP [84] or to other models like the µνMSSM [27,85].

Background and coincidence counting
So far we have completely neglected the background of Standard Model and Supersymmetric particles, a clearly optimistic working hypothesis. The most important SM background comes from the top-pair production at LHC. In fact, if the top quark decays into W + and bottom b and at last, the weak boson W + decays into antilepton + and neutrino ν , we obtain in the final state the same visible particles as in the RPC and RPV stop decay. On the other hand, the top decays even before hadronizing and much faster that the stop discussed here and therefore one can avoid all this SM background just by requiring a displaced vertex. More difficult is to eliminate another source of background coming from bbZ → bb + − , where the b-decay happens naturally away from the primary vertex and the lepton tracks are mis-reconstructed, as originating away from the interaction point. Moreover also underlying events can give rise to particles pointing to a secondary vertex, faking the presence of a long-lived particle.
In general a good strategy to eliminate such kind of reducible background is to consider the  presence of two displaced vertices in the same event, both consistent with the same decay time.
Indeed we can see that in many cases one expects to have both stop and anti-stop to decay in the same part of the detector and give a clear signal for the production of two long-lived particles. We give in the table below the percentage of decays of the stop and anti-stop in the different detector parts.
We see from Tables 1 and 2, that even in the unfortunate case of lifetime around 10 −9 s, which is at the boundary of the metastable particle searches, we have quite a large statistics of coincident events in the pixel and tracker, reaching approximately 50% of the displaced vertex events, irrespectively of the stop mass. For longer lifetimes, the coincidence of two metastable particles in the same event takes over, reaching quickly a large statistics. We can therefore conclude that requiring two coincident events does not reduce the signal statistics substantially, while it would certainly suppress the background from misidentification or underlying events.

Conclusion
We have studied in this paper the reach of the LHC in models with a stop NLSP and gravitino LSP, both for the case of R-parity conservation and violation. In both cases we expect the stop to have a long lifetime leading to the possibility of displaced vertices or metastable particles at the LHC.
From the cosmological perspective, in case of high reheat temperature above the electroweak scale, the RPC scenario seems to prefer gravitino masses above 1 GeV and therefore lifetimes giving metastable stops, while for the RPV case only indirect detection of gravitino DM decay gives an upper bound on the RPV parameter , still consistent with stop decays within the detector.
We have performed a model-independent analysis of the two expected signals: displaced vertices in pixel or tracker detectors and metastable tracks. Our analysis is based on the MadGraph event generator, which allows us to simulate both the stop-anti-stop production and its subsequent decays. We assumed here that most of the other colored supersymmetric particles are much heavier than the stop and outside of the LHC reach and therefore we considered only the stop direct production, compute by MadGraph at LO and corrected to include NLO with a constant k-factor. For the decay length distribution in the detector we also devised a simple analytical estimate based on the decay formula and the stop momentum distribution, that matches very well the MadGraph data-points and allows for a much easier exploration of the parameter space.
We have seen that even for relatively long lifetimes a substantial number of decays can happen in the tracker or pixel parts of the CMS detector. They surprisingly present similar number of events and can both be exploited to measure the short lifetime region, while the metastable particle search extends the reach quite strongly at long lifetimes. We have therefore shown that the two search strategies are complementary and that both are needed in order to cover all the macroscopic stop lifetimes up to a certain stop mass. It is interesting to note that both searches, either displaced vertices or metastable states run out of steam at a similar value of the stop mass, where the production cross-section becomes too small. In particular we obtain, neglecting the background, in the most conservative case a mass reach up to 1400 (2000) GeV at LO and 1600 (2200) GeV at NLO for LHC at 14 TeV and with an integrated luminosity of 25f b −1 and 3000f b −1 respectively. As expected this reach is much larger than the minimal one for a metastable stau from Drell-Yan production as given e.g. in [31], but smaller to the reach expected for a metastable gluino. Of course a full detector simulation is needed to confirm our findings, but the prospects seem to be very favorable for a combined analysis of both displaced vertices and metastable particle searches as discussed here.
We have translated the LHC reach in the parameter space of two models with gravitino DM, either RPC or RPV. In the first model the region compatible with high reheating temperature leads only to the signal of HSCP and the observation of displaced vertices in the first run II fase with 25 f b −1 would directly point to T R ∼ 10 4 − 10 3 GeV. In the second model instead the search for displaced vertices allows to close the gap between metastable particle and indirect detection bounds at low stop masses. In case displaced vertices are seen, the visible particle kinematics allows to distinguish between RPC and RPV decays from the presence or absence of any missing energy. The distributions and the flavour of the final lepton, together with the value of the stop lifetime will be crucial to disentangle the particular model.
In general, due to the pair production of stop and anti-stop, we expect to see in a large fraction of the events the coincident presence of two displaced vertices or two metastable particles in the same event, allowing to disentangle the signal from reducible backgrounds connected to b-decays or misidentification of overlapping events.