Quasi-stable neutralinos at the LHC

We study supersymmetric extensions of the Standard Model with small R-parity and lepton number violating couplings which are naturally consistent with primordial nucleosynthesis, thermal leptogenesis and gravitino dark matter. We consider supergravity models where the gravitino is the lightest superparticle followed by a bino-like next-to-lightest superparticle (NLSP). Extending previous work we investigate in detail the sensitivity of LHC experiments to the R-parity breaking parameter zeta for various gluino and squark masses. We perform a simulation of signal and background events for the generic detector DELPHES for which we implement the finite NLSP decay length. We find that for gluino and squark masses accessible at the LHC, values of zeta can be probed which are one to two orders of magnitude smaller than the present upper bound obtained from astrophysics and cosmology.


Introduction
Supersymmetric extensions of the Standard Model with broken R-parity have a rich phenomenology [1][2][3][4][5]. In most models rather large R-parity violating couplings are considered, which lead to prompt decays of the lightest superparticle in the detector. In models where small R-parity violating interactions generate neutrino masses, macroscopic decay lengths up to 1 mm are obtained [6]. In the case of gauge mediated supersymmetry breaking, R-parity violating decays then compete with R-parity conserving decays where the final state contains a gravitino [7].
If the gravitino is the lightest superparticle, its decays are doubly suppressed by the Planck mass and the small R-parity breaking parameter. Hence, its lifetime typically exceeds the age of the universe by many orders of magnitude, so that it remains a viable dark matter candidate [8]. In the case of very small R-parity breaking couplings, as they occur if R-parity is spontaneously broken at the grand unification scale, one obtains a consistent cosmology including primordial nucleosynthesis, thermal leptogenesis and gravitino dark matter [9]. In the following we shall study the collider phenomenology of these models extending our previous work [10].
Decaying gravitino dark matter leads to a diffuse gamma-ray flux [8,9,11]. For a binolike NLSP the matrix elements for gravitino decay and NLSP decay are directly related [10]. Together with the lower bound on the gravitino lifetime, and the corresponding upper bound on the R-parity breaking parameter ζ, which is derived from the diffuse gamma-ray flux observed by the Fermi-LAT collaboration [12,13], one obtains a lower bound on the NLSP decay length. The estimates in [10] led to ζ 3 × 10 −8 , and in the more detailed analysis of gamma-ray lines in [14] a slightly stronger bound was found, ζ 2 × 10 −8 . For gravitino masses O(100 GeV) the upper bound becomes more stringent: ζ ∼ 10 −9 [10]. This range in the parameter ζ corresponds to NLSP decay lengths varying from O(50 cm) to O(500 m).
Large macroscopic decay lengths are of great help in the search for decaying NLSPs. This remains true if the decay length is larger than the size of the detector since a sizeable fraction of NLSPs may still decay inside the detector. This has been studied for neutral [15] as well as charged [16] NLSPs. Neutralino decay lengths varying from 0.1 mm to 100 m also arise in models with generalized gauge mediated supersymmetry breaking [17]. Alternatively, charged [18] and neutral [19] NLSP decays have been studied for models where the decay lengths are so small that no displaced vertices are observed and R-parity breaking Yukawa couplings determine the hierarchy of decay channels. In this case multilepton events, and their flavour structure, are of crucial importance.
The subject of this paper is a quantitative analysis of neutralino NLSP decays at the LHC in the case of very small R-parity breaking, ζ 3×10 −8 . The goal is the determination of the sensitivity in ζ for varying gluino and squark masses. To achieve this we perform a simulation for the generic detector DELPHES [20] of signal and background events. We focus on events with a clean signature: cascade processes with jets where one of the produced neutralino NLSPs decays into Z-boson and neutrino, with the subsequent decay of the Z-boson into a muon pair. This allows us to determine a conservative 5σ discovery range.
Finally, we estimate the discovery reach of the LHC if all NLSP decays are taken into account.
The paper is organized as follows. In Section 2 we recall the main ingredients of the considered model of bilinear R-parity breaking, as well as the dominant NLSP decays into gauge bosons and leptons. In Section 3 we discuss qualitative signatures of the events: NLSP βγ distribution, / p T spectrum, and the number of final state leptons. In Section 4, five representative benchmark points are defined and the simulation of signal and background events is described with emphasis on the muon reconstruction. The numerical results for the chosen NLSP decay channel Z(µ + µ − )ν are given in Section 5, together with estimates for the discovery reach if all NLSP decays are included. We conclude with a summary in Section 6.

