R-parity Violation and Light Neutralinos at CODEX-b, FASER, and MATHUSLA

The $LQ\bar D$ operator in R-parity-violating supersymmetry can lead to meson decays to light neutralinos and neutralino decays to lighter mesons, with a long lifetime. Since the high-luminosity LHC is expected to accumulate as much as 3/ab of data, several detectors proposed to be built at the LHC may probe unexplored regions in the parameter space, for long-lived neutralinos. We estimate the sensitivity of the recently proposed detectors, CODEX-b, FASER, and MATHUSLA, for detecting such light neutralinos singly produced from $D$- and $B$-meson decays in a list of benchmark scenarios, and discuss the advantages and disadvantages of the proposed detectors in this context. We also present our results in a model independent fashion, which can be applied to any long-lived particle with mass in the GeV regime.


I. INTRODUCTION
The discovery of a Standard-Model (SM) like Higgs boson in 2012 has been a highlight of the Large Hadron Collider (LHC) [1,2]. The small Higgs mass, m h = 125.09 GeV [3], has since, however, consolidated the hierarchy problem [4,5]. Supersymmetric (SUSY) theories offer an elegant solution, for reviews see Ref. [6,7]. All searches for the new fields predicted by SUSY however, have been unsuccessful yet. This leads to lower limits on the masses of squarks and gluinos in various supersymmetric models of order 1 TeV and above [8][9][10][11][12]. On the other hand, the lightest neutralino is not similarly constrained. In fact, if we drop the assumption for the gaugino masses, M 1 = 5 3 tan 2 θ W M 2 , which is motivated by grand unified theories (GUTs), and do not require the lightest neutralino to comprise the dark matter of the Universe, light neutralino masses below ∼ 10 GeV are still allowed [13][14][15][16][17][18]. For neutralino massess between about 1 eV and 10 GeV the relic energy density of the neutralinos would overclose the Universe [11], thus such neutralinos must decay. Since they are light they will typically have long lifetimes.
In R-parity-violating (RPV) SUSY (for reviews see Ref. [19][20][21]), the lightest neutralino is no longer stable and decays via RPV couplings. 1 Thus such a neutralino can be light. For small couplings and small mass these neutralinos can be long-lived enough to escape the reach of the LHC detectors. Moreover, the RPV couplings can induce single production of neutralinos via rare meson decays. Such scenarios have been investigated in various fixed-target set-ups [13,[26][27][28]. More recently they have also been studied in the context of the proposed SHiP experiment [29][30][31]. Ref. [31] studied the expected LHC sensitivity to such scenarios assuming an integrated luminosity of 250/fb. Ignoring differences in the reconstruction efficiency the sensitivity in the R-parity violating couplings at ATLAS was lower than at SHiP by roughly a factor of 2.
It is expected that the high-luminosity LHC (HL-LHC) will deliver up to 3/ab of luminosity in the coming 20 years [32]. As the cross sections for producing long-lived particles (LLPs) are typically small, such a large amount of data is required for high sensitivity to LLPs. Unsurprisingly, there have appeared several proposals to build new detectors near the interaction points (IPs) at the LHC, exploiting the projected large luminosity: CODEX-b [33], FASER [34] and MATHUSLA [35]. In this study, we estimate the sensitivity reach of these detectors for discovering singly produced light neutralinos from Dand B-mesons via RPV LQD couplings, and compare them with each other. We also interpret our studies in a model 1 Incidentally in such RPV SUSY models the dark matter can be composed of axinos [22][23][24][25]. arXiv:1810.03617v2 [hep-ph] 14 Mar 2019 independent way, independently of the RPV couplings. Instead we set bounds on the product of the branching ratios of the production of an LLP from a meson decay and the decay of the LLP to a meson and charged lepton in terms of the neutralino decay length cτ . This can be applied to any potential LLP.
CODEX-b is a comparatively small cubic detector making use of a shielded space near the LHCb IP that is expected to be free soon. Since it is to be installed at LHCb instead of ATLAS or CMS, CODEX-b will have an expected luminosity of 300/fb, one order of magnitude smaller than ATLAS or CMS, if LHCb runs until 2035 with upgrades to a Phase-II [36]. FASER was proposed as a small cylindrical detector to be built in the very forward region several hundred meters downstream of the ATLAS or CMS IP. In comparison, MATHUSLA would be built as a massive surface detector above the ATLAS IP. The details of the experimental setups are summarized in Sec. III.
The CODEX-b physics proposal [33] examined two benchmark models, i.e. Higgs decay to dark photons, and B-meson decays via a Higgs mixing portal. Several FASER papers have respectively studied dark photons produced through light meson decays and photon bremsstrahlung [34], dark Higgs bosons produced through Band K-mesons [37], heavy neutral leptons [38] and axion-like particles [39]. There are also studies investigating MATHUSLA with dark Higgs [40], exotic Higgs decays to LLPs [35,41], and the Dynamical Dark Matter framework [42]. Recently a MATHUSLA white paper [43] appeared, where the theory community presented detailed studies of MATHUSLA's potential of detecting LLPs in many different models. Ref. [44] studied all these three detectors with heavy neutral leptons in the Type-I Seesaw model, and the lightest neutralino pair-produced from Z bosons with the RPV-SUSY model. Very recently Ref. [45] investigated inelastic dark matter models at various existing and proposed LHC experiments including CODEX-b, FASER and MATHUSLA. We extend this work to consider the production of supersymmetric neutralinos via both Dand B-mesons, as well as the decays of the neutralinos to a charged meson and a charged lepton. RPV-SUSY is a complete model and we thus also consider the full kinematic constraints due to phase space. The mass differences between the mesons, the neutralino and a potential tau-lepton strongly affect the search sensitivities.
It is the purpose of this paper to investigate the discovery potential of light neutralinos at the detectors CODEX-b, FASER, and MATHUSLA. The primary motivation of this scenario is that supersymmetry is a potential solution to the hierarchy problem. Light neutralinos are consistent with all laboratory and astrophysical data [16,44,46,47]. Thus this is an allowed supersymmetric parameter range, and should be investigated. Such a light neutralino is only consistent with the observed dark matter density if it decays on time-scales much shorter than the age of the universe. This is the case for R-parity violating scenarios. R-parity violating supersymmetry naturally obtains light neutrino masses, without introducing a super heavy see-saw Majorana mass of order 10 10 GeV or higher [48,49]. The scenario of an O(1 GeV) neutralino does not in itself resolve any discrepancy between the Standard Model and current data.
This paper is organized as follows. In Sec. II we briefly introduce the model of RPV-SUSY. In Sec. III we describe the three experiments for which we estimate the sensitivities, and explain the details of the numerical simulation. In Sec. IV we present results for various benchmark choices for RPV couplings. We summarize and conclude in Sec. V.

II. SUPERSYMMETRY WITH RPV, PRODUCTION AND DECAY OF LIGHT NEUTRALINOS
We give a brief introduction to the RPV-SUSY model, and describe the production and decay of light neutralinos via RPV couplings. Compared to the R-parity conserving (RPC) supersymmetric theories, RPV-SUSY has extra terms in the superpotential: 2 where the first three terms are lepton number violating (LNV) and the last is baryon number violating (BNV). The co-existence of LNV and BNV terms would lead to too fast proton decays, so in our study we choose to be exclusively interested in the LQD operators. With nonvanishing RPV couplings, the lightest supersymmetric particle (LSP) is not stable and can decay to SM particles. If the lightest neutralino is sufficiently light, it can be the LSP. We assume this is the case in our study. Neutralinos that are produced from charm and bottom meson decays are necessarily lighter than 10 GeV and are dominantly bino-like to avoid existing bounds, see Ref. [16]. Formulas for the partial widths of heavy meson decays and for the partial widths of neutralino decays via LQD couplings can be found in Refs. [13,16,26,31].
In principle, one single LQD coupling introduces several effective SM operators and hence may simultaneously induce both meson decays to neutralinos and neutralino decays to lighter mesons. However, as Ref. [31] points out, due to kinematic constraints only the coupling λ 112 may lead to such a complete decay chain: Moreover, since the mass difference between K 0 L/S and K ± is only 4 MeV, the kinematically allowed neutralino mass range is very small and this case is not worth studying. Therefore, we only consider scenarios with two distinct non-vanishing RPV operators, one for the production of the neutralinos and the other for the decay.
The couplings λ ijk for the operator L i Q jDk have strict bounds from different sources, though the bounds can be substantially weakened for heavy sfermion masses above 1 TeV. For reviews, see Refs. [19,[53][54][55][56]. Since we investigate the same benchmark scenarios as in Ref. [31], we only list the relevant bounds, reproduced from Ref. [56]: Some pairs of operators have even stricter product bounds. We take the relevant bounds from Ref. [53]: Throughout this work, we assume that all sfermions have degenerate masses mf . This allows us to directly compare the above bounds to our results even though the respectively relevant operators depend on the masses of possibly different SUSY particles. Note that results for significantly non-degenerate SUSY spectra may therefore differ significantly and can change the relative importance of bounds from different sources.

III. EXPERIMENTAL SETUPS AND SIMULATION
In this section we summarize the setups of the detectors, and explain in detail our simulation procedure. For more information on the proposed detectors we refer to Refs. [33][34][35]. The main difference between them lies in the projected luminosity, the geometry and installed position (large or small pseudorapidity η). Moreover, whether installed underground or above the ground, these proposals all argue that the background influence, for example from cosmic rays, can be well controlled. Therefore, we do not discuss it and always assume 100% detector efficiency.
In order to estimate the number of LLP decays inside the respective detector's chamber, we take into consideration both the production of the mesons and hence the neutralinos, and the decay of the neutralinos via mesons. On the production side, since we study neutralinos produced from meson decays, we use results published by the LHCb collaboration [57,58] for estimating N M , the total number of the meson type "M" produced at the LHC. For D-mesons, we consider only neutralinos produced from D + -and D s -mesons which are relevant for the benchmark scenarios we consider. Ref. [57] gives the cross section for producing D + , and D * + . The latter decays to D + -, and D s -mesons at the 13-TeV LHC for a certain kinematic range: 0 < p T < 15 GeV/c and 2.0 < y < 4.5, where p T denotes transverse momentum and y rapidity. We use the computer program FONLL [59][60][61][62] to extrapolate these numbers to the whole kinematic range, and, after taking into account the decay of the D * + -to D + -mesons, we obtain the total numbers of D + and D s produced over the hemisphere for L = 3/ab: N Ds = 5.11 × 10 15 .
At the LHC the mesons are produced over the full 4π. However the detectors we are considering here for LLPs are always off to one side of the collision point. We thus at first only consider the forward or backward hemisphere (2π) which contains the respective detector. We then later impose the necessary geometric cuts corresponding to a specific detector. Among the B-mesons, B 0 and B ± are of interest here. Ref. [58] presents the experimentally measured b-quark production cross section at the 13-TeV LHC for 2 < η < 5, and the corresponding number after extrapolation over the full η range with the numerical tool Pythia 8 [63,64]. We take the fragmentation factors of B-mesons directly from Ref. [31], which were obtained by simulating 1 M events of HardQCD:hardbbbar in Pythia 8 [63,64]. We obtain over a hemisphere for L = 3/ab. The branching ratios of these mesons decaying to neutralinos are easily calculated with the formulas given in Ref. [31]. We arrive at the following expression for the total number of neutralinos produced in a hemisphere N prod where l is the associated lepton in the meson decay, which can be charged or neutral. We then apply Monte Carlo (MC) techniques to determine the average probability of the neutralinos decaying inside the detector chamber, where P [(χ 0 1 ) i in d.r.] is the probability for a given generated neutralino to decay in the "detectable region", "d.r.". Dividing by the total number of simulated neutralinos produced, N M C χ 0 1 , gives the average. We explain how to calculate P [χ 0 1 in d.r.] for each detector in detail below.
Since it is difficult to experimentally reconstruct the trajectory of the neutral final-state particles of the neutralino decays, we consider only charged decay products to be detectable. (See Ref. [31] for a discussion of the potential influence of decays to K 0 's.) The final number of observed neutralino decays is expressed as N obs We use Pythia 8.205 [63,64] to perform the MC simulation in order to calculate P [χ 0 1 in d.r.] in Eq. (14). We use two matrix element calculators of Pythia, namely HardQCD:hardccbar and HardQCD:hardbbbar, to generate initial Dand B-mesons, respectively. Note that the differential cross section of producing heavy flavor mesons in the very forward direction, where FASER sits, is not validated in Pythia. In order to solve this problem, we reweigh the Pythia meson production cross section at different ranges of transverse momentum and pseudorapidity by the corresponding more reliable numbers calculated by using FONLL. We simulate 20,000 events for each benchmark scenario and extract the kinematical information of each neutralino (χ 0 1 ) i from Pythia: Here the z-direction is along the beam pipe, p z i is the component of the 3-momentum along the z-axis, E i is the total energy of the neutralino, and θ i , φ i are the polar and azimuthal angles, respectively. With this kinematical information we derive the relativistic quantities as follows: where Γ tot (χ 0 1 ) is the total decay width ofχ 0 1 and can be calculated with formulas given in Ref. [31], λ i is the decay length of (χ 0 1 ) i along the direction of its movement in the lab frame and λ z i is the z-component of λ i . These quantities are used to calculate P [(χ 0 1 ) i in d.r.] for each detector. We now discuss the detectors in turn.

A. CODEX-b
CODEX-b ("Compact detector for Exotics at LHCb") [33] was proposed as a cubic detector with dimension 10 3 m 3 , sitting inside an underground cavity at a distance L = 25 m from the LHCb IP. The differential production distribution is flat in the azimuthal angle and the azimuthal coverage of the detector is about 0.4/2π ≈ 6%. The polar angle range of the CODEX-b experiment at the appropriate azimuthal angle is between 11.4 • and 32.5 • . This corresponds to the pseudo-rapidity range η ∈ [0.2, 0.6]. For this narrow range, and at the precision of this analysis, we also treat the polar angle differential production distribution as flat. As we mentioned earlier, LHCb is expected to have a total integrated luminosity of 300/fb, smaller by one order of magnitude than ATLAS or CMS. We calculate P [(χ 0 1 ) i in d.r.] with the following expression: where we approximately treat the box detector as a spherical shell segment with the volume length L d = 10m. η i is the pseudorapidity of (χ 0 Fig. 1. FASER ("ForwArd Search ExpeRiment") [34] proposes to build a small cylindrical detector placed a few hundred meters downstream of the ATLAS or CMS IP in the very forward region. In a series of papers [34,[37][38][39] several different variants of FASER have been proposed. In this paper, we focus on a recent setup, which would sit at a particularly promising location in the side tunnel TI18 [39]. We denote the distance from the IP to the near end of the detector as L = 470m, the radius of FASER as R = 1m, and the detector length as L d = 10m. Following is the expression for calculating the probability for a given neutralino to decay inside FASER: There is no azimuthal angle suppression because the FASER detector is cylindrical. Here the three cases cor-respond respectively to 1) the extended potential neutralino trajectory misses the decay chamber, 2) the extended potential neutralino trajectory passes through the entire length of the detector, and 3) the extended neutralino trajectory exits through the side of the detector. In practice, we treat the third case as negligible. It corresponds to the very narrow angular range And furthermore the decay products of the neutralinos may exit through the side and may thus miss the detector. These neutralinos hence would not be detected. A sketch of the geometric configuration of FASER is shown in Fig. 2. In Ref. [35] it has been proposed to construct a surface detector 100 m above the ATLAS IP called MATHUSLA ("MAssive Timing Hodoscope for Ultra Stable neutraL pArticles"). The detector should be horizontally offset by 100 m from the ATLAS IP and with a massive dimension of 200m×200m×20m, MATHUSLA is expected to have excellent sensitivity for detecting LLPs. Below we show the formulae for calculating P [(χ 0 1 ) i in d.r.] in MATHUSLA: Here, L h and L v are the horizontal and vertical distance from the IP to the near end of MATHUSLA, and they both equal 100m. L d = 200m is the horizontal length of MATHUSLA and H = 20m is its height. The factor 1/4 comes from the azimuthal angle coverage. Both MATHUSLA and FASER expect to have 3/ab luminosity of data by ∼ 2035. We show the schematic plot of MATHUSLA in Fig. 3.

IV. RESULTS
We present our numerical results in this section. In Ref. [31] a series of benchmark scenarios representative of LQD couplings were investigated. In these scenarios, both the light lepton flavor (electron/muon) and the heavy tau flavor are considered, as the τ lepton leads to large phase space suppression effects. Also, different neutral or charged Dand B-mesons which would decay to the neutralino are considered; this is important because the cross sections of producing these mesons substantially differ, cf. Eqs. (9)- (12). In the present study, as a follow-up work to Ref. [31], we choose to focus only on one key benchmark scenario which features the important characteristics for a comparison of the proposed LHC(b) detectors' sensitivities, while only briefly discussing the other scenarios. We first consider the explicit RPV model and then also discuss the model-independent case.
Since the operators for production and decay scale with λ /m 2 f , we have three free parameters in the theory, after assuming that all SUSY fermionsf have degenerate masses, cf. Sec. II, namely: λ P /m 2 f , λ D /m 2 f , and mχ0 1 . Here λ P/D is the LQD coupling giving rise to the production/decay of theχ 0 1 , and m 2 f is the sfermion mass relevant for the production/decay process, respectively. 3 We therefore show model-dependent plots in two separate planes for the aforementioned benchmark scenarios: mχ0 1 vs. (λ P /m 2 f = λ D /m 2 f ) and λ P /m 2 f vs. λ D /m 2 f . For the latter plane, we present results for three different values of mχ0 1 . In addition, we present model-independent results in the plane BR vs. cτ for a generic LLP. Here cτ is the decay length of the LLP, and BR is the product of the branching ratios of the respective meson decaying to the LLP and of the LLP decaying to a charged meson and a charged lepton. These results can be interpreted in terms of any LLP which has the same or similar reaction chain.
For the key benchmark scenario, we choose to show all three types of plots while for the others we select only one single type of plot, where the distinctive features of the scenario may be best emphasized. Depending on the 3 The explicit formulae including the dependence on the relevant sfermion masses are given in Ref. [31]. λ P for production λ 122 λ D for decay λ 112 produced meson(s) Ds visible final state(s) K ± e ∓ , K * ± e ∓ invisible final state(s) via λ P (η, η , φ) + (νe,νe) invisible final state(s) via λ D (K 0 L , K 0 S , K * ) + (νe,νe) exact construction of the detectors, they can possibly also track neutral mesons. We thus show sensitivity estimates for two cases: 1) only charged final states can be tracked, and 2) both neutral and charged ones.
A. Benchmark Scenario 1 We begin with the RPV scenario we consider in detail in this study with D s -mesons produced at the LHC, which decay to a neutralino, which in turn travels for a macroscopic distance before decaying to a kaon and a lepton. In this scenario we assume λ 122 and λ 112 are the only non-vanishing LQD couplings. λ 122 gives rise to the production ofχ 0 1 via (27) and to the invisible neutralino decaỹ On the other hand, λ 112 leads to both visible and invisible decaysχ (decay via λ 112 ) (29) The invisible decays are important to take into account in the evaluation because they affect the total width of χ 0 1 . We summarize this scenario in Tab. I. We now present our results. In Fig. 4 we show modeldependent sensitivity estimates for the three detectors: CODEX-b, FASER and MATHUSLA. In the left column, plots are presented in the plane mχ0 1 vs. (λ P /m 2 f = λ D /m 2 f ). We have set λ D = λ P , and vary their values and the mass ofχ 0 1 . In order to see how the number of neutralino decay events change with varying mass and RPV couplings, we show the light blue, blue, and dark blue areas corresponding to the parameter space where respectively ≥ 3,  We do not show the product bound from Eq. (7), as for a 1 TeV sfermion mass it coincides almost exactly with the 5 TeV bound on the single couplings. The bound on λ /m 2 scales linearly with the sfermion mass, when taking the scaling of the bound on λ into account.
The 3-event dashed contour isocurve is extended to the lighter shaded region, bounded by a dotted line; this is obtained when we assume that invisible decays of the neutralinos can be detected as well. Whether this will be possible is an outstanding experimental question. In any case, we observe that for this benchmark scenario this would only give a very small extension in the sensitivity reach.
The range of sensitivity in the neutralino mass mχ0 1 is strictly determined by the kinematics of the production and decay and is thus identical for the three experiments. The range in sensitivity in λ /m 2 is determined by the experimental set-up. Comparing the results for the three detectors, we find that for this model CODEX-b and FASER reach similar values of λ /m 2 f , while MATHUSLA is more sensitive by a factor ∼ 5. Furthermore they can all extend well beyond existing low-energy limits on the R-parity violating couplings.
On the right in Fig. 4, we show plots in the plane λ P /m 2 f vs. λ D /m 2 f for three values of mχ0 1 : 600 MeV (light blue region), 1200 MeV (blue region), 1800 MeV (dark blue region). In this benchmark scenario, λ P = λ 122 and λ D = λ 112 . For these results, the requirement that λ P = λ D is lifted, so we observe an interplay between the production and decay ofχ 0 1 . We may compare each detector's sensitivity range in different parameters. For example, the λ P /m 2 f reach of FASER is only weaker than that of MATHUSLA by a factor ∼ 3, even though FASER is more than 25,000 times smaller than MATHUSLA. This arises because FASER exploits very well the advantage of receiving the light D-mesons (and the produced neutralinos) boosted in the very forward direction, where the differential production cross section is significantly higher. As for the reach in λ D /m 2 f , MATHUSLA shows again the strongest potential. Here we include single coupling bounds as solid hashed lines for three different sfermion masses (250, 1000 and 5000 GeV) and now also the product bound as a dashed hashed line for a 1 TeV sfermion mass. Again all experiments are sensitive well beyond existing limits.
We next consider a model-independent description, where we interpret our results in terms of the physical observables, BR= BR P · BR D , instead of the RPV-SUSY parameters. Here and we allow for any LLP. The results are shown in where βγ is the average boost of the neutralinos flying in the direction of the detector and L is the distance from the IP to the middle of the respective detector. We estimate βγ of the neutralinos that fly inside each detector by simulating 10,000 events in each case, which agrees with Fig. 5. The BR position of the valleys is determined by the luminosity of the experiment, the cross section of producing D s -mesons, the pseudorapidity coverage, the volume of the detector and the product of the branching ratios. The BR reach of CODEX-b is roughly one order of magnitude larger than that of FASER. This is mainly due to the fact that LHCb has a one order of magnitude lower projected luminosity than that of ATLAS/CMS. Perhaps more importantly, in spite of the huge volume difference between MATHUSLA and CODEX-b/FASER, the BR reach in MATHUSLA is only one order of magnitude stronger than that in FASER. For large cτ values MATHUSLA performs far better than CODEX-b, but for shorter neutralino lifetimes the detectors perform equally well. The reason is that the distance traveled to MATHUSLA is about ten times larger than for CODEX-b, such that less neutralinos reach the former detector for short-lived neutralinos. This leads to a similar sensitivity despite the larger integrated luminosity and the larger detector size of MATHUSLA.
Note that Fig. 5 is very similar to the first plot of Fig. 1 in Ref. [44], the result of which was obtained in λ P for production λ 121 λ D for decay λ 112 produced meson(s) D ± visible final state(s) K ± e ∓ , K * ± e ∓ invisible final state(s) via λ P (K 0 L , K 0 S , K * ) + (νe,νe) invisible final state(s) via λ D (K 0 L , K 0 S , K * ) + (νe,νe) the context of a Type-I Seesaw model, where the righthanded neutrino is the LLP with a mass of 1 GeV produced from D-meson decays. This illustrates the modelindependence of the results shown in the BR-cτ -plane.

