Searches for Atmospheric Long-Lived Particles

Long-lived particles are predicted in extensions of the Standard Model that involve relatively light but very weakly interacting sectors. In this paper we consider the possibility that some of these particles are produced in atmospheric cosmic ray showers, and their decay intercepted by neutrino detectors such as IceCube or Super-Kamiokande. We present the methodology and evaluate the sensitivity of these searches in various scenarios, including extensions with heavy neutral leptons in models of massive neutrinos, models with an extra $U(1)$ gauge symmetry, and a combination of both in a $U(1)_{B-L}$ model. Our results are shown as a function of the production rate and the lifetime of the corresponding long-lived particles.


Introduction
The Standard Model (SM) of particle physics has demonstrated to be an extremely predictive theory, which however cannot be complete. In particular, it lacks a mechanism to generate neutrino masses, dark matter, and the baryon asymmetry of the Universe that we observe today. While LHC upgrades and future colliders will keep pushing the energy frontier forward in order to search for new physics at high energies, it is important to keep in mind that the physics beyond the Standard Model (BSM) may very well involve weakly interacting particles at low scales instead, which in minimal extensions of the SM are expected to be relatively long-lived. In fact, the search for such long-lived particles (LLP) in colliders and beam-dump experiments has become a field of intense activity in recent years and several experiments have been proposed to conduct dedicated searches [1][2][3][4].
In this work, we study the potential of atmospheric neutrino detectors to search for these particles since they can also be copiously produced in atmospheric showers. In order to do this, we compute the expected flux of exotic LLPs being produced in proton-nucleus collisions in the atmosphere, either through proton bremsstrahlung or as a product of meson decays. After traveling typical distances of tens of kilometers through the atmosphere, such LLPs may subsequently decay within the detector volume. The signal would therefore be an excess of events where the interaction vertex is contained in the fiducial volume of the detector. These are a priori indistinguishable from an atmospheric neutrino interactions, which therefore constitute an irreducible background to this search.
The experimental setup considered in this work has two main advantages with respect to searches in laboratory experiments. First, the center-of-mass energy in proton-nucleus collisions in the atmosphere will typically be much higher than at beam-dump experiments, allowing to produce an abundant flux of heavy mesons and tau leptons. Secondly, the size of these huge detectors results in larger decay volumes and could be optimal in the searches for particles with long lifetimes, of about O(10) km in the laboratory frame. In this work, we will be considering two experimental setups: Super-Kamiokande, most sensitive to events with energies below 100 GeV; and Icecube, sensitive to events with energies at the TeV scale and above. Therefore, the results between the two will be largely complementary as they will be sensitive to different values of the lifetime of the particle in its rest frame, τ . As we will see, the best sensitivities will be achieved for τ /m ∼ O(1) km/GeV for Super-Kamiokande and τ /m ∼ O(10) m/GeV for Icecube.
In this work, for concreteness, we will focus on two minimal and theoretically wellmotivated extensions of the SM: (1) a scenario with heavy neutral leptons (HNL), that could be responsible for generating neutrino masses [5][6][7][8] and possibly the baryon asymmetry of the Universe [9,10], and (2) a model with an extra U (1) symmetry which, after being broken, leads to a massive dark photon and a portal to the dark sector [11,12]. In the first case we focus on LLPs in the mass range between the Kaon and D-meson masses, where existing laboratory bounds are weaker [13][14][15][16], while in the latter we will consider dark photons below the proton mass. The search for lighter, MeV range, atmospheric sterile neutrinos has been considered before in Super-Kamiokande [17,18] and IceCube [19]. Our analysis significantly improves the methodology of these earlier studies.
In the most minimal models with just HNL or dark photons, both production and decay of the LLP are highly correlated and, thus, longer lifetimes also mean lower production rates. However, in non-minimal extensions this may no longer be the case, as different mechanisms may be involved in production and decay processes. An example of this framework is provided by a B − L model with HNL [20]: in this case, production of the HNL could take place through the B − L gauge interaction [21,22], while the decay may occur via its mixing with the SM neutrinos. Such scenario will also be discussed in our work. Our results will be presented from a model-independent perspective, considering a more general BSM scenario where the production and decay rates could be decoupled, and we will discuss the interpretation of our results in the three scenarios outlined above.
The paper is organized as follows. In Section 2 we discuss the production of LLPs and describe the method to compute LLP fluxes from two main sources: meson decays (two-and-three body) and proton bremsstrahlung. In Section 3 we present the general procedure to compute LLP signals in a neutrino detector which requires defining detector effective areas for decaying events, detector efficiencies, backgrounds and the data samples that will be used to extract limits on LLPs. In Section 4 we focus on three BSM scenarios and use the results of the previous sections to compute the required fluxes. In sections 5 we present the sensitivities to LLPs of searches in IceCube and Super-Kamiokande, and our conclusions are drawn in section 6.