Connecting neutralino and gravitino decays 2.1 Bilinear R-parity breaking
We consider supersymmetric extensions of the Standard Model with bilinear R-parity breaking (cf. [4,5]) as they are obtained if the spontaneous breaking of B−L, the difference of baryon and lepton number, is related to the spontaneous breaking of R-parity [9]. Mass mixing terms between lepton and Higgs fields then appear in the superpotential 1 , as well as the scalar potential induced by supersymmetry breaking, These mixing terms, together with the R-parity conserving superpotential the scalar mass terms and the standard SU(3) × SU(2) × U(1) Y gauge interactions define the supersymmetric standard model with bilinear R-parity breaking. Note that the Higgs mass terms m 2 u and m 2 d contain the contributions both from the superpotential (3) and the soft supersymmetry breaking terms. For simplicity, we have assumed flavour diagonal mass matrices in (4).
As discussed in [10], it is convenient to work in a basis of SU(2) doublets where the mass mixings µ i , B i and m 2 id in Eqs. (1) and (2) are traded for R-parity breaking Yukawa couplings. This can be achieved by field redefinitions: the standard rotation of the superfields H d and l i , followed by a non-supersymmetric rotation involving all scalar SU(2) doublets, where [10]. The R-parity breaking Yukawa terms contain couplings between gauginos, lepton doublets and Higgs doublets. After electroweak symmetry breaking, where we have defined Here g, g ′ and h e ij are the SU(2) and the U(1) Y gauge couplings and the charged lepton Yukawa couplings, respectively.
The diagonal mass terms together with the mixing terms in Eq. (7) represent the 7 × 7 neutralino mass matrix of the gauginos b, w 3 , the higgsinos h 0 u , h 0 d and the three neutrinos ν i , and also the 5 × 5 chargino mass matrix of gaugino, higgsino and charged leptons. Both mass matrices have to be diagonalized to obtain the CKM-type matrix elements of neutral (V (χ,ν) ), charged (V (χ,e) ) and supercurrents (U (χ,ν) ). For the standard supergravity mass spectrum with a bino-like NLSP χ 0 1 one obtains the R-parity breaking matrix elements (s 2β = 2s β c β ) [10]: 2 Our notation for gauge fields and left-handed gauginos reads: B µ , b etc. where the photino matrix element is defined as Note that the charged and neutral current matrix elements agree up to the isospin factor at leading order in m 2 1i LO /2.

Gravitino and NLSP decays
The partial width for gravitino decay into photon and neutrino (cf. Fig. 1) is given by [8] Inserting the matrix element (12) one obtains for the gravitino lifetime to leading order in m Z /µ [10]: where Based on the Fermi-LAT search for dark matter decaying into two photons, upper bounds on the R-parity breaking parameter ζ were derived: ζ 3×10 −8 [10] and ζ 2×10 −8 [14]. The parameter ζ also controls the lifetime of the bino-like NLSP χ 0 1 . For a NLSP mass larger than 100 GeV decays into charged lepton and W -boson or neutrino and Z-boson (c.f. Fig. 2) are dominant [21]. The partial decay widths read 1i LO are the charged and neutral current matrix elements at leading order, which are given in Eqs. (10) and (11), respectively. The function f W,Z is a phase space factor which becomes important for neutralino masses close to the lower bound of 100 GeV (cf. Fig. 3), For large NLSP masses, m χ 0 1 ≫ m Z , one has whereas in the region m χ 0 1 ≃ 100 GeV The total neutralino NLSP decay width is given by the sum which corresponds to the decay length [10] cτ χ 0 Note that both, the gravitino and the neutralino NLSP lifetimes are functions just of ζ and the masses, without any further parameters. This is the case for the standard supergravity mass spectrum with a bino-like NLSP. This direct connection between the gravitino and NLSP lifetimes is the basis of our analysis.

Decaying neutralino NLSP at the LHC
In this section we classify the main LHC signatures of decaying neutralinos. We show that for small values of the R-parity breaking parameter ζ usual SUSY searches are insufficient to find the signal. However, R-parity violation leads to new signals including striking secondary vertices at large distances from the primary interaction point.