B. Benchmark Scenario 2
Now we briefly study the other the benchmark scenarios. In Benchmark Scenario 2, λ P = λ 121 instead of λ 122 , so that a D ± , instead of a D s , decays to the lightest neutralino. Correspondingly, the invisible final states due to λ P are now kaons, instead of η, η , φ. The relevant information is summarized in Tab. III.
The model-dependent results are very similar to those shown in Fig. 4. We do not show them again. One difference is that the low-energy product bound is stricter in this case, cf. Eq. (6). It is due to K 0 −K 0 -mixing and scales linearly with the sneutrino mass. As pointed out FIG. 6. Model-independent sensitivity estimate of CODEX-b, FASER and MATHUSLA for Benchmark Scenario 2. The format is the same as in Fig. 5. The plot corresponds to visible decay products only. The two LLP mass values are 600 and 1200 MeV. λ P = λ 121 , λ D = λ 112 . λ P for production λ 131 λ D for decay λ 112 produced meson(s) B 0 ,B 0 visible final state(s) K ± e ∓ , K * ± e ∓ invisible final state(s) via λ P none invisible final state(s) via λ D (K 0 L , K 0 S , K * ) + (νe,νe) in Ref. [31], in the case of SHiP, also here, if the sneutrino mass is equal to the relevant squark mass of production and decay, the sensitivity reach of CODEX-b and FASER is excluded by these low-energy bounds. If there is a strong hierarchy and the sneutrinos are (unexpectedly) significantly heavier than the relevant squarks, this scenario is still viable. All the same, we present the modelindependent results in Fig. 6 for the same LLP mass values as in Benchmark Scenario 1: 600 and 1200 MeV. We again drop the curve for the 1800 MeV neutralino mass. The main difference between Fig. 6 and Fig. 5 is the BR reach. Benchmark Scenario 1 has a weaker BR reach mainly because N D + 3 · N Ds , cf. Eqs. (9), (10), and the neutralinos have a smaller branching ratio to charged particles.