LLP production in atmospheric showers
The flux of any SM particle produced in the atmosphere is usually given as a function of the slant depth, X, related to the integral of the density in the direction of the particle from the top of the atmosphere to the production point at distance [23]: where ρ(h) is the atmospheric density at height h, and θ the zenith angle defined by the trajectory of the particle. The distance is related to the height h and to the zenith angle θ as: where R ⊕ 6370km is the Earth radius. In our calculations, we consider h( max , θ) = 80 km as the maximum vertical height of the atmosphere.
In our calculations, we have used the Matrix Cascade Equation (MCEq) Monte Carlo software [24,25] to compute the fluxes for the parent mesons and protons in the atmosphere, with the SYBILL-2.3 hadronic interaction model [26], the Hillas-Gaisser cosmic-ray model [27] and the NRLMSISE-00 atmospheric model [28]. With this procedure, meson and proton fluxes are obtained as a function of X, E and cos θ (see Appendix B for details). At this point, it is worth mentioning that these fluxes are subject to significant theoretical uncertainties related to the choice of cosmic-ray and hadronic interaction models, see e.g. Fig. 9 in Ref. [25] or Ref. [29], which can be of the order of 20-25%. Additional uncertainties may arise from seasonal and geomagnetic effects. As discussed in Ref. [30], while seasonal effects are negligible at SK, they could be potentially large for the fluxes observed at Icecube. However, they tend to cancel out when you take data over a whole year (the impact of the seasonal effect for Icecube can be seen e.g. in Ref. [31]). On the other hand, geomagnetic effects mostly affect the azimuthal distribution and horizontal events [32]. Since we only study time-averaged zenith distributions, their effects in the azimuthal distributions are averaged out. For SK, the three-dimensional effects change the height distribution of atmospheric neutrinos predominantly in the horizontal direction rather than the vertical direction where our signal dominates; see Figs. 18 and 19 in Ref. [33]. In the case of the IceCube analysis effects on the horizontal events are negligible as geomatic effects are only relevant for neutrino energies below 20 GeV (see Fig. 11 in Ref. [33]). We will assume that the boost is large enough E m P so that the zenith angle of the LLP is that of the parent particle and the integration over the parent zenith angle becomes trivial. In the following we will consider models where parent particles of LLP are protons and mesons. The differential fluxes of these particles as a function of the energy are shown in Fig. 1 at a fixed height, and as a function of the height at fixed energy in Fig. 2.
In the rest of this section we discuss separately the different production mechanisms of LLPs, from SM mesons and τ lepton decays, as well as from proton bremsstrahlung. In the case of production via meson decays, we will further distinguish between two-and three-body decays, as the outgoing LLP energy distributions will be different in the two cases.

Production from decays of SM particles
Let us first consider that a LLP, A, is produced in the decay of a SM meson P . Once the flux of the parent meson has been computed, the production profile of the LLP in the decay P → AY or P → AY Z is given by [23] where the sum goes over all channels (labeled as ch) contributing to the LLP production from P decays. Here, dn ch /dE is the fraction of decaying parents that produce an A with energy in the bin [E, E + dE] in the decay channel ch, while γ P , β P and τ P are the boost, velocity and lifetime in the rest frame of the parent particle P . Both the distribution dn ch /dE and the integration limits for E P (E min P , E max P ) depend on the decay kinematics and, in particular, they depend on whether the LLP is produced in a two-or in a three-body decay.

Two-body decays
In two-body decays P → AY , the differential distribution dn/dE is flat in E: The kinematical limits for the energy are: where E max is the energy of the LLP A when the decay takes place in the meson's rest frame, that is: (2.8) and γ P = E P /m P . Alternatively if E is fixed, the kinematical limits for E P m P are: . (2.9)

Three-body decays
In the case of three-body decay, P → AY Z, we have: The distribution is no longer flat and the kinematical limits, neglecting masses of Y and Z, are: Fixing E, the kinematical limits for E P are therefore: (2.13)

Bremsstrahlung
An alternative production mechanism for vectors and scalars is through bremsstrahlung, in the scattering of SM particles with air nuclei. Let us consider the case of a dark photon V of mass m V emitted by a proton p. In analogy to Eq. (2.3), the production profile of the dark photon can be written in this case as [23]: where Φ p is the proton flux and λ p is the interaction thickness of protons in air, i.e., the average amount of atmosphere (in g/cm 2 ) traversed between successive collisions with air nuclei N . It is common to parametrize the coupling of the dark photon to the proton as the QED coupling multiplied by a small number . In this case, the differential distribution dn/dE takes the form [34][35][36]: where |p max ⊥ | = 1GeV and Here, m p is the mass of the proton and with z ≡ p V /p p , where p V and p p are the dark photon and proton momenta, respectively.

Expected number of LLP decays in atmospheric neutrino detectors
Once the production profile of the LLP has been computed, the flux of particles that arrive at the detector Φ A is simply obtained integrating over all LLPs produced at different distances from the detector, weighted by their corresponding survival probabilities, as where dΠ A is the flux of LLP produced at a distance of the detector within the interval [ , + d ], computed as outlined in Sec. 2, and decay is the decay length of the LLP in the lab frame, related to its lifetime (τ A ) as where β A , γ A are the boost parameters of the LLP, and we assume E A m A . We will assume that, once the particle is produced, since it is very weakly coupled it does not interact any further in the atmosphere and only decays. Also, note that the upper limit of the integral in Eq. (3.1) is in practice a function of θ (see Eq. (2.2)): max ≡ (h max , θ) and h max 80km is the maximum height of the atmosphere where cosmic showers are produced.
Using the flux in Eq. (3.1), the number of decays inside the detector within a given time window ∆T , for LLPs with energies and trajectories in the intervals [E A , E A + dE A ], [cos θ, cos θ + d cos θ], can be computed as where A eff is an effective area which accounts for the probability that a decay takes place inside the detector. This area can be estimated integrating the surface of the detector normal to the flux direction, weighted by the decay probability of the LLP inside the detector: .

(3.4)
Here, ∆ det is the length of the segment of the LLP trajectory that cuts into the detector. Its calculation, which just depends on the zenith angle defining the trajectory of the LLP and on the geometry of the detector, is outlined in App. A. Figure 3 shows the effective decay areas for the IceCube and Super-Kamiokande detectors as a function of the lifetime of the LLP in the lab frame.