Decay signatures
Consider for simplicity, the following cascade process: where q is a squark, g is a gluino, and j denotes a jet. The final state neutralinos decay in a secondary vertex into W -bosons and leptons as well as into Z-bosons and neutrinos. Fig. 4 shows an example of a decay cascade with muons in final state. The distance between the collision point and the secondary vertex depends on the decay width of the neutralino (17) and hence on the R-parity breaking parameter ζ. Table 1 summarizes the LHC signatures for sufficiently large values of ζ such that it is probable that both neutralinos decay inside of the tracker volume. All the signatures contain at least three jets from the antecedent supersymmetric decays. The signatures are classified according to the final states in the neutralino decays: leptonic signatures involving only leptons in the final state, semi-leptonic signatures involving at least two charged leptons and jets, single lepton signatures containing only one lepton, all-hadronic signatures where only jets accompanied by neutrinos are present, and finally invisible channels where both neutralinos decay solely to neutrinos. Additionally we single out channels having a considerable amount of missing transverse energy / one of the main features searched for in usual searches for new physics. Furthermore the channels labeled as opposite sign could be found in usual supersymmetry (SUSY) searches as they include a considerable amount of / E T , many jets and one isolated lepton pair with different signs. However, some searches remove events with muon pairs having invariant mass around the Z pole in order to dispose of Drell-Yan Z/γ * → ll processes. Note that in the model presented in this work this cut would lead to a suppression of the signal. Other leptonic and semi-leptonic channels also contain opposite-sign lepton pairs but only small amount of / E T and therefore they are not considered in the usual searches, (c.f. [22,23]). Neutralino decays lead also to signatures containing same-sign lepton pairs but since no / E T is present in these channels they are usually discarded in order to suppress various backgrounds [24].
If the value of ζ is rather small one of the neutralinos will decay outside of the detector leading to signatures with large amount of / E T as shown in Table 2. The leptonic decays of one of the neutralinos inside the detector lead to a perfect opposite-sign signature. As mentioned above this signature can be hidden if one rejects events where the invariant mass distribution of the lepton pair is in the range of the Z-boson mass. Another strategy is the search for single lepton events with large amount of missing transverse energy.
Thus the applicability and the reach of the usual SUSY searches applied to the model presented in this work depends crucially on the size of R-parity breaking. In order to further evaluate this statement we investigated a number of characteristic variables in supersymmetric events. The events were generated with PYTHIA as described in the next section, with the mSUGRA boundary conditions m1 /2 = m 0 = 270, tan β = 10, a 0 = 0, and µ > 0. R-parity violating neutralino decays were taken into account. Fig. 5 shows the distribution of the βγ factors of the neutralinos. This factor enters the formula for the neutralino decay length and one sees from the plot that analytic results in the literature, which have been computed with βγ = 1, are correct within one order of  Figure 5: βγ distribution of neutralinos at generator level for benchmark point HH27 (see Table 4). The number of neutralinos corresponds to twice the number of the events scaled to 10 fb −1 at √ s = 7 TeV.  Table 4) and different values of the R-parity breaking parameter ζ. Generator level / p T is defined as sum over the p T of i) neutralinos decaying outside of the detector (see Section 4.4) and ii) all neutrinos produced inside of the detector. The number of events is scaled to 10 fb −1 at √ s = 7 TeV. magnitude. The most important kinematic property connected with the neutralino decay length is the amount of missing transverse momentum / p T which is shown in Fig. 6 for different values of the R-parity violation parameter ζ. The missing transverse momentum was computed as the sum of the transverse momenta of all neutrinos produced in the detector before the hadronic calorimeter (r < 1800 mm, |z| < 3700 mm) and the transverse momenta of the neutralinos decaying outside the hadronic calorimeter. The / p T distribution of the R-parity conserving model ζ = 0 cannot be distinguished from the model with ζ = 1 × 10 −9 . However, the distribution is significantly different for ζ = 3 × 10 −8 since in this case most events have only very little missing transverse momentum due to early neutralino decays. This suggests that our model could only hardly be discovered in usual searches relying on / E T . A further analysis with full detector simulation is needed in order to properly evaluate the discovery potential of usual SUSY searches.
Another general feature of models with relatively large ζ is the large possible number of leptons in the final state, illustrated in Fig. 7. The generator level particles selected for this plot had to fulfill the criteria shown in Table 3 imposed in order to select leptons from hard processes which could be reconstructed in a realistic detector. The cuts on the vertex position represent a pessimistic estimate of the reconstruction efficiency (see Section 4.5).

Search strategies
As mentioned in the previous section one of the striking features of the presented model are events with secondary vertices and possibly many leptons in the final state. The search for a secondary vertex is crucial in order to ensure the R-parity violating nature of the decays. Possible search strategies can be optimized in order to find some of the channels described in Tables 1 and 2. It is remarkable that many channels allow for the full reconstruction of the neutralino mass: all decay chains including Z-bosons or hadronically decaying W -bosons. The reconstruction of the neutralino mass from the particles produced in the Z-boson decay depends crucially on the full reconstruction of the secondary vertex, which is beyond the   Table 3. The color code for the curves in both plots is given in Fig. 7b. The number of events is scaled to scope of this work 3 . This method of neutralino mass reconstruction works also in R-parity conserving models where the neutralino decays into Z-boson and gravitino [17]. For example, one promising search strategy working for all ζ values considered in this work is based on single lepton events with some number of hard jets and missing transverse energy larger than 90 GeV. After the preselection one could look for events where the lepton is coming from a secondary vertex and try to reconstruct the W -boson mass from a jet pair. In the final step one could try to reconstruct the neutralino mass from the jets selected in the previous step and the lepton. However such study depends crucially on the knowledge of the detector response in the case of late decaying particles. A neutralino can decay in various detector components and lead to unusual signals. Furthermore the mass resolution is limited by the uncertainty in the jet energy scale and by the uncertainty in the determination of the jet momentum direction.
We will focus our study on leptonic final states, which have a particularly clean signature, and reconstruct the Z-boson coming from a secondary vertex. We will use only muon and track objects for which we assume to have modeled a realistic detector response (see Section 4.5). A possible background for this search are cosmic muons leaving no track in the detector. It is important to note that one would miss the signal in this channel entirely if one imposes a cosmic muon veto which rejects all events with muon pairs having no associated tracks (c.f. [22]). particle transverse momentum pseudorapidity vertex position electron p T > 7 GeV |η| < 2.5 r < 400 mm |z| < 1300 mm muon p T > 6 GeV |η| < 2.5 r < 4000 mm |z| < 6000 mm Table 3: Cuts for the generator level particle selection for the study of particle multiplicity.