C. Benchmark Scenario 3
We now study several scenarios where bottom mesons decay to a neutralino. Since the bottom mesons are much heavier than the charm mesons, the mass reach improves compared to the previous scenarios. In the present Benchmark Scenario 3, as before, we have λ D = λ 112 giving both invisible and visible neutralino decays. For the neutralino production we have λ P = λ 131 such that B 0 (andB 0 ) decay to a neutralino. This is summarized in Tab. IV. Kinematically we can thus probe For this scenario we only show the model-dependent plots in the plane λ P /m 2 f vs. λ D /m 2 f in Fig. 7. Here the sensitivity reach in λ P/D /m 2 f is only slightly weaker than the previous D-meson scenarios, mainly because at the LHC, the cross section of producing B-mesons is only smaller than that of D-mesons by a factor of ∼ 10 − 50, cf. Eqs. (9)- (12). In the previous scenarios, CODEX-b shows similar sensitivity reach in λ P/D /m 2 f to that of FASER, but now the former exceeds the latter, despite the fact that its projected luminosity is smaller by one order of magnitude. This is because the B-meson mass λ P for production λ 131 λ D for decay λ 121 produced meson(s) B 0 ,B 0 visible final state(s) D ± e ∓ , D * ± e ∓ invisible final state(s) via λ P none invisible final state(s) via λ D (K 0 L , K 0 S , K * ) + (νe,νe) is more than twice the D-meson mass, and hence the produced B-mesons are not as much boosted in the very forward direction as the D-mesons. For the same reason, we also have a larger sensitive mass range than in the previous benchmark scenarios. MATHUSLA again has the most extensive sensitivity range.

