Searching for light long-lived neutralinos at Super-Kamiokande

Light neutralinos could be copiously produced from the decays of mesons generated in cosmic-ray air showers. These neutralinos can be long-lived particles in the context of R-parity violating (RPV) supersymmetric models, implying that they could be capable of reaching the surface of the earth and decay within the instrumental volume of large neutrino detectors. In this letter, we use atmospheric neutrino data from the Super-Kamiokande experiment to derive novel constraints for the RPV couplings involved in the production of long-lived light neutralinos from the decays of charged $D$-mesons and kaons. Our results highlight the potential of neutrino detectors to search for long-lived particles, by demonstrating that it is possible to explore regions of parameter space that are not yet constrained by any fixed-target nor collider experiments.


I. INTRODUCTION
The discovery of the Higgs boson at the Large Hadron Collider (LHC) in 2012 [1,2] provides not only conclusive evidence of the Standard Model (SM), but also consolidates the hierarchy problem as one of the main theoretical puzzles in modern physics [3]. In this context, supersymmetry (SUSY) remains as one of the most compelling possibilities to address this problem [4,5]. At the same time, supersymmetry also provides a rich and complex phenomenology which has lead to an intensive search program at collider experiments [6]. Conventional SUSY theories assume a discrete symmetry called R-parity, which avoids conflict with experimental data on the non-observation of baryon and lepton number violating processes, such as proton decay [7] and neutrinoless double beta decay [8]. Within the context of Rparity conserving SUSY theories, the lightest neutralino as the lightest supersymmetric particle (LSP) provides a natural candidate for fermionic dark matter because of its stability and lack of electromagnetic interactions (see [9,10] for seminal articles and [11] for a review). Yet, it is possible to assume R-parity violation (RPV) [12] while respecting the bounds on the proton lifetime, as long as baryon number or lepton number are preserved. In such a case, the neutralino LSP is no longer stable and can decay into SM particles. The smallness of the R-parity violating couplings can make the decay macroscopic, making the neutralino,χ 0 1 , a long-lived particle (LLP) [13]. This implies that neutralinos can leave a variety of exotic signatures at colliders such as a displaced vertex (reviews on possible long-lived particle signals can be found in [14,15]). On the other hand, unlike the strong interacting sparticles whose masses have a lower bound around 1 TeV [16][17][18][19][20][21][22][23][24], the neutralino mass is less constrained, and can in principle be massless [25]. For a fraction of the R-parity violating model parameter space, the AT-LAS and CMS experimental collaborations at the LHC have searched for a long-livedχ 0 1 with a mass O(100) GeV with a displaced vertex signature [26][27][28], with null results. The most up-to-date constraint comes from AT-LAS with leptonic displaced decays, excluding long-lived neutralinos between 50 − 500 GeV [27].
Phenomenological prospects demonstrate that in extended regions in the R-parity violating couplings -leading to decays with different leptonic or hadronic final state particles -and lighter masses below 100 GeV, a long-livedχ 0 1 has the potential to be discovered at dedicated LLP experimental facilities that could operate at the LHC, such as FASER, MATHUSLA, CODEX-b and AL3X [29][30][31][32][33][34], future colliders at the intensity frontier [34] or beam-dump experiments as SHiP [32]. In particular, the absence of experimental constraints in the range from few MeV to O(1) GeV makes it a good candidate to be studied in scenarios were they can be produced from the mesons that are abundantly created in both colliders and cosmic-ray air showers.
So far, prospects for light, long-lived neutralinos have been performed from the decays of D and B mesons in references [30,32,[35][36][37], from Z boson decays in [29,31,34,38], and from the decay of τ leptons in Belle-II [33]. Most of these searches are based in experiments that are either in construction like FASER [39], or even in earlier stages as they are subject to funding/approval (e.g. SHiP [40], MATH-USLA [41] and others [42,43]). In contrast to this situation, large Cherenkov based neutrino detectors such as IceCube [44] or Super-Kamiokande (SK) [45] are already built and have been taking data for years, which can be used to search for LLPs, as the decay of these particles would generate a signal that is indistinguishable from the Cherenkov radiation measured in association with the neutrino interactions in the medium.
Searches for long-lived particles at large neutrino detectors were considered in the literature in references [46][47][48][49][50][51]. Of these studies, reference [49] used public data from IceCube and Super-Kamiokande to search for long-lived particles produced in atmospheric showers, using detailed numerical simulations. The results found demonstrated that the atmospheric neutrino data from the Super-Kamiokande experiment can be used to place stringent constraints to models predicting LLPs (see also reference [50]). The main reason for this is because the atmospheric neutrino data focused on the low energy regime, where the flux of the atmospheric particles peaks. In this work we build upon the general results of [49], and apply the same strategy to search for neutralinos produced from meson decays, including an appropriate treatment of the uncertainties originated from the hadronic interaction models used to simulate the production of mesons in the shower.
We consider the possibility of searching for neutralinos produced in the decays of D-mesons and kaons that are generated in the sky, when highly energetic cosmic rays collide with the upper layers of the atmosphere. The former process allows us to compare with the previously mentioned studies in the literature, particularly references [30] and [32], where the sensitivity reach at collider experiments was estimated for selected benchmarks. On the other hand, neutralino production from kaon decay was considered long before in the context of spontaneously broken supersymmetry [52]. In this work, we explore a novel RPV channel through kaon production that is particularly well suited for our setup because the production of kaons in the atmosphere greatly exceeds the one from D-mesons [53].
The rest of the paper is organized as follows. In Sec. II we summarize the relevant phenomenological aspects of the RPV theory, and describe the neutralino production from the decay of pseudoscalar mesons generated in cosmic-ray air showers. In Sec. III, we describe the expected signal from the visible decay of the neutralino in the Super-Kamiokande detector. The main results obtained for the two benchmark scenarios considered in this work are presented in section Sec. IV, along with the current best constraints from previous studies. Finally, we draw our main conclusions in section Sec. V.

A. Neutralinos in RPV
The simplest realistic realization of a supersymmetric model is the so-called minimal supersymmetric standard model (MSSM) [54]. As mentioned above, in this model a Z 2 symmetry called R-parity is imposed to avoid proton decay. However, it is still possible to break R-parity and have a stable proton by imposing a different discrete symmetry like the baryon triality B 3 symmetry [55]. Without imposing R-parity, the most general Lagrangian that respects gauge and space-time symmetries contains the term [56] where the hatted symbols denote gauge multiplets of superfields that include leptons i and sleptons˜ i in L i , left-handed quarks q jL and squarksq jL in Q j , and righthanded antiquarks d kR and antisquarksd * kR in D c k . The family indices i, j and k run from 1 to 3, leading to 27 independent parameters. This semileptonic operator violates lepton number by one unit, and can generate contributions to the mass and magnetic moment of neutrinos [57], neutrinoless double beta decay [58,59], and in the parameter space of our interest, neutralino production from meson decays [60][61][62]. Comprehensive reviews of RPV and its phenomenological implications can be found in [12,[63][64][65].
In this work we focus on the trilinear couplings λ 121 and λ 112 , which allow to produce neutralinos from the decays of D-mesons and kaons, respectively. These decays are accompanied always with an electron, as detailed in our two benchmarks in table II. The introduction of RPV parameters also allows neutralino two-body decays into a meson and a lepton. In our first benchmark, we study neutralinos with masses in the range from the mass of the kaon to the mass of the D-meson, while in the second benchmark case we restrict the neutralino mass to the interval between the mass of the pion and the mass of the kaon. Therefore, for the second benchmark, we include the trilinear λ 111 , as it is the only coupling associated with pions, and is therefore necessary to make the neutralino decay when its mass is lower than the kaon mass. The relevant Feynman diagrams involved in the production of neutralinos from these mesons are shown in figure 1.
To calculate the decay rate of the mesons to neutralinos, we closely follow reference [32] and consider the effective interactions involving the relevant mesons, leptons and neutralinos. Assuming that the sfermion masses involved are degenerate and large enough so that they can be integrated out, the relevant decay rates are and Γ(K ± →χ 0 1 + e ± ) = where K(x, y, z) ≡ x 2 + y 2 + z 2 − 2xy − 2xz − 2yz, is the Källén function [66] and f D , f K are the D-meson and kaon decay constants, respectively. The effective couplings G S ijk are given by where mf represents the common value that we assume for the masses of sleptons and squarks, which are the sfermions involved in neutralino production (details can be found in reference [32]). In equation (4), θ W ≡ tan −1 g 1 /g 2 is the weak mixing angle [21], and g 1 , g 2 are the gauge coupling constants associated to U(1) Y and SU(2) L , respectively. The symbols m K and m D denote the masses of the K ± and D ± mesons, and m c , m d , m u and m s are the masses of the charm, down, up and strange quarks, respectively. In our calculations, we use the values f K 156 MeV and f D 213 MeV [21].

B. Neutralinos from D-meson and kaon decays in atmospheric showers
Cosmic rays hitting our atmosphere provide us with a beam of protons (and other species) that is constantly switched on. A single cosmic-ray can produce an extensive cascade of charged particles and radiation called an air shower [67][68][69]. Mesons in the shower, including D-mesons and kaons, decay to charged leptons and neutrinos, among other particles. The flux of leptons can be measured both at the surface of the earth and in underground experiments, while their careful reconstruction allows to estimate the expected mesonic contributions to the spectrum [53,[70][71][72][73][74][75][76][77]. Similarly as with the case of neutrinos, we simulate the production of light neutralinos in the shower by solving the cascade equation involving only source terms from meson decays [74] where the sum runs over all possible mesons that can decay to neutralinos when a given trilinear coupling λ ijk is switched on. In equation (5), ρ is the density of the atmosphere at a column depth X, and λ M ≡ γ M β M cτ M is the decay length of the meson, which includes the boost factor γ M β M and its proper lifetime τ M . The differential production rate of mesons in the shower per unit of solid angle is given by dΦ M dE M dΩ . The number of neutralinos with energies between Eχ0 1 and Eχ0 1 + dEχ0 1 produced in the decay of the meson M is given by dn dEχ0 1 . For two-body decays, this last quantity is given by where Br(M →χ 0 1 +e) is the branching fraction of meson decays to neutralinos and p M is the meson momentum. In equation (5), both the density of the atmosphere and the differential production rate of mesons are extracted using the Matrix Cascade Equation (MCEq) software package [75,77]. Here, we choose the NRLMSISE-00 atmospheric model [78], while for the hadronic interaction model we consider the different event generators which are updated with LHC data (see for instance, references [79,80]). In particular, we focus on the SYBILL-2.3 [53], QGSJET-II-04 [81], EPOS-LHC [82], and DPMJET-III [83] models. These models might yield non-negligible differences in the production rate of mesons and can be a relevant source of uncertainty in our calculations 1 .
As an illustration of our numerical simulations, we show in figure 2 the production rate for a 0.31 GeV neutralino, at a height h of 15.51 km and with a vertical direction of 25.84º degrees from the zenith, produced in the decay of both D-mesons and kaons. The spread of the orange band in the plot covers the most pessimistic and optimist cases given the uncertainties associated with the election of an event generator. Note that there is an important difference in the amount of neutralinos produced Neutralino production diagrams. We display the specific case where production occurs via the decay of a charged meson M + kj composed of quarks d C k and uj. For the cases of our interest (see table II), we set the flavour indices of the initial states in the diagram such that M + kj corresponds to D + (k = 1, j = 2) for benchmark B1 and K + (k = 2, j = 1) for benchmark B2.
from the decay of D-mesons and kaons, as the latter are expected to be orders of magnitude more abundant in the atmosphere [53]. Finally, it is important to remark that we were not able to assess the uncertainty that pertains the simulations from D-mesons, as currently the only updated hadronic interaction model for charmed hadrons is SYBILL-2.3 [53,[85][86][87]. The uncertainties are estimated following reference [88] (see also reference [89] for a similar discussion). For a given meson, we calculate the total expected flux dφ dE at the surface of the earth with different hadronic models. In order to do so, we let the mesons propagate through the atmosphere without decaying, which can be accomplished by switching off their decay in MCEq. The impact of the variation in the meson production rate for different models is then quantified with respect to a benchmark model, which we take here to be the SYBILL-2.3 model, by defining the ratio ∆ j (M ) between the meson flux predicted by the event generators that are being compared: where M is the meson of interest, j is the index used to denote the model that is being compared with the benchmark model BM, E min = 1.6 GeV is the minimum energy available in MCEq, and Λ = 10 3 GeV is the upper energy cutoff that we use in order to obtain the neutralino production rate from a given parent meson. Table I shows the ∆ coefficients obtained for the production of kaons with different event generators. The uncertainty as quantified by equation (7) reaches a maximum of around 66% for QGSJET. We calculate the neutralino production rate using different hadronic interaction models, and assess their impact on the region of parameter space that can be probed. We stress that dedicated efforts are needed in order to reduce the uncertainties involved in meson production in the forward direction, a problem that is also crucial for precise modeling of neutrino fluxes at the LHC [90][91][92].

III. NEUTRALINO SIGNALS AT SUPER-KAMIOKANDE
After light neutralinos are produced from the decay of mesons in the atmospheric shower, they decay to SM particles as they propagate. We assume that all the particles in the decay chain are highly boosted, and hence approximately collinear in their trajectories. The degree of attenuation of the flux that arrives at the detector depends on the relation between the lifetime of the neutralino and the distance it travels before reaching the detector. In particular, it is expected that the detector signal from up-going events (i.e. the neutralinos that arrive from below the detector), will be suppressed in comparison with down-going events (i.e. the neutralinos that reach the detector from above). Assuming that the production rate of neutralinos is symmetric with respect to the azimuthal component, the expected differential flux is obtained by integrating over the column depth X, according to where corresponds to the traveled distance of the neutralino, measured from the production point, which can be obtained from the height h and the zenith angle θ using the geometrical relation with R ⊕ the earth's radius. In equation (8), the lifetime of the neutralino enters via the exponential decay factor that depends on the decay length λχ0 The neutralinos that reach the Super-Kamiokande detector can decay to SM particles inside its instrumental volume, which we model as a cylinder of 20 meters in radius and 40 meters in height [45]. Among the possible decay products of the neutralinos it is possible to find charged leptons and mesons, as well as neutral pseudoscalar and vector mesons. Different decay products correspond to different kinds of signals in the detector, which at Super-Kamiokande can be classified in two: showering (or e-like), and non-showering (or µ-like events). The former originates from electromagnetic and hadronic showers in the detector, while the latter is primarily associated to muons, which leave a distinctive Cherenkov ring with crisp edges. In this work, we focus on the showering signals originated from the decay of neutralinos described in table II. In particular, for the case of neutral mesons in the final states, we assume that they decay to a showering signal inside the volume of the detector. The parameter choice of the benchmark B1 also allows K 0 L as a final state for neutralino decays, but we do not consider this as part of our signal since this particle will typically decay outside the detector due to its large decay length. For our benchmark scenarios, the important decay width formulas of neutralinos to pseudoscalar and vec-RPV coupling Production Decay modẽ tor mesons accompanied with a lepton are [32]. where In the equation above, M ij (M * ij ) represents one of the final state pseudoscalar (vector) mesons listed in the neutralino decays of table II. The indices i, j designate the family of the meson's valence quarks. For the final state pseudoscalar mesons, we use the decay constants f π 130 MeV [21], , where f K was defined after equation (4). For vector mesons, on the other hand, f V M * represents the vector meson decay constant, and for K * 0 , K * + we use the approximate value f T K * 230 MeV [32,93]. As mentioned before, although for the production of neutralinos from D-mesons and kaons we only need the trilinear couplings λ 121 and λ 112 , the decay of neutralinos produced from kaons into visible showers in SK proceeds via the coupling λ 111 . The two benchmark scenarios considered in this article can be seen in table II. For benchmark B1 we consider neutralino masses in the range m K + + m e ≤ mχ0 1 ≤ m D + − m e , while for B2 the neutralino mass lies in the range m π 0 + m e ≤ mχ0 1 ≤ m K + − m e . The election of B1 allows for a direct comparison with references [30,32], where the sensitivity reach for future experiments aiming to explore the lifetime frontier was estimated. The election of B2 allows us to study a novel production channel where neutralinos can be abundantly produced in atmospheric showers. Furthermore, it can be seen from equations (4) and (12) that it is possible to combine the dependence of the decay widths on the RPV parameters and the sfermion mass in the ratio λ ijk /m 2 f , which we set as a free parameter.
To obtain the expected number of events with an energy in the range Eχ0 1 and Eχ0 1 + dEχ0 1 , and with trajectories within cos θ and cos θ + d cos θ, we include the effective surface S eff for a decay to take place inside the SK detector, such that for a given time window ∆T the event rate is where the effective surface is obtained by integrating the surface of the detector that is perpendicular to the incoming direction of the neutralino flux, weighted by the probability that the particle decays inside the detector: Here, ∆ det is the segment of the particle trajectory that traverses the detector, for which explicit analytical expressions can be found in the appendix of reference [49]. The computation of the effective surface for decay is a purely geometrical problem. Figure 3 shows the effective surface as a function of the decay length of the neutralino, with trajectories fixed by different values of the cosine of the zenith angle given by 0.9, 0.5 and 0.1. The event distribution in equation (13) can be integrated in energies and trajectories within a given bin determined by the resolution of the experiment. As mentioned before, we use atmospheric neutrino data reported by the Super-Kamiokande experiment in reference [94]. This data contains the angular distribution of events involving electron and muon neutrinos, from different energy regimes; the Sub-GeV and Multi-GeV sample of events with energies below and above 1330 MeV, respectively. Taking into account the trilinear couplings considered in this article, as well as the minimum energy available in MCEq, we chose the showering (or e-like) event sample in the Multi-GeV energy window. The total events reported in reference [94] correspond to the SK-I up to the SK-IV data taking periods with a total run of 5,326 days. The number of events contained in the Multi-GeV range for the i-th bin in the cosine of the zenith angle is com- where is the detection efficiency, which for the Multi-GeV sample is flat in energy with a value of 0.75, while E min and E max are equal to 1.5 and 90.5 GeV, respec-  [32]. On the right panel, the green contour indicates the excluded parameter space when λ 121 = 0 and λ 112 = 0. Solid lines indicate constraints from neutral kaon oscillations [95] in green and DY processes at the LHC (labeled as 'Colliders') [96] in black, both evaluated at a sfermion mass of 1 TeV (for details, see text).
tively [94]. The background for our search corresponds to the expected e-like events from electron neutrinos at SK as reported in reference [94]. As an illustration, in figure 4 we show an example of the expected event signal generated from neutralino decays, when the couplings λ 112 and λ 111 are both fixed to a value of 0.002, and the neutralino mass is 0.31 GeV.

IV. RESULTS AND DISCUSSION
The signal computed in equation (15) depends on the neutralino mass and RPV couplings through its lifetime, branching fraction of production from mesons, and the fraction of neutralinos that decay to an e-like signal in the detector. We adopt the distribution for Poisson events, which implies a chi-squared statistic of the form where the sum runs over the 10 bins in cosine of the zenith angle, and B ci , D ci are the background and data in the i-th bin, respectively. We apply this statistical test to the SK data and derive constraints at 90% confidence level. The results for our two benchmark scenarios are displayed in figures 5 and 6.
There are numerous phenomenological constraints on the lepton number violating operator λ ijk L i Q j D c k (see [93] and references therein). In particular, as pointed out in [98], RPV parameters can be subject to strong constraints due to their contribution to meson oscillation observables. This is the strongest bound that applies to the benchmark B1 specified in table II, where the treelevel contributions to kaon oscillations induced by that choice of parameters imply the sneutrino mass dependent bound [95] |λ 112 λ 121 | ≤ 2.2 × 10 −8 mν e 1 TeV where mν e is the sneutrino mass. Our results for B1 are shown figure 5, where we evaluated the limits for a sneutrino and squark mass of 1 TeV. As it can be seen on the left panel, the limits imposed by kaon oscillations for degenerated sneutrino and squark masses exclude all the parameter space within the projected sensitivity reach at SHiP [32], our limits from SK, and the limit from the Drell-Yan (DY) monolepton process pp → ν at the LHC (labeled as 'Colliders'). The latter is a single coupling bound for the parameter λ 112 given by [96] λ 112 ≤ 0.16 where ms R is the strange squark mass, which is also set to 1 TeV to compare with our results. However, since  [95] in green and neutrinoless double beta decay (0νββ) [97] in black. The bounds are evaluated for degenerated sfermion masses set to a value of 1 TeV.
in Super-Kamiokande all the neutralino decay products listed in table II are visible as e-like events, we can still have a signal when the parameter λ 112 is set to zero and we are left with neutral kaons in the final states. In such a case, the limits from kaon oscillations do not apply (see figure 5, right panel), and we are left with the collider constraint from the DY dilepton process pp → + − for the parameter λ 121 , which is given by [96].
where mq is the squark mass. In this case, we find that the limit obtained from Super-Kamiokande data is better than the current constraint for the corresponding mass range.
On the other hand, in the case of benchmark B2, there is also a stringent limit from kaon oscillations on the product of the parameters, namely [95] |λ 111 λ 112 | ≤ 1.
Nevertheless, the most stringent constraint comes from its contribution to 0νββ. The current limit for the half life of this process as determined by the KamLAND-Zen experiment [99], imposes the upper limit [97] λ 111 ≤ 2.2 × 10 −3 mq 1 TeV where mq and mg are the masses of the squarks and gluinos, respectively. To contrast these bounds with our results, again we set both mass parameters to the value of 1 TeV (see figure 6). Remarkably enough, due to the higher kaon flux with respect to D-mesons, the limits from Super-Kamiokande in this parameter space turn out to be more stringent than the bounds that come from both kaon oscillations and 0νββ. The excluded parameter space for this benchmark scenario changes marginally with different choices of hadronic interaction models, as it can be seen from the lines that indicate the exclusion region obtained with each model (see the left panel of figure 6). The difference between the kaon flux obtained with EPOS-LHC and with our benchmark model SYBILL-2.3 does not translate to any visible impact in the excluded parameter space, therefore we omit its contour in the figure. As a caveat, we note that the bound shown in equation (21) relies upon the gluino dominance assumption, where the neutralino mass is a fraction of the order of 10 −2 times the gluino mass [59]. Finally, we emphasize that our results do not depend on the value of the sfermion masses. However, the bounds shown in equations (17)(18)(19)(20)(21) become more strict as sfermion masses increase and less strict when they are lowered. Moreover, in the case where the squarks and sneutrino masses are kept fixed and the gluino mass is increased, the bound from 0νββ becomes less stringent. Limits at 90% confidence level in the plane defined by the total branching ratio and the proper decay length of the long-lived neutralino. On the left panel, the purple contour shows the limits for SK derived in this article for benchmark B1, and the coloured lines displays the projected sensitivity reach at future experiments for the same benchmark scenario derived in references [30,32]. On the right panel, the limit obtained using Super-Kamiokande data for the benchmark scenario B2 is shown. There are no other studies for this case in the literature.
The excluded combinations of values for the neutralino masses and RPV couplings allow us to determine the region of the parameter space that can be probed in the plane defined by the total branching fraction, which is the production branching ratio of the neutralino times the branching ratio of its decay to a visible signal, and the lifetime of the neutralino. The results can be seen in figure 7. Given that for benchmark B1 we are considering the same production and decay channels as in references [32] and [30], we can compare the current region that can be excluded by Super-Kamiokande, with the sensitivity reach expected for FASER, and other possible future experiments including CODEX-b, MATHUSLA and SHiP (left panel of figure 7).
In the case where the neutralinos are produced in kaon decays, there are no studies on the expected sensitivity reach at the next generation of detectors. The results obtained in this case demonstrate the great capacity of the Super-Kamiokande neutrino detector to place stringent limits for long-lived particles, owed to its ability to probe particles with lifetimes with a peak sensitivity around 10 km and total branching fractions of order 10 −9 .

V. CONCLUSIONS
The lifetime frontier has emerged as a powerful line of exploration to search for beyond the Standard Model physics, specially in the absence of new signals at the LHC. R-parity-violating supersymmetry with light, longlived neutralinos constitutes a well-motivated scenario for new physics, with a rich phenomenology that has been gathering attention in recent years. In this work, we demonstrate a new way to search for long-lived neutralinos with masses of order 0.1 − 1 GeV that could be produced from the decay of charged D-mesons and kaons in cosmic-ray air showers, and whose visible decay can take place within large neutrino detectors such as Super-Kamiokande. We have analyzed two benchmark scenarios that include the couplings λ 121 , λ 112 , and λ 111 . In both cases, it is possible to improve the excluded region in parameter space when compared to existing constraints from colliders and neutrinoless double beta decay. Note, however, that this comparison requires to fix the mass of the sfermions, which we set to 1 TeV. An interesting feature about benchmark B1 is that in this case it is possible to compare the lifetime range that can be probed using SK data, against the expected sensitivity reach of nextgeneration, long-lived particle detectors. We find that the sensitivity for Super-Kamiokande peaks for lifetimes of the order of 1.0 km. In the case of the benchmark scenario B2, the advantage is twofold. On the one hand, there are currently no searches for light neutralinos produced from kaon decays. On the other hand, since kaons have a considerably higher production rate in air showers (as it can be inferred from figure 2), we find that it is possible to probe neutralinos with lifetimes of the order of hundreds of kilometers, while also being able to achieve a limit in λ 111 /m 2 f and λ 112 /m 2 f of order 10 −9 GeV −2 , assuming these RPV couplings are equal.
When possible, we quantified the uncertainty that results from the choice of different hadronic interaction models. This was not possible for the benchmark scenario B1, since SYBILL-2.3 is the only event generator that provides a state-of-the-art simulation of charmed mesons. For the case of the benchmark scenario B2, we find no strong dependence on the hadronic interaction model chosen to simulate the production of kaons in the atmospheric shower.
Overall, the results found in this study demonstrate the potential of large neutrino detectors to place limits on beyond the SM scenarios predicting long-lived particles, specially when contrasted with supersymmetric searches at colliders, where signals must be carefully chosen, and the reinterpretation of results is usually a complicated task. Finally, we stress that further scenarios can be pursued systematically along the direction presented in this work.