Detector efficiencies and datasets used
The number of events will also depend on the detector efficiencies and reconstruction effects. In this work, two atmospheric neutrino detectors have been considered: IceCube (IC) [37] and Super-Kamiokande (SK) [38]. Since they observe events in two very different energy regimes, we expect their results to be complementary and to probe different regions in parameter space. In this section we describe the assumptions and data sets used for each of the two detectors separately.

IceCube
The IceCube Neutrino Observatory [37] is a ∼ 1 km 3 neutrino detector at the Geographic South Pole, optimized for detecting neutrinos with energies above 100 GeV. For the effective area, Eq. (A.1), we consider a simplified cylindrical geometry with H = 1 km and volume 1 km 3 . The effective area as a function of the decay length in the laboratory frame is shown in Fig. 3.
We consider the data sample corresponding to the analysis presented in Ref. [39], for which effective areas for interacting neutrino events are publicly available in the IceCube collaboration webpage [40]. The statistics corresponds to 641 days of contained events, and the reconstructed data samples are divided into tracks and cascade events, in the northern and southern hemispheres, and are binned in reconstructed energy. Effective areas for each of these samples are provided for each particle, interaction type and true neutrino energy. From these effective areas we can infer the detector efficiencies by dividing by the corresponding interaction cross sections. We therefore estimate the detector efficiencies as where α selects the reconstructed sample: cascade/track and north/south, β selects the true neutrino interaction type, CC-e, NC, etc. In the denominator we have the total cross section in the IceCube detector σ ν β , which is estimated in terms of the nucleon-neutrino cross section, times the number of nucleons in the detector: Avogadro number N A , times the ice density (ρ ice 0.92 gr cm −3 ) times the fidutial volume V IC . The required neutrinonucleon cross sections are extracted from Ref. [41].
We will make the assumption that our LLP decay events have the same efficiencies as charged-current (CC) neutrino interactions with a similar set of observed particles in the final state. We also assume that the relationship between the energy of the decaying particle and the deposited energy of all particles produced in the detector is similar to the relationship between true and reconstructed neutrino energies in CC interactions. We believe that this is a reasonable as long as all particles produced in the LLP decay are visible in the detector. More concretely we will consider only LLP visible decays to electrons and charged hadrons, mainly because the expected neutrino background is smaller, they are more likely to be contained, and their the energy can be better determined. We assume that LLP decays of this type have an efficiency similar to those of CC e-like neutrino interactions: In our IC calculations, we use the efficiencies corresponding to α =cascades. Therefore, we weight the number of signal events by the branching ratio into e-like events, in order to consider all decay channels which have electrons, photons, τ -leptons or hadronic showers in the final state, which would be observed as cascade-like events. An important point is that these effective areas do not include the muon veto applied in the event pre-selection. This veto reduces events with a muon signal within a 3µs time window of the neutrino signal. This veto is designed to optimize the search for astrophysical neutrinos, as it reduces significantly the overwhelming background from atmospheric muons, and the atmospheric neutrino background in southern sky events [42]. At sufficiently high energies it is quite likely that muons from the same atmospheric shower that produce the LLPs will also reach the detector within this large time window, and therefore it also reduces significantly our signal. Since this veto has been applied to the background and the data, we need to impose it also on our signal to extract bounds, however it is quite likely that an optimized LLP search can be carried out, where the time window for a coincident muon track is reduced significantly, since the LLP has an atmospheric origin and an accompanying muon is not unexpected. The veto is expected to affect similarly our signal and the prompt neutrino signal, which is the dominant background, so it does not improve the signal-to-noise ratio of the LLP search. To implement this veto on our signal we use the passing fractions of the muon veto, P pass (E A , θ) obtained in Ref. [43], which depend on the true energy and zenith of the LLP. Given that we are concentrating on cascade events, and in most LLP production there are no associated muons, the most appropriate passing fractions to use are those corresponding prompt ν e .
The number of signal events in the i-th energy bin at IC can be finally obtained as where 0.2 ≤ cos θ rec ≤ 1 corresponds to the southern sky, while α = cascade and β =CC e-like. It might be useful to also include DeepCore atmospheric data. However, the detector volume in this case is much smaller than for IceCube, and the energy range considered for atmospheric neutrino analyses for DeepCore is comparable to that of the multi-GeV data in SK. Therefore, we expect similar results for DeepCore as those obtained for SK. It may however be interesting to search for a signal using atmospheric neutrinos with energies between 100 GeV and 1 TeV, which may be sensitive to a different range of values of cτ . This is unfortunately not possible using currently available public data.