Simulation of signal and background
In this section we define a set of representative points in the parameter space of our model and describe the generation of the signal and dominant standard model (SM) background samples. In particular we examine the simulation of detector effects using the generic detector simulation DELPHES 1.8 [20] on signal and background in the presence of secondary vertices.

Benchmark points
A typical set of boundary conditions for the supersymmetry breaking parameters of the MSSM at the grand unification scale is given by equal scalar and gaugino masses, m 0 = m1 /2 . These boundary conditions lead to a bino-like neutralino χ 0 1 as NLSP. We choose a representative value of tan β and set the scalar trilinear couplings to zero, Thus the universal gaugino mass remains the only independent supersymmetry breaking parameter which will be varied in the present study. Electroweak precision tests lead to lower bounds on the supersymmetric particle spectrum. The LEP lower bound for the lightest Higgs at tan β = 10 is m H 95 GeV [25]. In the present study the lightest superparticle spectrum corresponds to the choice m 0 = m1 /2 = 270 GeV (HH27). At this benchmark point the NLSP is a neutralino with mass m χ 0 1 = 105.8 GeV and the lightest Higgs boson has a mass m h = 110.4 GeV. In order to probe the region of gluino and squark masses accessible at the LHC [26] we increase the gaugino mass parameter in four steps: m1 /2 = 350, 500, 650, 800 GeV. Some particle masses at these points are shown in Table 4. For the different benchmark points the production cross sections, calculated with PROSPINO 2.1 [27] at √ s = 7 TeV, are given in Table 5.

Major backgrounds
Neutralino decays always have W -and Z-bosons in the final state (c.f. Fig. 2). In our study we focus on the reconstruction of Z-boson decays to muon pairs. Therefore we only  consider SM backgrounds which lead to at least two muons in the final state originating from W -or Z-bosons: • tt production: W -bosons from top quark decays.
• Z production • Di-boson production (W W , W Z, ZZ) • Tri-boson production (W W W , W W Z, ZZW , ZZZ) Table 6a gives an overview of the background samples used in our analysis. We assume that pure QCD background can be efficiently suppressed in multi-lepton final states with high transverse momentum, particularly after imposing lepton isolation criteria (c.f. [23,28]).

Event simulation
All Monte Carlo samples were generated using parton distribution functions given by CTEQ6L1 [33]. For the simulation of the background we used MADGRAPH 4.4.44 [34] interfaced with PYTHIA 6.4.22 [35]. Our simulation of the signal events relied on the following procedure. First, supersymmetric mass spectra were calculated with a modified version of SOFTSUSY 3.1.5 [36] assuming mSUGRA boundary conditions and R-parity conservation. The latter assumption is justified due to the tiny amount of R-parity breaking in our model. The SOFTSUSY version was modified in order to produce additionally to the spectrum the R-parity violating neutralino decay width and branching ratios according to Equation (17). The SOFTSUSY mass spectra were fed into SDECAY [37] via the MADGRAPH homepage [38] in order to calculate the decay widths of the SUSY particles (besides the neutralino LSP). In the next step neutralino decay information was included into the SDECAY output. The signal process (production of g g, g q, q q and q q) was simulated with MADGRAPH and then given to PYTHIA for computation of all subsequent decays according to the SDECAY output as well as for parton showering and hadronization.  The generic detector simulation DELPHES, tuned to the CMS detector, was used in order to account for effects of event reconstruction at the detector level. However, DELPHES describes the detector geometry solely in terms of angular variables, i.e. the detector is stretched infinitely in the radial direction. This approximation is sufficient for most studies involving prompt decays but is untenable in the case of late decaying particles. We overcome this obstacle by adding vertex information from particles at the generator level to objects at the detector level. Usually, this information is provided by the detector simulation. Our procedure is described in detail in the following section. We emphasize that a full detector simulation, which includes vertex reconstruction, needs to be done to improve our analysis.