D. Benchmark Scenario 4
In this scenario, we use the same λ P as in the previous scenario, but change λ D from λ 112 to λ 121 . Correspondingly the decay mode of the neutralino changes from the decay to a K ± to a D ± , though the invisible decay mode remains the same. We summarize the relevant information in Tab. V. The kinematic reach in the neutralino mass is for the charged decay modes. It is extended when the invisible modes are included by replacing M D ± → M K 0 . We present the results of this scenario in the plane (λ P /m 2 f = λ D /m 2 f ) vs. mχ0 1 in Fig. 8. The lower mass sensitivity is now raised up to the D-meson mass, if only the visible decays of the neutralinos are considered. If we consider the detectors able to track neutral final states, the lower mass sensitivity is dramatically extended, as expected, down to m K ∼ 500 MeV. For large values of (λ P /m 2 f = λ D /m 2 f ) we produce many more neutralinos, but they now mostly decay before reaching the detector. That is why there is no sensitivity here. For very small values of (λ P /m 2 f = λ D /m 2 f ) we produce too few neutralinos and the neutralinos decay well after the detector. Otherwise the results are similar to those in Benchmark Scenario 3. λ P for production λ 313 λ D for decay λ 312 produced meson(s) B 0 ,B 0 , B ± (+ τ ∓ ) visible final state(s) K ± τ ∓ , K * ± τ ∓ invisible final state(s) via λ P none invisible final state(s) via λ D (K 0 L , K 0 S , K * ) + (ν,ν) While the previous benchmark scenarios concern only the light electron, we here explore the effect of the heaviest lepton τ on the sensitivity estimates. We consider λ P = λ 313 and λ D = λ 312 . Then both B 0 (B 0 ) and B ± may decay to a neutralino. In particular, B ± decays then include a τ lepton along with a neutralino. This gives a large suppression in phase space and a correspondingly lower neutralino mass sensitivity. While λ P does not induce any invisible decay of the lightest neutralino, λ D = λ 312 leads to both visible decays to a kaon and a τ , and invisible decays to a kaon and a ν. We summarize the information in Tab. VI. The mass sensitivity is given by We present a plot for this benchmark scenario in the plane BR vs. cτ in Fig. 9. We restrict ourselves to the production via B 0 . As can be seen in Tab. II, between the B 0 and the B ± cases the βγ values, which determine the cτ sensitivity, differ between 15 and 20%. This is below the resolution of our logarithmic plot. For the BR sensitivity the dominant contribution is the B 0 vs. B ± production rate, however they are almost identical, cf. Eqs. (11), (12). The only real difference is that for B ± we must have mχ0 1 < ∼ 3500 MeV, i.e. the medium mass case (3750 MeV) is not possible.
In Fig. 9 the labeling is similar to the previous scenarios. For each detector, i.e. FASER (blue), CODEX-b (pink) and MATHUSLA (red), the light/medium colors correspond to the lightest(2500 MeV)/medium(3750 MeV) mχ0 1 . Fig. 9 is very similar to the right panel of Fig. 1 of Ref. [44] where the sensitivity to sterile neutrinos was   After having presented results of different benchmark scenarios in the previous subsections, we supplement our results by showing in Fig. 10 the decay branching ratios of theχ 0 1 to visible, i.e. charged meson final states, as a function of mχ0 1 in the kinematically allowed range for all the scenarios. The curves can be well understood by considering the kinematic thresholds for the various neutralino decays. For Benchmark Scenarios 1, 2 and 3 the first decay channel to open is to the charged kaon: χ 0 1 → K ± e ∓ . Thus the visible branching ratio starts at 1. It rapidly drops as the K 0 -threshold is crossed. The asymptotic value of the branching ratios in Benchmark Scenarios 1, 2 and 3 is simply determined by the number of charged or neutral decay channels. For example, in Benchmark Scenario 2 above all thresholds there are 4 visible final states and 8 invisible final states, giving a branching ratio to visible of 4/(4+8)=1/3. The bump in Benchmark Scenario 1 (red curve) is due to the extra threshold of the η-meson below the K + * and K 0 * masses. Note that it is the vector meson K 0 * which is relevant for the neutralino decay, not the pseudoscalar [31].
In Benchmark Scenarios 4 and 5 the first kinematically accessible final state is invisible: K 0 ν e . Thus for low values of the neutralino mass the visible branching ratio vanishes. In Benchmark Scenario 4 the first visible decay mode to open up is: D + e − , in Benchmark Scenario 5 it is: K + τ − . We also point out that we divide Scenario 5 into two cases: the neutral bottom meson decaying to a neutralino, and the charged bottom meson decaying to a neutralino accompanied by a tau lepton. This leads to the overlap of the two corresponding curves, the blue and the black ones, below a mass of ∼ 3.5 GeV. Above that mass value, only B 0 may decay to such a neutralino, because of the large mass of the tau lepton.