Super-Kamiokande
Super-Kamiokande is a significantly smaller water Cherenkov detector, but it has a lower energy threshold where the atmospheric flux is much higher and therefore might have a comparable or event better sensitivity than IceCube to the scenarios considered. In the case of SK, in the computation of the effective decay areas (Eq. (A.1)) we have assumed a cylindrical geometry with H = 0.04 km and R = 0.02 km. The result is shown in Fig. 3.
We will again only consider electron-like events, that is, events with an e-like topology as defined above. We have used in this case the fully-contained multi GeV e-like sample from Ref. [44] (including single and multi-ring) with reconstructed zenith angle 0 ≤ cos θ ≤ 1. We extract data as well as neutrino background from Fig. 5 of [44]. In this case, while we do have information regarding the zenith angle of the event sample (see Fig. 5 in Ref. [44]) we do not have available information regarding the energy of the events. Therefore, for SK we bin the data only in cos θ and integrate over all neutrino energies between 1 GeV and 90 GeV, which is the range of parent neutrino energies for multi-GeV e-like contained events (see Fig. 6 of Ref. [44]). We think this is conservative, as SK may be sensitive to events outside this range. We also assume that the efficiencies for the decay of the LLP are similar to those of electron neutrino CC events in the multi-GeV range. Our only assumption regarding detection efficiencies for Super-Kamiokande is that they are flat as a function of energy, which we think is a reasonable assumption for fully contained events. We note that migration matrices for SK are not publicly available.
The number of events in the i-th bin in cos θ can therefore be computed as: where cos θ min i and cos θ max i are the lower and upper limits of the bin. In this case, we have assumed a flat detection efficiency SK = 0.75, in line with the values quoted in Ref. [44] for the multi-GeV ν e event sample.

Atmospheric LLP in selected scenarios
We now consider three simple BSM scenarios of very weakly interacting sectors that can lead to LLPs. In the first two examples, both production and decay are controlled by the same couplings, in such a way that the requirements of long enough lifetimes and large enough production go in opposite directions. Instead in the third example we consider a scenario where it is possible to decouple production from decay.