Muon reconstruction process
Particles produced in the late decay of the neutralino will not be properly reconstructed in a real detector if the position of their vertex is beyond or even within the crucial detector component responsible for the respective identification. For example, an electron produced inside of the electromagnetic calorimeter will leave no track in the tracker and will therefore be identified as a photon or jet. In order to simulate the detector response to such events we use a detector geometry in the (r, z) coordinates, which is inspired by the CMS detector at the LHC (see Fig. 8). The angular position of the detector components is given by the CMS tune of DELPHES.
In order to be as conservative as possible we only use muon and track objects for the present analysis, since these objects allow a simple simulation of detection efficiency losses due to the finite size of the detector. Namely, we assume that a muon can be reconstructed ζ events HH27 = 1 × 10 −9 22280 1 × 10 −9 222800 HH35 = 1 × 10 −9 10000 1 × 10 −9 100000 HH50 = 1 × 10 −9 10000 1 × 10 −9 100000 HH65 = 1 × 10 −9 10000 1 × 10 −9 100000 HH80 all ζ 10000 (b) Samples of signal events for different benchmark points (see Table 4) and ζ = α × 10 −9 (α = 0.1, 0.5, 1, 5, 10, 20, 30). as long as its vertex is in front of the muon chambers, and analogously a track can be reconstructed if it originates approximately in the first third of the tracker (This region is called pixel detector in Fig. 8). For the matching between generator level particles and objects reconstructed by DELPHES we use the distance in pseudorapidity η and azimuthal angle φ, defined as ∆R = (∆φ) 2 + (∆η) 2 . In the following we will call generator level muons, produced by PYTHIA, GenMuons, muons reconstructed initially by DELPHES muon candidates, and track objects reconstructed by DELPHES RecoTracks. Only GenMuons and RecoTracks have the coordinates of their vertex.
First, we perform the following p T cuts on muon candidates and RecoTracks: • p T (µ) > 20 GeV, These cuts are guided by our SUSY search strategy (c.f. Section 5), since we expect that muons coming from Z-boson decay have high p T , and a sufficiently high p T cut can effectively suppress QCD fake leptons. Furthermore, DELPHES itself reconstructs only muons with p T above 10 GeV. Additionally, these cuts were optimized in order to get a realistic muon reconstruction efficiency (see Section 4.5).
In the second step vertex information is added to the muon candidates by matching with GenMuons: • A GenMuon is selected for matching with muon candidates if its vertex lies in front of the muon system : r µ = √ x 2 + y 2 < 4000 mm, |z µ | < 6000 mm (see Fig. 8). • The ∆R distance between each selected GenMuon and all muon candidates is computed.
• A GenMuon vertex is added to the muon candidate closest in ∆R, if ∆R < 0.1 and GenMuon and muon candidate have the same charge.
• Muon candidates with added vertex information are called RecoMuons.
In the final step, muons with or without signal in the tracker are distinguished: • A RecoTrack is selected for matching with RecoMuons if the track vertex lies in the following range: r T < 400 mm , |z T | < 1300 mm.
• Each selected RecoTrack is matched with the RecoMuon closest in ∆R, if ∆R < 0.1.

• Matched RecoTracks and RecoMuons are called tracker muons.
RecoMuons which cannot be matched with RecoTracks are called chamber muons. Each RecoMuon is therefore either a tracker muon or a chamber muon.
After the reconstruction procedure one is left with two kinds of muon objects: (i) chamber muons which have no track in the tracker and are therefore reconstructed solely by the muon chambers, and (ii) tracker muons which have a track. The muon reconstruction process is depicted in Fig. 9. The ∆R matching condition has been optimized in order to get a realistic muon reconstruction efficiency (see next section).

Muon reconstruction efficiency
In order to test our method of obtaining physically sensible objects we compute the muon reconstruction efficiency in the following way: • Muons are created as described above.
• GenMuons are matched with RecoMuons without any constraints on the position of the GenMuon vertex.
• The number of successfully matched objects is compared binwise (in bins of r and |z|) with the number of all GenMuons.
The second condition is necessary in order to see whether the assignment between Reco-Muons and GenMuons is correct. Since the matching procedure only relies on angular variables, it is possible that a RecoMuon originally matched with a GenMuon created in front of the muon chamber belongs in fact (i.e. has smaller angular distance) to a Gen-Muon coming from a decay inside the muon chamber or even outside of the detector. Such wrong matchings would be seen in the efficiency plot as efficiencies not equal to zero in regions where muons could not be detected by the detector defined above (r µ > 4000 mm, |z µ | > 6000 mm). Fig. 10 shows the computed muon efficiency in bins of r and |z|. As expected one sees a sharp decline in efficiency in the r plot at r µ = 4000 mm, where the hard cut applies. The decline in the z plot is gradually, since physical particles have to fulfill both r and |z| criteria. The particles originating at small values of r and large values of |z| are not reconstructed due to the limited pseudorapidity coverage of the muon detector. The efficiency stays at zero beyond r = 4000 mm and |z| = 6000 mm as expected, confirming our method of muon reconstruction. We expect that the computed muon efficiency agrees within 15 % with efficiencies of present LHC detectors including losses due to muon-jet separation requirements.