G. Summary of our Results
We finish this section with some conclusions drawn from the above results for specific benchmarks. For all scenarios, we observe similar reach in λ /m 2 f for the two experiments CODEX-b and FASER and the strongest sensitivity for MATHUSLA. Even though FASER takes good advantage of the boost of the Dand B-mesons in the very forward direction, MATHUSLA overcomes this disadvantage by virtue of its much larger volume. Compared to earlier results determined for the SHiP experiment [31], both FASER and CODEX-b have a smaller expected reach in λ /m 2 f . Even MATHUSLA cannot outperform SHiP in scenarios with D-meson dependence, because SHiP's centreof-mass energy of ≈ 27 GeV results in very high sensitivity. For models with B-meson decays, however, we expect MATHUSLA's sensitivity to be comparable or even better than SHiPs.
We also translated our results into sensitivity limits on the meson's branching ratio BR and here we observe differences in the experimental sensitivities of FASER and CODEX-b, depending on the meson flavour: while for the B-meson's BR we observe a similar ranking of the two experiments in the reach for λ /m 2 f , FASER's reach is expected to be slightly stronger than CODEX-b's in case of scenarios with D-meson decays. In all cases MATHUSLA shows the largest sensitivity reach.
As discussed above, our results in terms of cτ vs. BR can be used to estimate experimental sensitivities for long-lived particles different from RPV neutralinos. We used our results for Benchmark Scenario 1 as an example by pointing out the strong resemblance of Fig. 5 and the first plot of Fig. 1 in Ref. [44] determined for a different BSM model albeit with similar decay topology. Though we use this observation to label these limits as model independent, our combined set of results still points to several sources of model-dependence. An important degree of freedom to mention here is the LLP mass. For fixed cτ and branching ratio, changing the mass of the LLP has an important impact on its kinematic parameters βγ and η which, as we have shown, play a non-negligible role as they are strongly connected to the preferred distance L of the detector to the primary interaction vertex. Another important aspect neither covered in cτ nor the LLP BR is the overall topology which leads to the production of the LLP. Though all our benchmarks share the same topology pp → meson+X, meson → LLP+Y , we observe sizable differences in the experimental coverages depending on the flavor of the produced meson, here D or B. Not only do these have different total production cross sections but also their own kinematics -and with that the kinematics of the LLPs they decay into -differ. A full "model-independent" analysis would require the consideration of several additional degrees of freedom, some of which cannot be formulated as a continuous parameter like the overall production-and-decay-topology of the LLP. This results in an unfeasible, if not impossible, exercise. Nevertheless, although the dependence on these additional parameters may not be explicitly covered in our chosen degrees of freedom cτ and BR, our results can still be applied to a large class of LLP models different from RPV as long as they share similarities to the topologies discussed here.

V. CONCLUSIONS
We have investigated the sensitivity of three recently proposed detectors at the LHC: CODEX-b, FASER, and MATHUSLA with respect to the detection of light longlived neutralinos in RPV-SUSY scenarios. The neutralinos are produced and decay via the RPV LQD oper-ator with coupling λ . We studied five representative benchmark scenarios of the RPV couplings proposed in Ref. [31] where a similar sensitivity study for SHiP was completed. In general CODEX-b and FASER show similar reach in λ /m 2 f , where mf is the mass of supersymmetric fermion partners, while MATHUSLA performs better by approximately one order of magnitude. Comparing MATHUSLA results with SHiP estimates, we find that MATHUSLA shows a better sensitivity in scenarios involving B-meson while it provides only comparable or slightly weaker results than SHiP for models in which neutralinos interact with D-mesons.
We also want to point out that cosmic rays may provide an additional argument to choose an underground experiment like FASER over a surface experiment like MATHUSLA. Although we ignored this allegedly control-lable background contamination in our analysis, the required workload to fully control this background source may be significantly different between the experiments discussed here.