Heavy Neutral Leptons
The existence of heavy neutral leptons (HNL) is a generic prediction of Type-I seesaw models of neutrino masses [5][6][7][8]. The model is the simple extension of the SM with at least two heavy Majorana singlets N j . In the basis where the Majorana mass terms are diagonal, the Lagrangian reads: whereΦ ≡ iσ 2 Φ * is the complex conjugate of the Higgs field Φ, L α is the SM lepton doublet with flavor α, Y j is a Yukawa coupling and m N j is the Majorana mass of the singlet N j , which in principle is a free parameter of the model. After spontaneous symmetry breaking, the heavy Majorana states mix with the standard neutrinos resulting in a spectrum of almost standard light states (with masses m ν ∝ (Y v) 2 m N ) which correspond to the light neutrinos, and almost singlet heavy states corresponding to the HNL (with masses ∝ m N ). In the following, the full leptonic mixing matrix will be denoted as U and the mixing between the charged leptons and the heavy states is given by |U αj | 2 , which naively scales . It is through this mixing that the heavy singlets could be produced either through CC or neutral-current (NC) processes and also how they would decay back to SM particles.
A priori, the Majorana masses of the HNL are free parameters of the theory and therefore have to be probed experimentally. The GeV range is interesting from the theoretical point of view, since it has been shown that models with HNL at the GeV scale could explain neutrino masses and mixing parameters, as well as the matter-antimatter asymmetry of the Universe [9,10] without conflicting cosmological bounds. The search for GeV scale HNL is a very active field and several future experiments such as SHIP [45], FCCee [46], or  DUNE [47] have the potential to improve significantly over present bounds [13,14]. On the phenomenological front, new avenues to constrain these models are being proposed and/or further explored using current and past data [15,16,[48][49][50][51][52].
Here we consider instead the production of HNL in atmospheric showers, where their dominant production mechanism is through meson decays. Depending on their mass, the leading production channel is π, K, D (s) -meson decays (or in the decays of even heavier resonances, such as B mesons). Particularly interesting is the mass range slightly above the kaon mass, where existing laboratory bounds are weaker [13][14][15][16]. We have therefore considered the production via D-and D s -meson decays, as well as from τ decays if one considers the possibility that the production is dominated by the poorly constrained coupling to taus. The production in the lower-mass range from kaon and pion decays will be considered elsewhere. For simplicity, hereafter we use the notation U α ≡ U αj (α = e, µ, τ ) to refer to the elements of the leptonic mixing matrix that control both the production and decay of the heavy state N .
For D (s) meson decays, and for HNL between 0.5 and 1.5 GeV the dominant decay is two-body: D ± → N l ± α , where the flavour of the lepton α = e, µ (note that, if m N ≥ m K , this decay with α = τ is kinematically forbidden). The N fluxes are obtained from Eqs. (2.3), (2.4) and (2.5). The production profile of N from D and D s decays are shown in Fig. 4 in units of the corresponding branching ratios into this channel.
In the case of taus, the decays to produce N are necessarily three-body. Assuming that the dominant mixing is to τ leptons and the mass of the HNL is in the GeV range, the dominant decay channels are τ ± → N ν α l ± α , where l α ≡ e, µ. In this case it is more convenient to perform the change of variables E → z ≡ E/γ τ , where γ τ = E τ /m τ is the boost factor of the τ . Equation (2.10) then becomes The differential decay width in the rest frame of the τ lepton can be found in Ref. [53], or in Appendix B of Ref. [54]. Neglecting all lepton masses except those of the τ and N , after boosting to the laboratory frame we get 1 The contribution to the flux from τ decays is also shown in Fig. 4 in units of the τ → N l α ν α branching ratio. The partial decay widths of D, D s and τ to N can be obtained from Refs. [53], [54], and [55]. For the D and D s mesons, the total branching ratio for two-body decays with a HNL in the final state is obtained adding the contributions from D (s) → N e and D (s) → N µ. Assuming that the partial decay width of the parent meson to N is very small compared to its SM decay width, the branching ratio reads (4.5) where Γ D (s) is the total decay width of the D (s) meson (taken from Ref. [56]), y i ≡ m i /m D (s) and λ is defined in Eq. (2.6). In Eq. (4.5) G F is the Fermi constant, f D (s) is the decay constant of the parent meson and V D (s) is the mixing matrix element in the CKM matrix that participates in the decay vertex. The values used for the meson masses and decay constants (f D = 212 MeV; f Ds = 249 MeV) have been taken from Ref. [56]. For the CKM matrix elements, we use V D ≡ V cd = 0.22 and V Ds ≡ V cs = 0.995. In the case of decays from taus, we assume that the τ branching ratio into N is dominated by the value of U τ since it is less constrained by other experiments. In this case, we get: where, again in this case, for the total value of the width of the τ we use the SM value of Γ τ (i.e., ignoring the small contribution of the HNL). The expected number of decays inside the volume of the detector will also depend on the lifetime of the N , as outlined in Sec. 3. The total width is computed adding the partial widths into two-body decay channels (that is, into charged mesons and charged leptons, or neutral mesons and neutrinos), the three-body decay channels into charged leptons and neutrinos, and the three-body invisible decay into three light neutrinos. The partial widths for these decays can be found in several references [13,47,54,57] although they do not all agree. We have re-computed the two-body decay widths into mesons and leptons in the effective Lagrangian at low energies, and found good agreement with the results from the recent review in Ref. [57]. The total width can be written as where the dependence on the mixing matrix elements goes as Γ α ∝ |U α | 2 . In Fig. 5 we show the ratios |U α | −2 Γ α as a function of the N mass. From these, it is easy to reconstruct the total width and lifetime for any given values of the mixing matrix elements |U α | 2 . As already outlined, in this model both the production and decay rates (that is, the lifetime of the LLP) are controlled by the same set of parameters and therefore are highly correlated. Figures 6 show the regions allowed on the Br(P → N ) vs cτ N plane (for P = D, D s , τ ) allowing the mixing matrix elements to vary within the presently allowed regions [13], for two values of the HNL mass. As can be seen from both figures, in this model it is possible to achieve values of the lifetime in the right ballpark needed to obtain a signal in neutrino detectors (cτ ∼ O(10) km, after boosting to the lab frame), although at the price of very small production branching ratios. Br(τ ± → N l ± α ν α ) (right panel) versus cτ N , for m N = 0.5 GeV (dashed line, blue region) and 1 GeV (solid line, red region). The lines correspond to 10 −10 ≤ |U e | 2 = |U µ | 2 = |U τ | 2 ≤ 10 −6 (10 −7 ) for the lighter (heavier) mass. The ranges correspond to 10 −10 ≤ |U e | 2 = |U µ | 2 ≤ 10 −6 (10 −7 ) and 10 −10 ≤ |U τ | 2 ≤ 10 −2 (10 −3 ).
Finally, the number of events will also be proportional to the branching ratio of the decays leading to a cascade-like signature (at IC), or leading to e-like events (at SK). In practice, this amounts to adding the branching ratios into decay channels with either electrons, taus or hadronic resonances in the final state: where the processes involving light neutrinos include the diagrams for all active neutrino flavors, and the sum over M includes all charged or neutral pseudoscalar and vector mesons below the mass of the N (π ± , K ± , ρ ± , π 0 , η, η , ρ 0 , etc). The result is shown in Fig. 7 as a function of the mass of the HNL.

Dark Photons
Extensions of the Standard Model with an extra secluded U (1) gauge boson are ubiquitous in the BSM landscape. We first consider a model where the dark photon V is coupled to the visible sector via kinetic mixing [12]: where B µν is the hypercharge field strength tensor and V µν is analogously defined for the dark U (1), which is assumed spontaneously broken. The Lagrangian in Eq. (4.9) implies that the dark photon couples universally to all charged particles, like the SM photon, although with a coupling that is reduced by the factor . An enormous amount of work has been done in recent years to derive bounds on this scenario from beam dump experiments, e + e − and pp colliders, neutrino scattering and other intensity frontier experiments. A recent summary of present bounds can be found in [62].
The mechanisms leading to the production of dark photons in the atmosphere are the same as in a proton beam dump, which has been extensively studied in the literature, see e.g., Refs. [34,63,64]. For m V ≤ m π 0 (m η ), the dominant production channel is the two-body decay π 0 (η ) → γV . For heavier masses and m V ≤ m p , it can be produced via bremsstrahlung. The flux of V from neutral meson decays can be obtained from Eqs. (2.3), (2.4) and (2.5), with the production branching ratios [65]: for P = π 0 , η. The flux from bremsstrahlung is obtained from Eqs. (2.14) and (2.15). The production profile of dark photons from these processes in units of −2 are shown in Fig. 8.
The next ingredient we need in order to compute the expected number of decays in the detector is the lifetime of the dark photon. The partial decay widths into leptonic and hadronic channels can be found, for example, in Refs. [34,63,64]: Production profile of dark photons produced from π 0 , η decay and proton bremsstrahlung, for m V = 10 −2 GeV. The profile is shown at approximately 15 km from the Earth's surface, since this is where the parent meson fluxes are expected to be maximal, see Fig. 2, and for cos θ = 0.9 since we expect the signal to be largest as cos θ → 1.
where l stands for a SM charged lepton and α QED ≡ e 2 /(4π) stands for the SM fine structure constant. Adding up all possible decay channels below the hadronic threshold, and neglecting the electron mass, we get: (4.13) Conversely, above the hadronic production threshold we can approximate the lifetime by the fitted formula: Again in this case, both the production rate and the decay length depend on just two parameters: and m V . Figure 9 shows the correlation of and the lifetime of the dark photon, for three different values of m V in the range considered in this work. Finally, in our sensitivity calculations we focus on decays of the LLP leading to an electron-like signal in the detectors (i.e., decays with electrons or hadrons in the final state). Therefore, in order to compute the number of events in this sample we need the corresponding branching ratio. The result is shown in Fig. 10 as a function of the mass of the dark photon.

Heavy Neutral Leptons and U (1) B−L gauge symmetry
We have finally considered an extension of the SM that includes both HNL and a new U (1) gauge symmetry associated to the B − L number. Models with HNL and extra U (1) have been recently discussed in the context of the MiniBoone anomaly [51,[66][67][68], as well as in the context of explaining the matter-antimatter asymmetry or/and the observed dark matter abundance in the Universe [69]. The relevant terms in the Lagrangian are where g B−L is the gauge coupling of the dark photon, and the sum runs over all standard model fermions, which are charged under this new symmetry (Q f B−L = 1/3 for quarks, and Q f B−L = −1 for leptons). In this model, the origin of the B-L gauge boson mass typically requires a more complex hidden sector. However, our approach is phenomenological and therefore we will not address the possible origin of the the mass term which this is beyond the scope of this work. Finally, note that in principle a kinetic mixing term could also be added to Eq. (4.15). If the kinetic mixing is large, the phenomenological consequences of the model reduce to the dark photon case, which was already discussed in Sec. 4.2. Here we focus on the opposite limit instead, when the kinetic mixing term is negligible and the coupling dominates, since it leads to a different phenomenology. While a kinetic mixing term is always generated at loop level, this is expected to be subleading.
An interesting feature of this model is that the production of the N is dominated by the decays of the dark photon V → N N , provided m N ≤ m V /2, while its decay is controlled by the mixing with the light neutrinos. We might expect in this case to have an enhanced production if g B−L is not too small, while the N might be very long-lived as shown in Fig. 6. Present upper limits [62,70,71] on g 2 B−L are around ∼ 10 −7 in the mass range of interest, around m V ∼ O(GeV). From Fig. 6 we see that cτ N in the range [10 −2 , 10 5 ] km for m V ∼ 0.5 GeV are allowed.
Neglecting the decay of the dark photon, the flux of V can be obtained as in Sec. 4.2 simply replacing [62]: However, in the B − L model the dark photon decays very promptly. Its lifetime (for 2m N ≤ m V , and 2m µ ≤ m V ≤ m ω ) is approximately given by (4.17) Therefore, decay effects should be accounted for in the calculation of the flux, following the same arguments as in App. B. It is easy to show that, if we denote byφ V the dark photon flux obtained without accounting for the decay, then the solution to the cascade equation including decay effects is given by: The flux of dark photons taking into account the decay is shown in Fig. 11 for a three different values of g B−L in the experimentally allowed range. Interestingly, while at high energies there is clearly a suppression of the dark photon flux, which depends on the value of g B−L , at low energies the production remains practically independent from the value of g B−L . This can be understood as follows. In this part of the spectrum the lifetime of the dark photon is very short and the exponential in Eq. (4.18) goes to zero very rapidly, except for values of such that ∆ = − ∝ cτ V . Thus, at low energies the resulting flux will take the form (4.19) and since cτ V ∝ g −2 B−L this effectively renders the result independent of g B−L at low energies. From these fluxes, we can get the production profile of N as in Eq. (2.3), where now the parent particle P is the dark photon instead of SM mesons. The decay V → N N is two-body and therefore the decay distribution is the same as in Eqs. (2.4) and (2.5). Neglecting the electron mass, its partial decay width reads: Finally, the branching ratio in this channel (in the mass region of interest) can be well approximated by: since the decays into mesons are suppressed. The resulting production profile of N from dark photon decays is shown in Fig. 12.
Since the production and decay processes are effectively decoupled, the decay length of the N (which is controlled by the mixing) can now take a wide range of values (as shown in Fig. 6) without affecting the production rates.

Results
In this section we present our numerical results for the sensitivity to the scenarios described in Sec. 4. Using the data sets described in Sec. 3, a χ 2 fit to the data is performed. A Poissonian χ 2 function has been used, defined as: where the sum runs over the energy (angular) bins for the IC (SK) detector analyses, and we use α = cascade(e-like) for IC (SK). Here, n α i stands for the data observed in each bin while N α i is the predicted number of signal events and B α i is the background prediction, which includes the predicted number of atmospheric neutrino events in the SM (plus the contribution from the fitted astrophysical neutrino flux, in the case of IC).
A comment regarding the treatment of the astrophysical neutrino background for IC is in order, since the distribution of astrophysical neutrinos is obtained from a fit to the same data that we use in order to derive a limit to LLP models. Although this may seem inconsistent, it should be noted that due to the nature of the decay and the spectra of the parent mesons in the atmosphere (which obey a power law), at IC we expect most of our sensitivity to come from events observed at energies below 10 TeV, where the contribution  Figure 13. Expected number of decays for the IC Medium Energy Starting Events (MESE) sample, as a function of the deposited energy in the detector and for the down-going sample only (cos θ > 0), for two of the scenarios considered in this work. The shaded histograms show the background prediction, including a fitted distribution to the astrophysical neutrino signal, the cosmic muon background, and atmospheric neutrino events. The black dots correspond to the observed data with error bars, as in Ref. [39]. In the left panel, the signal has been computed for a HNL with a mass m N = 0.6 GeV, a cτ N = 10 −4 km and Br(D s → N ) Br(N → CC e-like) = 7.4 · 10 −3 . In the right panel, the signal has been computed for a dark photon with a mass m V = 25 MeV for = 2 · 10 −4 and an (uncorrelated) lifetime cτ V = 5 · 10 −6 km. from the astrophysical neutrinos is subdominant. Moreover, from the comparison between the left and right panels in Fig. 8 in Ref. [39], it seems that the fit to the astrophysical neutrino background is mostly driven by the data from the northern sky, whereas we expect LLP decays to contribute mostly to the southern sky data set, which is the only data we use. Therefore, although a more detailed analysis by the experimental collaboration would be necessary to do this analysis properly, we believe the results would not be very different. Of course, given the very characteristic zenith angle dependence of the LLP signal, a greater sensitivity is expected if 2D binned data (using both energy and angular information) were to be used instead.
For illustration, Fig. 13 shows the expected number of signal events for HNL (left panel) and dark photons (right panel), for an assumed value of the production branching ratio of the HNL and its lifetime. The shaded histograms show the background prediction, while the solid lines show the expected signal plus background event rates per bin. For reference, the black points show the observed data, as given in Ref. [39].
The zenith angle distribution for the neutrino background and LLP is signal events in Super-Kamiokande is shown in Fig. 14. As expected the signal is peaked at small zenith angles, that is, most of the decaying LLP are expected to come from trajectories entering the detector from above.
Before moving on to the discussion of our results, we should note that no systematic uncertainties have been included in our calculations, since this is an extremely challenging task using only publicly available information. Eventually, a detailed study including a proper implementation of systematic uncertainties and detector reconstruction effects should be performed by the experimental collaborations themselves.

Exclusion limits for IceCube and Super-Kamiokande
Following the procedure described above, we proceed to derive exclusion limits for the HNL and dark photon LLPs, for both IC and SK. In doing so, we consider that the data is chi-squared distributed with n degrees of freedom (d.o.f.), n being the number of bins in the data (that is, 5 bins in cos θ in the case of SK, as opposed to 15 energy bins in the case of IC). Our exclusion limits thus indicate the region in parameter space where the χ 2 value exceeds the corresponding one at 90% CL, regardless of where the best-fit point lies in the parameter space. Our results are shown in Fig. 15 for the HNL, assuming m N = 1 GeV, and in Fig In the case of the HNL (Fig 15), the solid red, dashed blue, and dot-dashed green lines show our results assuming that the LLP is mainly produced from D, D s or τ decays. Although our results are shown for m N = 1GeV, the regions do not change significantly for other values of m N between 0.5 and 1.5 GeV. As can be seen from this figure, the limits obtained for the two experiments are very similar, although slightly better for SK. However, in both cases the limits cannot probe the allowed regions of parameter space once we take into account the correlation between production and decay, as shown in Figs. 6. Therefore we conclude that, in the minimal scenario described in Sec. 4.1, these limits are not competitive (although this may not be the case in non-minimal BSM scenarios where the production and decay may be uncorrelated).
The exclusion limits for the dark photon scenario (Sec. 4.2) are shown in Fig. 16. In this case our results include all production mechanisms kinematically available for the production of the dark photon: π 0 decays, η decays, and p bremsstrahlung. The solid red, dashed green and dotted blue lines indicate the results obtained for three different values of the dark photon mass, as indicated in the legend. Again in this case, as for the HNL scenario, we find that this search cannot probe the region of parameter space shown in Fig. 9 for correlated production and decay rates. However, the limits might be useful in the context of more complex models where production and decay are uncorrelated. In particular, it is noticeable that the limits for SK are three orders of magnitude better than for IC. This is because, in the case of light mesons and protons, the flux follows a harder power law than for D (s) mesons and taus, and therefore reducing the energy of the events detected leads to huge enhancement in the expected number of events.  On the other hand, in the B − L model it is possible to uncorrelate production and decay since the LLP is not the dark photon, but the HNL to which it decays. Therefore, production rates in this model are controlled by the B − L gauge coupling, while the decay depends on the mixing with the light neutrinos. The limits for the B − L model

IceCube preferred regions
An interesting observation is that, in the case of IC, the minimization of the χ 2 on these two-dimensional planes gives a minimum of the χ 2 that lies at a non-zero value of the LLP production rate. This is the result from a small excess in the data at around 30 TeV, as can be seen by naked eye in Fig. 13. Although the source of the excess is unclear, and it might be related to an underestimation of the µ background or an unidentified source of systematic uncertainties, in the remainder of this section we explore possible explanations in the LLP models discussed in Sec. 4. Figure 18 shows the 1σ regions for 2 d.o.f., i.e., the region ∆χ 2 = χ 2 − χ 2 min ≤ 2.3, where χ 2 min is the global minimum, for a HNL produced from D, D s and τ decays. On the other hand, χ 2 (no LLP) -χ 2 and 4.2 cannot live on those regions, according to Figs. 6 (for the HNL) and Fig. 9 (for the dark photon).
Regarding the B − L model with HNL, as in the previous cases we get a best fit to IC data away from zero and, in particular, for Br(V → N )g 2 B−L ∼ 10 −6 − 10 −5 and cτ N ∼ 10 −4 km. Unfortunately, the required value of g B−L is excluded by collider constraints [62,70,71] in this dark photon mass range and, at the same time, such as small cτ N for a HNL at m N = 0.35 GeV is at least one order of magnitude too small, assuming that the mixing is below the present constraints from direct searches of HNL [13].
Before finalizing this section it should be stressed, however, that our IC analysis has been performed using just the energy information available on the data release of Ref. [39], and that in the case of decay of an LLP the zenith angle distribution of the signal events could differ significantly from the distributions observed for the data (see Fig. 9 of [39]). More stringent constraints might be attainable if information on the zenith angle distribution of the events were included in the analysis, which is not possible with the publicly available information.

Conclusions
Very weakly interacting particles beyond the Standard Model (SM) might be light enough to be produced in accelerators, but easy to miss in standard BSM search analyses if they are very long-lived. Many running accelerator experiments, both fixed-target and colliders, are actively exploring new strategies to improve the sensitivity to such exotic signals, and several future experiments are being proposed to target specifically such BSM scenarios.
In this paper we have presented the first detailed analysis of the sensitivity of large neutrino detectors, such as Super-Kamiokande and IceCube, to putative long-lived particles (LLP) that might be produced in atmospheric showers from the decay of standard model mesons and/or bremsstrahlung. We have presented the methodology to evaluate LLP fluxes from these processes, and the procedure to quantify the expected number of signal events. We have considered in particular two minimal and theoretically well-motivated scenarios where the LLP are either heavy neutral leptons produced primarily in meson decays, or dark photons produced from bremsstrahlung or/and radiative π 0 and η decays. We have evaluated bounds using Super-Kamiokande and IceCube data from Refs. [39,44], as a function of the branching ratios (or production rates) of SM → LLPs, and the LLP lifetime. The best sensitivity is obtained for lifetimes of O(10 −2 × m LLP 1GeV ) km in Icecube and O(10 × m LLP 1GeV ) km for Super-Kamiokande. However, in the minimal models considered in this work production and decay are usually strongly correlated since both are controlled by the same coupling. In this case, we find that atmospheric searches do not lead to competitive bounds. On the other hand, our bounds might be complementary to other searches in more complex models where production and decay are uncorrelated. As an example for such scenario we considered an extended model with a B − L gauge boson that couples to the HNL. In this case, the production is controlled by the gauge coupling, g B−L , while the LLP is the HNL whose lifetime is controlled by its mixing with the light neutrinos. We have also presented the sentitivity to this scenario.
Interestingly, in the case of Icecube the addition of an LLP to the SM provides a better fit to the observed data than the prediction obtained in the SM (with a ∆χ 2 ∼ 7), due to a small excess observed at around 30 TeV. This is obtained, for example, for a HNL produced in D meson decays with a Br(D → X) ∼ 2 · 10 −3 and a lifetime cτ ∼ 10 −4 . However, we find that the values of the parameters needed in order to fit the excess lie outside of the allowed regions of parameter space for the three models considered.
Finally, we stress again that no systematic uncertainties have been included in our calculations. A more detailed study, including a realistic implementation of systematic and reconstruction errors by the experimental collaborations, would be mandatory to validate our results.

A Computation of effective decay areas
Assuming the fidutial volume of the detector is a cylinder of radius R and height H, it is a simple geometrical exercise to obtain ∆ det , and the effective area: A eff decay (E A , cos θ) = | cos θ|A 1 (E A , cos θ) + | sin θ|A 2 (E A , cos θ), (A.1) where A 1 correspond to the flux entering from the top: and A 2 is that entering laterally:

B Computation of parent particle fluxes in the atmosphere
In all our calculations, the SM parent particle fluxes in the atmosphere have been computed using the MCEq software [24,25], with the SYBILL-2.3 hadronic interaction model [26], the Hillas-Gaisser cosmic-ray model [27] and the NRLMSISE-00 atmospheric model [28].
Obtaining the parent particle fluxes is however not straightforward, because only the flux for protons and unstable particles that are relatively long-lived compared to their interaction length can be directly extracted using MCEq. A possible strategy to circumvent this is to manually switch off the decay of the parent particle and extract the flux, which then needs to be corrected to account for the decay 2 .
To simplify notation, we denote the differential flux of the parent particle P by φ P : Let us assume thatφ P is the differential flux of parent particle P assuming that P is stable 3 . Assuming that P is directly produced in nucleon interactions, this flux satisfies a cascade equation of the form [23] where φ N is the flux of nucleons, and the spectrum-weighted momenta are defined as Here, λ k is the particle interaction length of hadron k and dn(kN →hY ) dE is the number of hadrons h produced with energy between E and E + dE, in the scattering of hadron k with energy E P on a nucleus N .
The second and third terms in Eq. B.2 are the production and regeneration terms, respectively, while the first term is the absorption term. Had the particle decayed, the correct cascade equation should have also included a decay term: where d P is the decay length of the parent particle in the laboratory frame and ρ is the column density at the point with slant depth X. It is easy to see that fromφ P we can get the solution to Eq. (B.4) as To the best of our knowledge, although MCEq allows to obtain the values of λ P it does not allow to extract the values of Z P P directly. Therefore, we assume Z τ τ = 0 while for the rest of the mesons considered we assume Feynmann scaling (that is, Z M M = 0.3) since, according to Fig. 4 in Ref. [23], this is a reasonable approximation for charged pions and kaons. In any case, we do not expect our results to change significantly if these assumptions are modified. With this procedure we get a reasonable agreement with the output of MCEq available at large energies, where decay can be ignored. We do a small rescaling of our approximate result by matching the two curves in the overlapping region.
Finally, in the case of D and D s mesons both the interaction and regeneration terms in Eq. B.4 can be neglected. In this case, the meson flux may be computed directly from the flux of protons in the atmosphere, following Ref. [23] (Eqs. (25)(26)(27) in that reference): where φ N (X) and the yields, Z N M , are extracted from MCEq. We have also checked that this method gives a good agreement with our fluxes, in the case of D and D s mesons.