Search for the neutralino decay χ 0 1 → Zν
As described in Section 3.2 our study is focused on the channel χ 0 1 → Zν → µ + µ − ν. This channel possesses certain physical and technical advantages. On the physical side reliable muon identification is possible already in the early stage of the LHC data taking and one can assume that QCD background can hardly fake two muons at the same time. Furthermore this signal leads to spectacular events and has no easily identifiable SM background at all, as shown in this section. Additionally, the muon chamber is the detector component which is farthermost away from the primary vertex and hence one can expect that it will be possible to detect a significant number of clean late time decays even for very small R-parity breaking. On the technical side, muons seem to be the simplest objects for which a realistic detector response can be modeled within DELPHES (see Section 4.4), due to the limitations of this simulation in the presence of secondary vertices.
The spectacular feature of this signal are opposite sign muon pairs with invariant mass close to the Z-boson mass, which have either associated tracks in the tracker with clearly visible secondary vertices or no associated tracks at all. Such muon pairs cannot be generated by usual SM background as will be shown in the following. However, a similar signal can arise from cosmic muons traversing the detector. We could not create a Monte Carlo background sample for cosmic muons, and we simply assume that such background can be suppressed by use of the full timing information of the event: cosmic muons will first cause a signal in the muon chamber which is closest to the ceiling of the experimental hall followed by a signal in the opposite direction.
An intrinsic background for the presented search are muon pairs from R-parity violating decays, where one muon is coming from the W -boson decay while the other muon is coming either from the neutralino decay into the W -boson in either of the two branches or from the W -or Z-boson decay in the second branch (c.f. Fig. 4). This background can be suppressed if one has access to the corresponding tracks by demanding that both of them originate from the same vertex. In the case of muons without tracks this background is irreducible. However it belongs itself to the signal one is looking for.

Event selection
In order to find the signal we now employ a series of simple cuts on the reconstructed objects (muons, tracker muons, and chamber muons). As we will see, already with an integrated luminosity of only 1 fb −1 at √ s = 7 TeV a discovery of the benchmark scenario HH27 with ζ = 3 × 10 −8 is possible. First, we perform a selection cut on the number of muons in the event: • N(muons) ≥ 2.
We define two event classes depending on the number of tracker muons: • Class 1: the event contains at least two tracker muons N(tracker muons) ≥ 2.
From the description of the signal presented above, we implement additionally two sets of cuts depending on the class of the event. The cuts for Class 1 events are: • All possible invariant masses of opposite sign tracker muons are computed. The event passes the cut if at least one invariant mass is in the range of the Z-boson mass: 80 GeV < M µ + µ − < 100 GeV. If the event contains more than one appropriate combination of the tracker muons then the muons from the combination with invariant mass closest to the Z-boson mass are selected for further analysis.
• d(Vertex) > 5 mm: Each of the tracks associated with the two selected tracker muons should have a vertex which is further than 5 mm away from the primary vertex. This value is approximately one order of magnitude larger than the current resolution of the inner tracker (c.f. [28,39]).
• ∆d(Vertex) ij < 5 mm: The distance between the two track vertices should be less than 5 mm.
• If the event fails one of the cuts it is classified as a Class 2 event.
The cuts for Class 2 events are: • N(chamber muons) ≥ 2: If an event has less than two tracker muons than it should have at least two chamber muons.
• All possible invariant masses of opposite sign chamber muons are computed. An event passes the cut if at least one invariant mass is in the range of the Z-boson mass: 80 GeV < M µ + µ − < 100 GeV.
Since each Class 1 event is classified as a Class 2 event if it fails one of the cuts, no signal event is discarded because of the presence of muons with tracks not coming from neutralino decay.  (ζ = 1 × 10 −9 ) corresponds to an integrated luminosity of 10 fb −1 (≈ 100 fb −1 ).
Most events will fall into the second class. The analysis is then very simple and amounts to the search for events with muons without associated track in which the invariant mass of a muon pair lies in the Z-boson mass range. The cut flow is given in Table 7. As expected, no background events survived the cuts, since no standard model process should produce secondary vertices so far away from the primary interaction point. Although our background estimate has an uncertainty due to the limited statistics, we assume on physical grounds that no background events will pass the cuts if we increase the number of simulated events. However, the major uncertainty in this study, the number of the background events from cosmic muons, cannot be estimated with the present software. Therefore a full fledged analysis with full detector simulation which takes into account the cosmic muon background is needed. In the following we assume that this background can be efficiently suppressed with the full timing information of the event as described in the introduction to Section 5. Furthermore, we only estimate the systematic uncertainty due to the background and neglect statistical errors and the uncertainty of the muon reconstruction efficiency.
The significance of the signal is computed with the profile likelihood method [40] incorporated in the SIGCALC code [41]. We assume an integrated luminosity of 1 fb −1 at √ s = 7 TeV LHC and a ten times higher Monte Carlo luminosity L M C = N b /σ b = 10 fb −1 for all the background events. At this integrated luminosity 17 signal events and no background events survive the cuts, which corresponds to a significance Z P L = 9.03. Instead, if one makes the pessimistic estimate that 1 background event from the cosmic muons passes the cuts one finds a significance Z P L = 6.39. Therefore we conclude that at the benchmark point HH27 with ζ = 3 × 10 −8 , R-parity breaking neutralino decays can be discovered with the first inverse femtobarn of LHC data.

Discovery reach at the LHC
In the previous section we have studied in detail the benchmark point HH27: m1 /2 = m 0 = 270 GeV, which yields the rather small superparticle masses m χ 0 1 = 106 GeV, mg ≃ 660 GeV and mq ≃ 650 GeV for the light quark flavors (cf. Table 4). From the decay rates given in Section 2.2 and the phase space factors shown in Fig. 3 one obtains for decay length and branching ratio into Z-boson/neutrino final states: Based on the production cross sections listed in Table 5 an integrated luminosity L = 10 fb −1 yields about 22000 events and therefore 44000 NLSPs. We have studied this benchmark point for two different values of the R-parity breaking parameter: ζ = 3 × 10 −8 and ζ = 1 × 10 −9 . For the larger value of ζ one has cτ χ 0 1 ≃ 50 cm. Hence, essentially all neutralinos decay inside the detector, most of them close to the origin. The spacial distribution of secondary vertices is displayed in the contour plot Fig. 11. Using BR(Z → µ + µ − ) ≃ 0.034 and the branching ratio given in Eq. all cuts (cf. Table 7). The locations of the secondary vertices of these events are shown in Fig. 12.
For the smaller value of the R-parity breaking parameter, ζ = 1×10 −9 , the decay length increases to cτ χ 0 1 ≃ 450 m. Now most neutralino NLSPs decay outside the detector. This is apparent from Fig. 13 where the total number of decays in the different subvolumina of the detector are given. Compared to ζ = 3 × 10 −8 , the number of decays inside the detector is smaller by a factor ∼ 200, which roughly corresponds to the ratio of the decay lengths, as suggested in [15].
According to the simulation described in the previous section, for ζ = 1 × 10 −9 an integrated luminosity of 100 fb −1 is needed to obtain 2 decays χ 0 1 → Zν → µ + µ − ν, inside the detector. This number is very small compared to the total number of about 2000 decays in the detector, which is a consequence of the tiny branching ratio into the chosen specific final state. It is likely that a substantially larger fraction of the events can be used in the search for a decaying neutralino. In [15] it has been argued that already 10 χ 0 1 decays inside the detector may be sufficient for the discovery of a decaying NLSP, which would require an integrated luminosity of only 1 fb −1 . It remains to be seen whether for events with a secondary vertex and jets, signal and background can be sufficiently well separated.
Let us now consider the benchmark point HH50: m1 /2 = m 0 = 500 GeV, which implies the heavier superparticle masses m χ 0 1 = 206 GeV and mg ≃ mq ≃ 1200 GeV for the light quark flavours (cf. Table 4). The phase space suppression is now negligible, and one obtains for decay length and branching ratio into Z-boson/neutrino final states: The total production cross section for these heavier gluino/squark pairs is about two orders of magnitude smaller (cf. Table 5), and therefore an integrated luminosity L = 10 fb −1 only yields 460 NLSPs. We have studied this benchmark point again for the two different values of the Rparity breaking parameter ζ = 3 × 10 −8 and ζ = 1 × 10 −9 . For the larger value of ζ one has cτ χ 0 1 ≃ 10 cm, and essentially all neutralinos decay inside the detector. The branching ratio into the considered final state is now somewhat larger, BR(χ 0 1 → Zν → µ + µ − ν) ≃ 0.01, so that one expects about 4 events with this final state, which is consistent with our simulation. Hence, for this larger value of the R-parity breaking parameter and this benchmark point, the discovery of a decaying NLSP appears feasible already in the early phase of the LHC.
For ζ = 1 × 10 −9 , the decay length is cτ χ 0 1 ≃ 80 m and most neutralino NLSPs decay outside the detector. The spacial distribution of secondary vertices inside the detector, in total 12 for 10 fb −1 , is shown in Fig. 14. Due to the 1 % branching ratio into the Zν → µ + µ − ν final state one then estimates that 1000 fb −1 will be needed for a discovery, which is consistent with our simulation.
In Fig. 15 we have summarized the results of our simulations for the decay chain χ 0 1 → Zν with Z → µ + µ − . The benchmark points HH27-HH80 correspond to gluino and squark masses ranging from 650 GeV to 1800 GeV (cf. Table 4). The bands reflect the different number of events required for a 5σ discovery depending on the simulated background. The central value corresponds to 6 signal events (with luminosity L) with no background events for a simulated luminosity of 10 × L; the lower (upper) boundary represents 3 (13) signal events (with luminosity L) with no (1) background event for a simulated luminosity of 100 × L (10 × L). We conclude that with 10 fb −1 a 5σ discovery of a quasi-stable neutralino is possible for squark and gluino masses of 830 GeV (cf. HH35) and an R-parity breaking parameter ζ = 3 × 10 −9 , which is one order of magnitude smaller than the present astrophysical bound [10,14].
We expect that the sensitivity in the parameter ζ can be significantly improved if also neutralino decays with jets are taken into account. Fig. 16 represents an estimate of the discovery reach for quasi-stable neutralino NLSPs at the LHC, assuming 10-20 decays inside the detector (cf. [15]). The parameter space, which can be probed, is now significantly extended. As an example, with 10 fb −1 and squark and gluino masses of 830 GeV (cf. HH35), one is now sensitive to ζ = 3×10 −10 , which lies two orders of magnitude below the present astrophysical bound. Correspondingly, for heavier gluinos and squarks, mg ≃ mq ≃ 1480 GeV (cf. HH65), one can probe values of the R-parity breaking parameter down to ζ = 3 × 10 −9 .   Figure 16: Estimate of the 5σ discovery reach in ζ for quasi-stable neutralino NLSPs at the LHC; the lower (upper) boundary of the bands corresponds to 10 (20) decays inside the detector. The different bench mark points correspond to gluino and squark masses between 650 GeV and 1800 GeV.

Summary and conclusion
We have studied the decays of quasi-stable neutralinos in supersymmetric extensions of the Standard Model with small R-parity and lepton number violating couplings which are consistent with primordial nucleosynthesis, thermal leptogenesis and gravitino dark matter.
As representative examples with neutralino NLSP we have considered five benchmark points, HH27-HH80, with gluino and squark masses ranging from 650 GeV to 1800 GeV. For the considered benchmark points the neutralino is bino-like, and the R-parity breaking parameter ζ, which governs the neutralino decays χ 0 1 → W ± l ∓ and χ 0 1 → Zν, is directly related to the partial gravitino decay width into photon and neutrino.
The main goal of the present work has been to determine the range of the parameter ζ which can be probed at the LHC, for varying superparticle masses. As a conservative starting point, we have focused on events with a clean signature: cascade processes with jets where one of the produced neutralino NLSPs decays into Z-boson and neutrino, with a subsequent decay of the Z-boson into a muon pair.
In Section 3 we have studied the qualitative signatures of these events, the βγ distribution of the produced NLSPs, the / p T -spectrum and the number of leptons in the final state. A detailed simulation of signal and background events for the generic detector DELPHES (with CMS tune) has been described in Section 4, with emphasis on the reconstruction of muons. The crucial new element of our analysis has been the implementation of the finite NLSP decay length.
The numerical results of our simulations have been summarized in Section 5. It is very instructive to look at the spacial distribution of decay vertices which have been given for different benchmark points and different values of ζ. One can clearly see that the sensitivity extends from decay lengths O(50 cm), with NLSP decays mostly inside the detector, to values O(500 m) where almost all NLSPs decay outside the detector.
The results for the discovery reach for quasi-stable neutralino NLSPs roughly agree with the simple estimates which one obtains from the branching ratios into the Z(µ + µ − )ν final state together with the assumption that these events are background free. It is remarkable that already with 10 fb −1 a 5σ discovery is possible for squark and gluino masses of 830 GeV and an R-parity breaking parameter ζ = 3 × 10 −9 , which is one order of magnitude smaller than the present astrophysical bound.
It is likely that the severe restriction of our simulation to Z(µ + µ − )ν final states can be relaxed and that a much larger fraction of events can be used for the analysis. Assuming optimistically that 10 decays inside the detector are sufficient for a discovery, one would be sensitive to values of ζ down to 3 × 10 −10 for squark and gluino masses of 830 GeV and 10 fb −1 . Alternatively, for gluino and squark masses of 1480 GeV one could probe values of the R-parity breaking parameter down to 3×10 −9 . For the same quark masses a luminosity of 100 fb −1 would improve the sensitivity in ζ by a factor of three. We conclude that for gluino and squark masses accessible at the LHC values of the R-parity breaking parameter can be probed which are far below the present upper bounds obtained from astrophysics and cosmology.